vendor/gonum.org/v1/gonum/lapack/testlapack/dlaqr23.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         "testing"
  10
  11         "golang.org/x/exp/rand"
  12
  13         "gonum.org/v1/gonum/blas"
  14         "gonum.org/v1/gonum/blas/blas64"
  15 )
  16
  17 type Dlaqr23er interface {
  18         Dlaqr23(wantt, wantz bool, n, ktop, kbot, nw int, h []float64, ldh int, iloz, ihiz int, z []float64, ldz int, sr, si []float64, v []float64, ldv int, nh int, t []float64, ldt int, nv int, wv []float64, ldwv int, work []float64, lwork int, recur int) (ns, nd int)
  19 }
  20
  21 type dlaqr23Test struct {
  22         wantt, wantz bool
  23         ktop, kbot   int
  24         nw           int
  25         h            blas64.General
  26         iloz, ihiz   int
  27
  28         evWant []complex128 // Optional slice with known eigenvalues.
  29 }
  30
  31 func Dlaqr23Test(t *testing.T, impl Dlaqr23er) {
  32         rnd := rand.New(rand.NewSource(1))
  33
  34         for _, wantt := range []bool{true, false} {
  35                 for _, wantz := range []bool{true, false} {
  36                         for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 18, 31, 100} {
  37                                 for _, extra := range []int{0, 11} {
  38                                         for cas := 0; cas < 30; cas++ {
  39                                                 var nw int
  40                                                 if nw <= 75 {
  41                                                         nw = rnd.Intn(n) + 1
  42                                                 } else {
  43                                                         nw = 76 + rnd.Intn(n-75)
  44                                                 }
  45                                                 ktop := rnd.Intn(n - nw + 1)
  46                                                 kbot := ktop + nw - 1
  47                                                 kbot += rnd.Intn(n - kbot)
  48                                                 h := randomHessenberg(n, n+extra, rnd)
  49                                                 if ktop-1 >= 0 {
  50                                                         h.Data[ktop*h.Stride+ktop-1] = 0
  51                                                 }
  52                                                 if kbot+1 < n {
  53                                                         h.Data[(kbot+1)*h.Stride+kbot] = 0
  54                                                 }
  55                                                 iloz := rnd.Intn(ktop + 1)
  56                                                 ihiz := kbot + rnd.Intn(n-kbot)
  57                                                 test := dlaqr23Test{
  58                                                         wantt: wantt,
  59                                                         wantz: wantz,
  60                                                         ktop:  ktop,
  61                                                         kbot:  kbot,
  62                                                         nw:    nw,
  63                                                         h:     h,
  64                                                         iloz:  iloz,
  65                                                         ihiz:  ihiz,
  66                                                 }
  67                                                 testDlaqr23(t, impl, test, false, 1, rnd)
  68                                                 testDlaqr23(t, impl, test, true, 1, rnd)
  69                                                 testDlaqr23(t, impl, test, false, 0, rnd)
  70                                                 testDlaqr23(t, impl, test, true, 0, rnd)
  71                                         }
  72                                 }
  73                         }
  74                 }
  75         }
  76
  77         // Tests with n=0.
  78         for _, wantt := range []bool{true, false} {
  79                 for _, wantz := range []bool{true, false} {
  80                         for _, extra := range []int{0, 1, 11} {
  81                                 test := dlaqr23Test{
  82                                         wantt: wantt,
  83                                         wantz: wantz,
  84                                         h:     randomHessenberg(0, extra, rnd),
  85                                         ktop:  0,
  86                                         kbot:  -1,
  87                                         iloz:  0,
  88                                         ihiz:  -1,
  89                                         nw:    0,
  90                                 }
  91                                 testDlaqr23(t, impl, test, true, 1, rnd)
  92                                 testDlaqr23(t, impl, test, false, 1, rnd)
  93                                 testDlaqr23(t, impl, test, true, 0, rnd)
  94                                 testDlaqr23(t, impl, test, false, 0, rnd)
  95                         }
  96                 }
  97         }
  98
  99         // Tests with explicit eigenvalues computed by Octave.
 100         for _, test := range []dlaqr23Test{
 101                 {
 102                         h: blas64.General{
 103                                 Rows:   1,
 104                                 Cols:   1,
 105                                 Stride: 1,
 106                                 Data:   []float64{7.09965484086874e-1},
 107                         },
 108                         ktop:   0,
 109                         kbot:   0,
 110                         iloz:   0,
 111                         ihiz:   0,
 112                         evWant: []complex128{7.09965484086874e-1},
 113                 },
 114                 {
 115                         h: blas64.General{
 116                                 Rows:   2,
 117                                 Cols:   2,
 118                                 Stride: 2,
 119                                 Data: []float64{
 120                                         0, -1,
 121                                         1, 0,
 122                                 },
 123                         },
 124                         ktop:   0,
 125                         kbot:   1,
 126                         iloz:   0,
 127                         ihiz:   1,
 128                         evWant: []complex128{1i, -1i},
 129                 },
 130                 {
 131                         h: blas64.General{
 132                                 Rows:   2,
 133                                 Cols:   2,
 134                                 Stride: 2,
 135                                 Data: []float64{
 136                                         6.25219991450918e-1, 8.17510791994361e-1,
 137                                         3.31218891622294e-1, 1.24103744878131e-1,
 138                                 },
 139                         },
 140                         ktop:   0,
 141                         kbot:   1,
 142                         iloz:   0,
 143                         ihiz:   1,
 144                         evWant: []complex128{9.52203547663447e-1, -2.02879811334398e-1},
 145                 },
 146                 {
 147                         h: blas64.General{
 148                                 Rows:   4,
 149                                 Cols:   4,
 150                                 Stride: 4,
 151                                 Data: []float64{
 152                                         1, 0, 0, 0,
 153                                         0, 6.25219991450918e-1, 8.17510791994361e-1, 0,
 154                                         0, 3.31218891622294e-1, 1.24103744878131e-1, 0,
 155                                         0, 0, 0, 1,
 156                                 },
 157                         },
 158                         ktop:   1,
 159                         kbot:   2,
 160                         iloz:   0,
 161                         ihiz:   3,
 162                         evWant: []complex128{9.52203547663447e-1, -2.02879811334398e-1},
 163                 },
 164                 {
 165                         h: blas64.General{
 166                                 Rows:   2,
 167                                 Cols:   2,
 168                                 Stride: 2,
 169                                 Data: []float64{
 170                                         -1.1219562276608, 6.85473513349362e-1,
 171                                         -8.19951061145131e-1, 1.93728523178888e-1,
 172                                 },
 173                         },
 174                         ktop: 0,
 175                         kbot: 1,
 176                         iloz: 0,
 177                         ihiz: 1,
 178                         evWant: []complex128{
 179                                 -4.64113852240958e-1 + 3.59580510817350e-1i,
 180                                 -4.64113852240958e-1 - 3.59580510817350e-1i,
 181                         },
 182                 },
 183                 {
 184                         h: blas64.General{
 185                                 Rows:   5,
 186                                 Cols:   5,
 187                                 Stride: 5,
 188                                 Data: []float64{
 189                                         9.57590178533658e-1, -5.10651295522708e-1, 9.24974510015869e-1, -1.30016306879522e-1, 2.92601986926954e-2,
 190                                         -1.08084756637964, 1.77529701001213, -1.36480197632509, 2.23196371219601e-1, 1.12912853063308e-1,
 191                                         0, -8.44075612174676e-1, 1.067867614486, -2.55782915176399e-1, -2.00598563137468e-1,
 192                                         0, 0, -5.67097237165410e-1, 2.07205057427341e-1, 6.54998340743380e-1,
 193                                         0, 0, 0, -1.89441413886041e-1, -4.18125416021786e-1,
 194                                 },
 195                         },
 196                         ktop: 0,
 197                         kbot: 4,
 198                         iloz: 0,
 199                         ihiz: 4,
 200                         evWant: []complex128{
 201                                 2.94393309555622,
 202                                 4.97029793606701e-1 + 3.63041654992384e-1i,
 203                                 4.97029793606701e-1 - 3.63041654992384e-1i,
 204                                 -1.74079119166145e-1 + 2.01570009462092e-1i,
 205                                 -1.74079119166145e-1 - 2.01570009462092e-1i,
 206                         },
 207                 },
 208         } {
 209                 test.wantt = true
 210                 test.wantz = true
 211                 test.nw = test.kbot - test.ktop + 1
 212                 testDlaqr23(t, impl, test, true, 1, rnd)
 213                 testDlaqr23(t, impl, test, false, 1, rnd)
 214                 testDlaqr23(t, impl, test, true, 0, rnd)
 215                 testDlaqr23(t, impl, test, false, 0, rnd)
 216         }
 217 }
 218
 219 func testDlaqr23(t *testing.T, impl Dlaqr23er, test dlaqr23Test, opt bool, recur int, rnd *rand.Rand) {
 220         const tol = 1e-14
 221
 222         h := cloneGeneral(test.h)
 223         n := h.Cols
 224         extra := h.Stride - h.Cols
 225         wantt := test.wantt
 226         wantz := test.wantz
 227         ktop := test.ktop
 228         kbot := test.kbot
 229         nw := test.nw
 230         iloz := test.iloz
 231         ihiz := test.ihiz
 232
 233         var z, zCopy blas64.General
 234         if wantz {
 235                 z = eye(n, n+extra)
 236                 zCopy = cloneGeneral(z)
 237         }
 238
 239         sr := nanSlice(kbot + 1)
 240         si := nanSlice(kbot + 1)
 241
 242         v := randomGeneral(nw, nw, nw+extra, rnd)
 243         var nh int
 244         if nw > 0 {
 245                 nh = nw + rnd.Intn(nw) // nh must be at least nw.
 246         }
 247         tmat := randomGeneral(nw, nh, nh+extra, rnd)
 248         var nv int
 249         if nw > 0 {
 250                 nv = rnd.Intn(nw) + 1
 251         }
 252         wv := randomGeneral(nv, nw, nw+extra, rnd)
 253
 254         var work []float64
 255         if opt {
 256                 work = nanSlice(1)
 257                 impl.Dlaqr23(wantt, wantz, n, ktop, kbot, nw, nil, h.Stride, iloz, ihiz, nil, z.Stride,
 258                         nil, nil, nil, v.Stride, tmat.Cols, nil, tmat.Stride, wv.Rows, nil, wv.Stride, work, -1, recur)
 259                 work = nanSlice(int(work[0]))
 260         } else {
 261                 work = nanSlice(max(1, 2*nw))
 262         }
 263
 264         ns, nd := impl.Dlaqr23(wantt, wantz, n, ktop, kbot, nw, h.Data, h.Stride, iloz, ihiz, z.Data, z.Stride,
 265                 sr, si, v.Data, v.Stride, tmat.Cols, tmat.Data, tmat.Stride, wv.Rows, wv.Data, wv.Stride, work, len(work), recur)
 266
 267         prefix := fmt.Sprintf("Case wantt=%v, wantz=%v, n=%v, ktop=%v, kbot=%v, nw=%v, iloz=%v, ihiz=%v, extra=%v",
 268                 wantt, wantz, n, ktop, kbot, nw, iloz, ihiz, extra)
 269
 270         if !generalOutsideAllNaN(h) {
 271                 t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
 272         }
 273         if !generalOutsideAllNaN(z) {
 274                 t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
 275         }
 276         if !generalOutsideAllNaN(v) {
 277                 t.Errorf("%v: out-of-range write to V\n%v", prefix, v.Data)
 278         }
 279         if !generalOutsideAllNaN(tmat) {
 280                 t.Errorf("%v: out-of-range write to T\n%v", prefix, tmat.Data)
 281         }
 282         if !generalOutsideAllNaN(wv) {
 283                 t.Errorf("%v: out-of-range write to WV\n%v", prefix, wv.Data)
 284         }
 285         if !isAllNaN(sr[:kbot-nd-ns+1]) || !isAllNaN(sr[kbot+1:]) {
 286                 t.Errorf("%v: out-of-range write to sr", prefix)
 287         }
 288         if !isAllNaN(si[:kbot-nd-ns+1]) || !isAllNaN(si[kbot+1:]) {
 289                 t.Errorf("%v: out-of-range write to si", prefix)
 290         }
 291
 292         if !isUpperHessenberg(h) {
 293                 t.Errorf("%v: H is not upper Hessenberg", prefix)
 294         }
 295
 296         if test.evWant != nil {
 297                 for i := kbot - nd + 1; i <= kbot; i++ {
 298                         ev := complex(sr[i], si[i])
 299                         found, _ := containsComplex(test.evWant, ev, tol)
 300                         if !found {
 301                                 t.Errorf("%v: unexpected eigenvalue %v", prefix, ev)
 302                         }
 303                 }
 304         }
 305
 306         if !wantz {
 307                 return
 308         }
 309
 310         var zmod bool
 311         for i := 0; i < n; i++ {
 312                 for j := 0; j < n; j++ {
 313                         if z.Data[i*z.Stride+j] == zCopy.Data[i*zCopy.Stride+j] {
 314                                 continue
 315                         }
 316                         if i < iloz || i > ihiz || j < kbot-nw+1 || j > kbot {
 317                                 zmod = true
 318                         }
 319                 }
 320         }
 321         if zmod {
 322                 t.Errorf("%v: unexpected modification of Z", prefix)
 323         }
 324         if !isOrthonormal(z) {
 325                 t.Errorf("%v: Z is not orthogonal", prefix)
 326         }
 327         if wantt {
 328                 hu := eye(n, n)
 329                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, test.h, z, 0, hu)
 330                 uhu := eye(n, n)
 331                 blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hu, 0, uhu)
 332                 if !equalApproxGeneral(uhu, h, 10*tol) {
 333                         t.Errorf("%v: Z^T*(initial H)*Z and (final H) are not equal", prefix)
 334                 }
 335         }
 336 }