vendor/gonum.org/v1/gonum/lapack/testlapack/dlarf.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 Dlarfer interface {
  18         Dlarf(side blas.Side, m, n int, v []float64, incv int, tau float64, c []float64, ldc int, work []float64)
  19 }
  20
  21 func DlarfTest(t *testing.T, impl Dlarfer) {
  22         rnd := rand.New(rand.NewSource(1))
  23         for i, test := range []struct {
  24                 m, n, ldc    int
  25                 incv, lastv  int
  26                 lastr, lastc int
  27                 tau          float64
  28         }{
  29                 {
  30                         m:   3,
  31                         n:   2,
  32                         ldc: 2,
  33
  34                         incv:  4,
  35                         lastv: 1,
  36
  37                         lastr: 2,
  38                         lastc: 1,
  39
  40                         tau: 2,
  41                 },
  42                 {
  43                         m:   2,
  44                         n:   3,
  45                         ldc: 3,
  46
  47                         incv:  4,
  48                         lastv: 1,
  49
  50                         lastr: 1,
  51                         lastc: 2,
  52
  53                         tau: 2,
  54                 },
  55                 {
  56                         m:   2,
  57                         n:   3,
  58                         ldc: 3,
  59
  60                         incv:  4,
  61                         lastv: 1,
  62
  63                         lastr: 0,
  64                         lastc: 1,
  65
  66                         tau: 2,
  67                 },
  68                 {
  69                         m:   2,
  70                         n:   3,
  71                         ldc: 3,
  72
  73                         incv:  4,
  74                         lastv: 0,
  75
  76                         lastr: 0,
  77                         lastc: 1,
  78
  79                         tau: 2,
  80                 },
  81                 {
  82                         m:   10,
  83                         n:   10,
  84                         ldc: 10,
  85
  86                         incv:  4,
  87                         lastv: 6,
  88
  89                         lastr: 9,
  90                         lastc: 8,
  91
  92                         tau: 2,
  93                 },
  94         } {
  95                 // Construct a random matrix.
  96                 c := make([]float64, test.ldc*test.m)
  97                 for i := 0; i <= test.lastr; i++ {
  98                         for j := 0; j <= test.lastc; j++ {
  99                                 c[i*test.ldc+j] = rnd.Float64()
 100                         }
 101                 }
 102                 cCopy := make([]float64, len(c))
 103                 copy(cCopy, c)
 104                 cCopy2 := make([]float64, len(c))
 105                 copy(cCopy2, c)
 106
 107                 // Test with side right.
 108                 sz := max(test.m, test.n) // so v works for both right and left side.
 109                 v := make([]float64, test.incv*sz+1)
 110                 // Fill with nonzero entries up until lastv.
 111                 for i := 0; i <= test.lastv; i++ {
 112                         v[i*test.incv] = rnd.Float64()
 113                 }
 114                 // Construct h explicitly to compare.
 115                 h := make([]float64, test.n*test.n)
 116                 for i := 0; i < test.n; i++ {
 117                         h[i*test.n+i] = 1
 118                 }
 119                 hMat := blas64.General{
 120                         Rows:   test.n,
 121                         Cols:   test.n,
 122                         Stride: test.n,
 123                         Data:   h,
 124                 }
 125                 vVec := blas64.Vector{
 126                         Inc:  test.incv,
 127                         Data: v,
 128                 }
 129                 blas64.Ger(-test.tau, vVec, vVec, hMat)
 130
 131                 // Apply multiplication (2nd copy is to avoid aliasing).
 132                 cMat := blas64.General{
 133                         Rows:   test.m,
 134                         Cols:   test.n,
 135                         Stride: test.ldc,
 136                         Data:   cCopy,
 137                 }
 138                 cMat2 := blas64.General{
 139                         Rows:   test.m,
 140                         Cols:   test.n,
 141                         Stride: test.ldc,
 142                         Data:   cCopy2,
 143                 }
 144                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMat2, hMat, 0, cMat)
 145
 146                 // cMat now stores the true answer. Compare with the function call.
 147                 work := make([]float64, sz)
 148                 impl.Dlarf(blas.Right, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
 149                 if !floats.EqualApprox(c, cMat.Data, 1e-14) {
 150                         t.Errorf("Dlarf mismatch right, case %v. Want %v, got %v", i, cMat.Data, c)
 151                 }
 152
 153                 // Test on the left side.
 154                 copy(c, cCopy2)
 155                 copy(cCopy, c)
 156                 // Construct h.
 157                 h = make([]float64, test.m*test.m)
 158                 for i := 0; i < test.m; i++ {
 159                         h[i*test.m+i] = 1
 160                 }
 161                 hMat = blas64.General{
 162                         Rows:   test.m,
 163                         Cols:   test.m,
 164                         Stride: test.m,
 165                         Data:   h,
 166                 }
 167                 blas64.Ger(-test.tau, vVec, vVec, hMat)
 168                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, hMat, cMat2, 0, cMat)
 169                 impl.Dlarf(blas.Left, test.m, test.n, v, test.incv, test.tau, c, test.ldc, work)
 170                 if !floats.EqualApprox(c, cMat.Data, 1e-14) {
 171                         t.Errorf("Dlarf mismatch left, case %v. Want %v, got %v", i, cMat.Data, c)
 172                 }
 173         }
 174 }