vendor/gonum.org/v1/gonum/lapack/testlapack/dormqr.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         "fmt"
   9         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/blas"
  14         "gonum.org/v1/gonum/floats"
  15 )
  16
  17 type Dormqrer interface {
  18         Dorm2rer
  19         Dormqr(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
  20 }
  21
  22 func DormqrTest(t *testing.T, impl Dormqrer) {
  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                                 {6, 7, 8, 0, 0},
  30                                 {6, 8, 7, 0, 0},
  31                                 {7, 6, 8, 0, 0},
  32                                 {7, 8, 6, 0, 0},
  33                                 {8, 6, 7, 0, 0},
  34                                 {8, 7, 6, 0, 0},
  35                                 {100, 200, 300, 0, 0},
  36                                 {100, 300, 200, 0, 0},
  37                                 {200, 100, 300, 0, 0},
  38                                 {200, 300, 100, 0, 0},
  39                                 {300, 100, 200, 0, 0},
  40                                 {300, 200, 100, 0, 0},
  41                                 {100, 200, 300, 400, 500},
  42                                 {100, 300, 200, 400, 500},
  43                                 {200, 100, 300, 400, 500},
  44                                 {200, 300, 100, 400, 500},
  45                                 {300, 100, 200, 400, 500},
  46                                 {300, 200, 100, 400, 500},
  47                                 {100, 200, 300, 500, 400},
  48                                 {100, 300, 200, 500, 400},
  49                                 {200, 100, 300, 500, 400},
  50                                 {200, 300, 100, 500, 400},
  51                                 {300, 100, 200, 500, 400},
  52                                 {300, 200, 100, 500, 400},
  53                         } {
  54                                 var ma, na, mc, nc int
  55                                 if side == blas.Left {
  56                                         ma = test.common
  57                                         na = test.adim
  58                                         mc = test.common
  59                                         nc = test.cdim
  60                                 } else {
  61                                         ma = test.common
  62                                         na = test.adim
  63                                         mc = test.cdim
  64                                         nc = test.common
  65                                 }
  66                                 // Generate a random matrix
  67                                 lda := test.lda
  68                                 if lda == 0 {
  69                                         lda = na
  70                                 }
  71                                 a := make([]float64, ma*lda)
  72                                 for i := range a {
  73                                         a[i] = rnd.Float64()
  74                                 }
  75                                 // Compute random C matrix
  76                                 ldc := test.ldc
  77                                 if ldc == 0 {
  78                                         ldc = nc
  79                                 }
  80                                 c := make([]float64, mc*ldc)
  81                                 for i := range c {
  82                                         c[i] = rnd.Float64()
  83                                 }
  84
  85                                 // Compute QR
  86                                 k := min(ma, na)
  87                                 tau := make([]float64, k)
  88                                 work := make([]float64, 1)
  89                                 impl.Dgeqrf(ma, na, a, lda, tau, work, -1)
  90                                 work = make([]float64, int(work[0]))
  91                                 impl.Dgeqrf(ma, na, a, lda, tau, work, len(work))
  92
  93                                 cCopy := make([]float64, len(c))
  94                                 copy(cCopy, c)
  95                                 ans := make([]float64, len(c))
  96                                 copy(ans, cCopy)
  97
  98                                 if side == blas.Left {
  99                                         work = make([]float64, nc)
 100                                 } else {
 101                                         work = make([]float64, mc)
 102                                 }
 103                                 impl.Dorm2r(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
 104
 105                                 // Make sure Dorm2r and Dormqr match with small work
 106                                 for i := range work {
 107                                         work[i] = rnd.Float64()
 108                                 }
 109                                 copy(c, cCopy)
 110                                 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
 111                                 if !floats.EqualApprox(c, ans, 1e-12) {
 112                                         t.Errorf("Dormqr and Dorm2r mismatch for small work")
 113                                 }
 114
 115                                 // Try with the optimum amount of work
 116                                 copy(c, cCopy)
 117                                 impl.Dormqr(side, trans, mc, nc, k, nil, lda, nil, nil, ldc, work, -1)
 118                                 work = make([]float64, int(work[0]))
 119                                 for i := range work {
 120                                         work[i] = rnd.Float64()
 121                                 }
 122                                 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
 123                                 if !floats.EqualApprox(c, ans, 1e-12) {
 124                                         t.Errorf("Dormqr and Dorm2r mismatch for full work")
 125                                         fmt.Println("ccopy")
 126                                         for i := 0; i < mc; i++ {
 127                                                 fmt.Println(cCopy[i*ldc : (i+1)*ldc])
 128                                         }
 129                                         fmt.Println("ans =")
 130                                         for i := 0; i < mc; i++ {
 131                                                 fmt.Println(ans[i*ldc : (i+1)*ldc])
 132                                         }
 133                                         fmt.Println("c =")
 134                                         for i := 0; i < mc; i++ {
 135                                                 fmt.Println(c[i*ldc : (i+1)*ldc])
 136                                         }
 137                                 }
 138
 139                                 // Try with amount of work that is less than
 140                                 // optimal but still long enough to use the
 141                                 // blocked code.
 142                                 copy(c, cCopy)
 143                                 if side == blas.Left {
 144                                         work = make([]float64, 3*nc)
 145                                 } else {
 146                                         work = make([]float64, 3*mc)
 147                                 }
 148                                 impl.Dormqr(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, len(work))
 149                                 if !floats.EqualApprox(c, ans, 1e-12) {
 150                                         t.Errorf("Dormqr and Dorm2r mismatch for medium work")
 151                                 }
 152                         }
 153                 }
 154         }
 155 }