vendor/gonum.org/v1/gonum/lapack/testlapack/dorml2.go

   1 // Copyright ©2015 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package testlapack
   6
   7 import (
   8         "testing"
   9
  10         "golang.org/x/exp/rand"
  11
  12         "gonum.org/v1/gonum/blas"
  13         "gonum.org/v1/gonum/blas/blas64"
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 type Dorml2er interface {
  18         Dgelqfer
  19         Dorml2(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64)
  20 }
  21
  22 func Dorml2Test(t *testing.T, impl Dorml2er) {
  23         rnd := rand.New(rand.NewSource(1))
  24         // TODO(btracey): This test is not complete, because it
  25         // doesn't test individual values of m, n, and k, instead only testing
  26         // a specific subset of possible k values.
  27         for _, side := range []blas.Side{blas.Left, blas.Right} {
  28                 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
  29                         for _, test := range []struct {
  30                                 common, adim, cdim, lda, ldc int
  31                         }{
  32                                 {3, 4, 5, 0, 0},
  33                                 {3, 5, 4, 0, 0},
  34                                 {4, 3, 5, 0, 0},
  35                                 {4, 5, 3, 0, 0},
  36                                 {5, 3, 4, 0, 0},
  37                                 {5, 4, 3, 0, 0},
  38
  39                                 {3, 4, 5, 6, 20},
  40                                 {3, 5, 4, 6, 20},
  41                                 {4, 3, 5, 6, 20},
  42                                 {4, 5, 3, 6, 20},
  43                                 {5, 3, 4, 6, 20},
  44                                 {5, 4, 3, 6, 20},
  45                                 {3, 4, 5, 20, 6},
  46                                 {3, 5, 4, 20, 6},
  47                                 {4, 3, 5, 20, 6},
  48                                 {4, 5, 3, 20, 6},
  49                                 {5, 3, 4, 20, 6},
  50                                 {5, 4, 3, 20, 6},
  51                         } {
  52                                 var ma, na, mc, nc int
  53                                 if side == blas.Left {
  54                                         ma = test.adim
  55                                         na = test.common
  56                                         mc = test.common
  57                                         nc = test.cdim
  58                                 } else {
  59                                         ma = test.adim
  60                                         na = test.common
  61                                         mc = test.cdim
  62                                         nc = test.common
  63                                 }
  64                                 // Generate a random matrix
  65                                 lda := test.lda
  66                                 if lda == 0 {
  67                                         lda = na
  68                                 }
  69                                 a := make([]float64, ma*lda)
  70                                 for i := range a {
  71                                         a[i] = rnd.Float64()
  72                                 }
  73                                 ldc := test.ldc
  74                                 if ldc == 0 {
  75                                         ldc = nc
  76                                 }
  77                                 // Compute random C matrix
  78                                 c := make([]float64, mc*ldc)
  79                                 for i := range c {
  80                                         c[i] = rnd.Float64()
  81                                 }
  82
  83                                 // Compute LQ
  84                                 k := min(ma, na)
  85                                 tau := make([]float64, k)
  86                                 work := make([]float64, 1)
  87                                 impl.Dgelqf(ma, na, a, lda, tau, work, -1)
  88                                 work = make([]float64, int(work[0]))
  89                                 impl.Dgelqf(ma, na, a, lda, tau, work, len(work))
  90
  91                                 // Build Q from result
  92                                 q := constructQ("LQ", ma, na, a, lda, tau)
  93
  94                                 cMat := blas64.General{
  95                                         Rows:   mc,
  96                                         Cols:   nc,
  97                                         Stride: ldc,
  98                                         Data:   make([]float64, len(c)),
  99                                 }
 100                                 copy(cMat.Data, c)
 101                                 cMatCopy := blas64.General{
 102                                         Rows:   cMat.Rows,
 103                                         Cols:   cMat.Cols,
 104                                         Stride: cMat.Stride,
 105                                         Data:   make([]float64, len(cMat.Data)),
 106                                 }
 107                                 copy(cMatCopy.Data, cMat.Data)
 108                                 switch {
 109                                 default:
 110                                         panic("bad test")
 111                                 case side == blas.Left && trans == blas.NoTrans:
 112                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 113                                 case side == blas.Left && trans == blas.Trans:
 114                                         blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 115                                 case side == blas.Right && trans == blas.NoTrans:
 116                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
 117                                 case side == blas.Right && trans == blas.Trans:
 118                                         blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
 119                                 }
 120                                 // Do Dorm2r ard compare
 121                                 if side == blas.Left {
 122                                         work = make([]float64, nc)
 123                                 } else {
 124                                         work = make([]float64, mc)
 125                                 }
 126                                 aCopy := make([]float64, len(a))
 127                                 copy(aCopy, a)
 128                                 tauCopy := make([]float64, len(tau))
 129                                 copy(tauCopy, tau)
 130                                 impl.Dorml2(side, trans, mc, nc, k, a, lda, tau, c, ldc, work)
 131                                 if !floats.Equal(a, aCopy) {
 132                                         t.Errorf("a changed in call")
 133                                 }
 134                                 if !floats.Equal(tau, tauCopy) {
 135                                         t.Errorf("tau changed in call")
 136                                 }
 137                                 if !floats.EqualApprox(cMat.Data, c, 1e-14) {
 138                                         isLeft := side == blas.Left
 139                                         isTrans := trans == blas.Trans
 140                                         t.Errorf("Multiplication mismatch. IsLeft = %v. IsTrans = %v", isLeft, isTrans)
 141                                 }
 142                         }
 143                 }
 144         }
 145 }