vendor/gonum.org/v1/gonum/lapack/testlapack/dgehrd.go

   1 // Copyright ©2016 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         "math"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/blas"
  15         "gonum.org/v1/gonum/blas/blas64"
  16 )
  17
  18 type Dgehrder interface {
  19         Dgehrd(n, ilo, ihi int, a []float64, lda int, tau, work []float64, lwork int)
  20
  21         Dorgqr(m, n, k int, a []float64, lda int, tau, work []float64, lwork int)
  22 }
  23
  24 func DgehrdTest(t *testing.T, impl Dgehrder) {
  25         rnd := rand.New(rand.NewSource(1))
  26
  27         // Randomized tests for small matrix sizes that will most likely
  28         // use the unblocked algorithm.
  29         for _, n := range []int{1, 2, 3, 4, 5, 10, 34} {
  30                 for _, extra := range []int{0, 13} {
  31                         for _, optwork := range []bool{true, false} {
  32                                 for cas := 0; cas < 10; cas++ {
  33                                         ilo := rnd.Intn(n)
  34                                         ihi := rnd.Intn(n)
  35                                         if ilo > ihi {
  36                                                 ilo, ihi = ihi, ilo
  37                                         }
  38                                         testDgehrd(t, impl, n, ilo, ihi, extra, optwork, rnd)
  39                                 }
  40                         }
  41                 }
  42         }
  43
  44         // These are selected tests for larger matrix sizes to test the blocked
  45         // algorithm. Use sizes around several powers of two because that is
  46         // where the blocked path will most likely start to be taken. For
  47         // example, at present the blocked algorithm is used for sizes larger
  48         // than 129.
  49         for _, test := range []struct {
  50                 n, ilo, ihi int
  51         }{
  52                 {0, 0, -1},
  53
  54                 {68, 0, 63},
  55                 {68, 0, 64},
  56                 {68, 0, 65},
  57                 {68, 0, 66},
  58                 {68, 0, 67},
  59
  60                 {132, 2, 129},
  61                 {132, 1, 129}, // Size = 129, unblocked.
  62                 {132, 0, 129}, // Size = 130, blocked.
  63                 {132, 1, 130},
  64                 {132, 0, 130},
  65                 {132, 1, 131},
  66                 {132, 0, 131},
  67
  68                 {260, 2, 257},
  69                 {260, 1, 257},
  70                 {260, 0, 257},
  71                 {260, 0, 258},
  72                 {260, 0, 259},
  73         } {
  74                 for _, extra := range []int{0, 13} {
  75                         for _, optwork := range []bool{true, false} {
  76                                 testDgehrd(t, impl, test.n, test.ilo, test.ihi, extra, optwork, rnd)
  77                         }
  78                 }
  79         }
  80 }
  81
  82 func testDgehrd(t *testing.T, impl Dgehrder, n, ilo, ihi, extra int, optwork bool, rnd *rand.Rand) {
  83         a := randomGeneral(n, n, n+extra, rnd)
  84         aCopy := a
  85         aCopy.Data = make([]float64, len(a.Data))
  86         copy(aCopy.Data, a.Data)
  87
  88         var tau []float64
  89         if n > 1 {
  90                 tau = nanSlice(n - 1)
  91         }
  92
  93         var work []float64
  94         if optwork {
  95                 work = nanSlice(1)
  96                 impl.Dgehrd(n, ilo, ihi, nil, a.Stride, nil, work, -1)
  97                 work = nanSlice(int(work[0]))
  98         } else {
  99                 work = nanSlice(max(1, n))
 100         }
 101
 102         impl.Dgehrd(n, ilo, ihi, a.Data, a.Stride, tau, work, len(work))
 103
 104         if n == 0 {
 105                 // Just make sure there is no panic.
 106                 return
 107         }
 108
 109         prefix := fmt.Sprintf("Case n=%v, ilo=%v, ihi=%v, extra=%v", n, ilo, ihi, extra)
 110
 111         // Check any invalid modifications of a.
 112         if !generalOutsideAllNaN(a) {
 113                 t.Errorf("%v: out-of-range write to A\n%v", prefix, a.Data)
 114         }
 115         for i := ilo; i <= ihi; i++ {
 116                 for j := 0; j < min(ilo, i); j++ {
 117                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
 118                                 t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j)
 119                         }
 120                 }
 121         }
 122         for i := ihi + 1; i < n; i++ {
 123                 for j := 0; j < i; j++ {
 124                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
 125                                 t.Errorf("%v: unexpected modification of A[%v,%v]", prefix, i, j)
 126                         }
 127                 }
 128         }
 129         for i := 0; i <= ilo; i++ {
 130                 for j := i; j < ilo+1; j++ {
 131                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
 132                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
 133                         }
 134                 }
 135                 for j := ihi + 1; j < n; j++ {
 136                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
 137                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
 138                         }
 139                 }
 140         }
 141         for i := ihi + 1; i < n; i++ {
 142                 for j := i; j < n; j++ {
 143                         if a.Data[i*a.Stride+j] != aCopy.Data[i*aCopy.Stride+j] {
 144                                 t.Errorf("%v: unexpected modification at A[%v,%v]", prefix, i, j)
 145                         }
 146                 }
 147         }
 148
 149         // Check that tau has been assigned properly.
 150         for i, v := range tau {
 151                 if math.IsNaN(v) {
 152                         t.Errorf("%v: unexpected NaN at tau[%v]", prefix, i)
 153                 }
 154         }
 155
 156         // Extract Q and check that it is orthogonal.
 157         q := eye(n, n)
 158         if ilo != ihi {
 159                 for i := ilo + 2; i <= ihi; i++ {
 160                         for j := ilo + 1; j < ihi; j++ {
 161                                 q.Data[i*q.Stride+j] = a.Data[i*a.Stride+j-1]
 162                         }
 163                 }
 164                 nh := ihi - ilo
 165                 impl.Dorgqr(nh, nh, nh, q.Data[(ilo+1)*q.Stride+ilo+1:], q.Stride, tau[ilo:ihi], work, len(work))
 166         }
 167         if !isOrthonormal(q) {
 168                 t.Errorf("%v: Q is not orthogonal\nQ=%v", prefix, q)
 169         }
 170
 171         // Construct Q^T * AOrig * Q and check that it is upper Hessenberg.
 172         aq := blas64.General{
 173                 Rows:   n,
 174                 Cols:   n,
 175                 Stride: n,
 176                 Data:   make([]float64, n*n),
 177         }
 178         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, aCopy, q, 0, aq)
 179         qaq := blas64.General{
 180                 Rows:   n,
 181                 Cols:   n,
 182                 Stride: n,
 183                 Data:   make([]float64, n*n),
 184         }
 185         blas64.Gemm(blas.Trans, blas.NoTrans, 1, q, aq, 0, qaq)
 186         for i := 0; i <= ilo; i++ {
 187                 for j := ilo + 1; j <= ihi; j++ {
 188                         qaqij := qaq.Data[i*qaq.Stride+j]
 189                         diff := qaqij - a.Data[i*a.Stride+j]
 190                         if math.Abs(diff) > 1e-13 {
 191                                 t.Errorf("%v: Q^T*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
 192                         }
 193                 }
 194         }
 195         for i := ilo + 1; i <= ihi; i++ {
 196                 for j := ilo; j < n; j++ {
 197                         qaqij := qaq.Data[i*qaq.Stride+j]
 198                         if j < i-1 {
 199                                 if math.Abs(qaqij) > 1e-13 {
 200                                         t.Errorf("%v: Q^T*AOrig*Q is not upper Hessenberg, [%v,%v]=%v", prefix, i, j, qaqij)
 201                                 }
 202                                 continue
 203                         }
 204                         diff := qaqij - a.Data[i*a.Stride+j]
 205                         if math.Abs(diff) > 1e-13 {
 206                                 t.Errorf("%v: Q^T*AOrig*Q and A are not equal, diff at [%v,%v]=%v", prefix, i, j, diff)
 207                         }
 208                 }
 209         }
 210 }