vendor/gonum.org/v1/gonum/lapack/testlapack/dormr2.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 Dormr2er interface {
  18         Dgerqf(m, n int, a []float64, lda int, tau, work []float64, lwork int)
  19         Dormr2(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 Dormr2Test(t *testing.T, impl Dormr2er) {
  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                                 ma := test.adim
  49                                 na := test.common
  50                                 var mc, nc int
  51                                 if side == blas.Left {
  52                                         mc = test.common
  53                                         nc = test.cdim
  54                                 } else {
  55                                         mc = test.cdim
  56                                         nc = test.common
  57                                 }
  58
  59                                 // Generate a random matrix
  60                                 lda := test.lda
  61                                 if lda == 0 {
  62                                         lda = na
  63                                 }
  64                                 a := make([]float64, ma*lda)
  65                                 for i := range a {
  66                                         a[i] = rnd.Float64()
  67                                 }
  68                                 ldc := test.ldc
  69                                 if ldc == 0 {
  70                                         ldc = nc
  71                                 }
  72                                 // Compute random C matrix
  73                                 c := make([]float64, mc*ldc)
  74                                 for i := range c {
  75                                         c[i] = rnd.Float64()
  76                                 }
  77
  78                                 // Compute RQ
  79                                 k := min(ma, na)
  80                                 tau := make([]float64, k)
  81                                 work := make([]float64, 1)
  82                                 impl.Dgerqf(ma, na, a, lda, tau, work, -1)
  83                                 work = make([]float64, int(work[0]))
  84                                 impl.Dgerqf(ma, na, a, lda, tau, work, len(work))
  85
  86                                 // Build Q from result
  87                                 q := constructQ("RQ", ma, na, a, lda, tau)
  88
  89                                 cMat := blas64.General{
  90                                         Rows:   mc,
  91                                         Cols:   nc,
  92                                         Stride: ldc,
  93                                         Data:   make([]float64, len(c)),
  94                                 }
  95                                 copy(cMat.Data, c)
  96                                 cMatCopy := blas64.General{
  97                                         Rows:   cMat.Rows,
  98                                         Cols:   cMat.Cols,
  99                                         Stride: cMat.Stride,
 100                                         Data:   make([]float64, len(cMat.Data)),
 101                                 }
 102                                 copy(cMatCopy.Data, cMat.Data)
 103                                 switch {
 104                                 default:
 105                                         panic("bad test")
 106                                 case side == blas.Left && trans == blas.NoTrans:
 107                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 108                                 case side == blas.Left && trans == blas.Trans:
 109                                         blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, cMatCopy, 0, cMat)
 110                                 case side == blas.Right && trans == blas.NoTrans:
 111                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, cMatCopy, q, 0, cMat)
 112                                 case side == blas.Right && trans == blas.Trans:
 113                                         blas64.Gemm(blas.NoTrans, blas.Trans, 1, cMatCopy, q, 0, cMat)
 114                                 }
 115                                 // Do Dorm2r ard compare
 116                                 if side == blas.Left {
 117                                         work = make([]float64, nc)
 118                                 } else {
 119                                         work = make([]float64, mc)
 120                                 }
 121                                 aCopy := make([]float64, len(a))
 122                                 copy(aCopy, a)
 123                                 tauCopy := make([]float64, len(tau))
 124                                 copy(tauCopy, tau)
 125                                 impl.Dormr2(side, trans, mc, nc, k, a[(ma-k)*lda:], lda, tau, c, ldc, work)
 126                                 if !floats.Equal(a, aCopy) {
 127                                         t.Errorf("a changed in call")
 128                                 }
 129                                 if !floats.Equal(tau, tauCopy) {
 130                                         t.Errorf("tau changed in call")
 131                                 }
 132                                 if !floats.EqualApprox(cMat.Data, c, 1e-14) {
 133                                         t.Errorf("Multiplication mismatch.\n Want %v \n got %v.", cMat.Data, c)
 134                                 }
 135                         }
 136                 }
 137         }
 138 }