vendor/gonum.org/v1/gonum/lapack/testlapack/dormlq.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/floats"
  14 )
  15
  16 type Dormlqer interface {
  17         Dorml2er
  18         Dormlq(side blas.Side, trans blas.Transpose, m, n, k int, a []float64, lda int, tau, c []float64, ldc int, work []float64, lwork int)
  19 }
  20
  21 func DormlqTest(t *testing.T, impl Dormlqer) {
  22         rnd := rand.New(rand.NewSource(1))
  23         for _, side := range []blas.Side{blas.Left, blas.Right} {
  24                 for _, trans := range []blas.Transpose{blas.NoTrans, blas.Trans} {
  25                         for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
  26                                 for _, test := range []struct {
  27                                         common, adim, cdim, lda, ldc int
  28                                 }{
  29                                         {0, 0, 0, 0, 0},
  30                                         {6, 7, 8, 0, 0},
  31                                         {6, 8, 7, 0, 0},
  32                                         {7, 6, 8, 0, 0},
  33                                         {7, 8, 6, 0, 0},
  34                                         {8, 6, 7, 0, 0},
  35                                         {8, 7, 6, 0, 0},
  36                                         {100, 200, 300, 0, 0},
  37                                         {100, 300, 200, 0, 0},
  38                                         {200, 100, 300, 0, 0},
  39                                         {200, 300, 100, 0, 0},
  40                                         {300, 100, 200, 0, 0},
  41                                         {300, 200, 100, 0, 0},
  42                                         {100, 200, 300, 400, 500},
  43                                         {100, 300, 200, 400, 500},
  44                                         {200, 100, 300, 400, 500},
  45                                         {200, 300, 100, 400, 500},
  46                                         {300, 100, 200, 400, 500},
  47                                         {300, 200, 100, 400, 500},
  48                                         {100, 200, 300, 500, 400},
  49                                         {100, 300, 200, 500, 400},
  50                                         {200, 100, 300, 500, 400},
  51                                         {200, 300, 100, 500, 400},
  52                                         {300, 100, 200, 500, 400},
  53                                         {300, 200, 100, 500, 400},
  54                                 } {
  55                                         var ma, na, mc, nc int
  56                                         if side == blas.Left {
  57                                                 ma = test.adim
  58                                                 na = test.common
  59                                                 mc = test.common
  60                                                 nc = test.cdim
  61                                         } else {
  62                                                 ma = test.adim
  63                                                 na = test.common
  64                                                 mc = test.cdim
  65                                                 nc = test.common
  66                                         }
  67                                         // Generate a random matrix
  68                                         lda := test.lda
  69                                         if lda == 0 {
  70                                                 lda = na
  71                                         }
  72                                         a := make([]float64, ma*lda)
  73                                         for i := range a {
  74                                                 a[i] = rnd.Float64()
  75                                         }
  76                                         // Compute random C matrix
  77                                         ldc := test.ldc
  78                                         if ldc == 0 {
  79                                                 ldc = nc
  80                                         }
  81                                         c := make([]float64, mc*ldc)
  82                                         for i := range c {
  83                                                 c[i] = rnd.Float64()
  84                                         }
  85
  86                                         // Compute LQ
  87                                         k := min(ma, na)
  88                                         tau := make([]float64, k)
  89                                         work := make([]float64, 1)
  90                                         impl.Dgelqf(ma, na, a, lda, tau, work, -1)
  91                                         work = make([]float64, int(work[0]))
  92                                         impl.Dgelqf(ma, na, a, lda, tau, work, len(work))
  93
  94                                         cCopy := make([]float64, len(c))
  95                                         copy(cCopy, c)
  96                                         ans := make([]float64, len(c))
  97                                         copy(ans, cCopy)
  98
  99                                         var nw int
 100                                         if side == blas.Left {
 101                                                 nw = nc
 102                                         } else {
 103                                                 nw = mc
 104                                         }
 105                                         work = make([]float64, max(1, nw))
 106                                         impl.Dorml2(side, trans, mc, nc, k, a, lda, tau, ans, ldc, work)
 107
 108                                         var lwork int
 109                                         switch wl {
 110                                         case minimumWork:
 111                                                 lwork = nw
 112                                         case optimumWork:
 113                                                 impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, -1)
 114                                                 lwork = int(work[0])
 115                                         case mediumWork:
 116                                                 work := make([]float64, 1)
 117                                                 impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, -1)
 118                                                 lwork = (int(work[0]) + nw) / 2
 119                                         }
 120                                         lwork = max(1, lwork)
 121                                         work = make([]float64, lwork)
 122
 123                                         impl.Dormlq(side, trans, mc, nc, k, a, lda, tau, c, ldc, work, lwork)
 124                                         if !floats.EqualApprox(c, ans, 1e-13) {
 125                                                 t.Errorf("Dormqr and Dorm2r results mismatch")
 126                                         }
 127                                 }
 128                         }
 129                 }
 130         }
 131 }