vendor/gonum.org/v1/gonum/lapack/testlapack/dorm2r.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 Dorm2rer interface {
  18         Dgeqrfer
  19         Dorm2r(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 Dorm2rTest(t *testing.T, impl Dorm2rer) {
  23         rnd := rand.New(rand.NewSource(1))
  24         for _, side := range []blas.Side{blas.Left, blas.Right} {
  25                 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
  26                         for _, test := range []struct {
  27                                 common, adim, cdim, lda, ldc int
  28                         }{
  29                                 {3, 4, 5, 0, 0},
  30                                 {3, 5, 4, 0, 0},
  31                                 {4, 3, 5, 0, 0},
  32                                 {4, 5, 3, 0, 0},
  33                                 {5, 3, 4, 0, 0},
  34                                 {5, 4, 3, 0, 0},
  35                                 {3, 4, 5, 6, 20},
  36                                 {3, 5, 4, 6, 20},
  37                                 {4, 3, 5, 6, 20},
  38                                 {4, 5, 3, 6, 20},
  39                                 {5, 3, 4, 6, 20},
  40                                 {5, 4, 3, 6, 20},
  41                                 {3, 4, 5, 20, 6},
  42                                 {3, 5, 4, 20, 6},
  43                                 {4, 3, 5, 20, 6},
  44                                 {4, 5, 3, 20, 6},
  45                                 {5, 3, 4, 20, 6},
  46                                 {5, 4, 3, 20, 6},
  47                         } {
  48                                 var ma, na, mc, nc int
  49                                 if side == blas.Left {
  50                                         ma = test.common
  51                                         na = test.adim
  52                                         mc = test.common
  53                                         nc = test.cdim
  54                                 } else {
  55                                         ma = test.common
  56                                         na = test.adim
  57                                         mc = test.cdim
  58                                         nc = test.common
  59                                 }
  60
  61                                 // Generate a random matrix
  62                                 lda := test.lda
  63                                 if lda == 0 {
  64                                         lda = na
  65                                 }
  66                                 a := make([]float64, ma*lda)
  67                                 for i := range a {
  68                                         a[i] = rnd.Float64()
  69                                 }
  70                                 ldc := test.ldc
  71                                 if ldc == 0 {
  72                                         ldc = nc
  73                                 }
  74                                 // Compute random C matrix
  75                                 c := make([]float64, mc*ldc)
  76                                 for i := range c {
  77                                         c[i] = rnd.Float64()
  78                                 }
  79
  80                                 // Compute QR
  81                                 k := min(ma, na)
  82                                 tau := make([]float64, k)
  83                                 work := make([]float64, 1)
  84                                 impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
  85                                 work = make([]float64, int(work[0]))
  86                                 impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
  87
  88                                 // Build Q from result
  89                                 q := constructQ("QR", ma, na, a, lda, tau)
  90
  91                                 cMat := blas64.General{
  92                                         Rows:   mc,
  93                                         Cols:   nc,
  94                                         Stride: ldc,
  95                                         Data:   make([]float64, len(c)),
  96                                 }
  97                                 copy(cMat.Data, c)
  98                                 cMatCopy := blas64.General{
  99                                         Rows:   cMat.Rows,
 100                                         Cols:   cMat.Cols,
 101                                         Stride: cMat.Stride,
 102                                         Data:   make([]float64, len(cMat.Data)),
 103                                 }
 104                                 copy(cMatCopy.Data, cMat.Data)
 105                                 switch {
 106                                 default:
 107                                         panic("bad test")
 108                                 case side == blas.Left && trans == blas.NoTrans:
 109                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 110                                 case side == blas.Left && trans == blas.Trans:
 111                                         blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 112                                 case side == blas.Right && trans == blas.NoTrans:
 113                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
 114                                 case side == blas.Right && trans == blas.Trans:
 115                                         blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
 116                                 }
 117                                 // Do Dorm2r ard compare
 118                                 if side == blas.Left {
 119                                         work = make([]float64, nc)
 120                                 } else {
 121                                         work = make([]float64, mc)
 122                                 }
 123                                 aCopy := make([]float64, len(a))
 124                                 copy(aCopy, a)
 125                                 tauCopy := make([]float64, len(tau))
 126                                 copy(tauCopy, tau)
 127                                 impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, c, ldc, work)
 128                                 if !floats.Equal(a, aCopy) {
 129                                         t.Errorf("a changed in call")
 130                                 }
 131                                 if !floats.Equal(tau, tauCopy) {
 132                                         t.Errorf("tau changed in call")
 133                                 }
 134                                 if !floats.EqualApprox(cMat.Data, c, 1e-14) {
 135                                         t.Errorf("Multiplication mismatch.\n Want %v \n got %v.", cMat.Data, c)
 136                                 }
 137                         }
 138                 }
 139         }
 140 }