vendor/gonum.org/v1/gonum/lapack/testlapack/dlaqr04.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 Dlaqr04er interface {
  19         Dlaqr04(wantt, wantz bool, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, iloz, ihiz int, z []float64, ldz int, work []float64, lwork int, recur int) int
  20
  21         Dlahqrer
  22 }
  23
  24 type dlaqr04Test struct {
  25         h            blas64.General
  26         ilo, ihi     int
  27         iloz, ihiz   int
  28         wantt, wantz bool
  29
  30         evWant []complex128 // Optional slice holding known eigenvalues.
  31 }
  32
  33 func Dlaqr04Test(t *testing.T, impl Dlaqr04er) {
  34         rnd := rand.New(rand.NewSource(1))
  35
  36         // Tests for small matrices that choose the ilo,ihi and iloz,ihiz pairs
  37         // randomly.
  38         for _, wantt := range []bool{true, false} {
  39                 for _, wantz := range []bool{true, false} {
  40                         for _, n := range []int{1, 2, 3, 4, 5, 6, 10, 11, 12, 18, 29} {
  41                                 for _, extra := range []int{0, 11} {
  42                                         for recur := 0; recur <= 2; recur++ {
  43                                                 for cas := 0; cas < n; cas++ {
  44                                                         ilo := rnd.Intn(n)
  45                                                         ihi := rnd.Intn(n)
  46                                                         if ilo > ihi {
  47                                                                 ilo, ihi = ihi, ilo
  48                                                         }
  49                                                         iloz := rnd.Intn(ilo + 1)
  50                                                         ihiz := ihi + rnd.Intn(n-ihi)
  51                                                         h := randomHessenberg(n, n+extra, rnd)
  52                                                         if ilo-1 >= 0 {
  53                                                                 h.Data[ilo*h.Stride+ilo-1] = 0
  54                                                         }
  55                                                         if ihi+1 < n {
  56                                                                 h.Data[(ihi+1)*h.Stride+ihi] = 0
  57                                                         }
  58                                                         test := dlaqr04Test{
  59                                                                 h:     h,
  60                                                                 ilo:   ilo,
  61                                                                 ihi:   ihi,
  62                                                                 iloz:  iloz,
  63                                                                 ihiz:  ihiz,
  64                                                                 wantt: wantt,
  65                                                                 wantz: wantz,
  66                                                         }
  67                                                         testDlaqr04(t, impl, test, false, recur)
  68                                                         testDlaqr04(t, impl, test, true, recur)
  69                                                 }
  70                                         }
  71                                 }
  72                         }
  73                 }
  74         }
  75
  76         // Tests for matrices large enough to possibly use the recursion (but it
  77         // doesn't seem to be the case).
  78         for _, n := range []int{100, 500} {
  79                 for cas := 0; cas < 5; cas++ {
  80                         h := randomHessenberg(n, n, rnd)
  81                         test := dlaqr04Test{
  82                                 h:     h,
  83                                 ilo:   0,
  84                                 ihi:   n - 1,
  85                                 iloz:  0,
  86                                 ihiz:  n - 1,
  87                                 wantt: true,
  88                                 wantz: true,
  89                         }
  90                         testDlaqr04(t, impl, test, true, 1)
  91                 }
  92         }
  93
  94         // Tests that make sure that some potentially problematic corner cases,
  95         // like zero-sized matrix, are covered.
  96         for _, wantt := range []bool{true, false} {
  97                 for _, wantz := range []bool{true, false} {
  98                         for _, extra := range []int{0, 1, 11} {
  99                                 for _, test := range []dlaqr04Test{
 100                                         {
 101                                                 h:    randomHessenberg(0, extra, rnd),
 102                                                 ilo:  0,
 103                                                 ihi:  -1,
 104                                                 iloz: 0,
 105                                                 ihiz: -1,
 106                                         },
 107                                         {
 108                                                 h:    randomHessenberg(1, 1+extra, rnd),
 109                                                 ilo:  0,
 110                                                 ihi:  0,
 111                                                 iloz: 0,
 112                                                 ihiz: 0,
 113                                         },
 114                                         {
 115                                                 h:    randomHessenberg(2, 2+extra, rnd),
 116                                                 ilo:  1,
 117                                                 ihi:  1,
 118                                                 iloz: 1,
 119                                                 ihiz: 1,
 120                                         },
 121                                         {
 122                                                 h:    randomHessenberg(2, 2+extra, rnd),
 123                                                 ilo:  0,
 124                                                 ihi:  1,
 125                                                 iloz: 0,
 126                                                 ihiz: 1,
 127                                         },
 128                                         {
 129                                                 h:    randomHessenberg(10, 10+extra, rnd),
 130                                                 ilo:  0,
 131                                                 ihi:  0,
 132                                                 iloz: 0,
 133                                                 ihiz: 0,
 134                                         },
 135                                         {
 136                                                 h:    randomHessenberg(10, 10+extra, rnd),
 137                                                 ilo:  0,
 138                                                 ihi:  9,
 139                                                 iloz: 0,
 140                                                 ihiz: 9,
 141                                         },
 142                                         {
 143                                                 h:    randomHessenberg(10, 10+extra, rnd),
 144                                                 ilo:  0,
 145                                                 ihi:  1,
 146                                                 iloz: 0,
 147                                                 ihiz: 1,
 148                                         },
 149                                         {
 150                                                 h:    randomHessenberg(10, 10+extra, rnd),
 151                                                 ilo:  0,
 152                                                 ihi:  1,
 153                                                 iloz: 0,
 154                                                 ihiz: 9,
 155                                         },
 156                                         {
 157                                                 h:    randomHessenberg(10, 10+extra, rnd),
 158                                                 ilo:  9,
 159                                                 ihi:  9,
 160                                                 iloz: 0,
 161                                                 ihiz: 9,
 162                                         },
 163                                 } {
 164                                         if test.ilo-1 >= 0 {
 165                                                 test.h.Data[test.ilo*test.h.Stride+test.ilo-1] = 0
 166                                         }
 167                                         if test.ihi+1 < test.h.Rows {
 168                                                 test.h.Data[(test.ihi+1)*test.h.Stride+test.ihi] = 0
 169                                         }
 170                                         test.wantt = wantt
 171                                         test.wantz = wantz
 172                                         testDlaqr04(t, impl, test, false, 1)
 173                                         testDlaqr04(t, impl, test, true, 1)
 174                                 }
 175                         }
 176                 }
 177         }
 178
 179         // Tests with known eigenvalues computed by Octave.
 180         for _, test := range []dlaqr04Test{
 181                 {
 182                         h: blas64.General{
 183                                 Rows:   1,
 184                                 Cols:   1,
 185                                 Stride: 1,
 186                                 Data:   []float64{7.09965484086874e-1},
 187                         },
 188                         ilo:    0,
 189                         ihi:    0,
 190                         iloz:   0,
 191                         ihiz:   0,
 192                         evWant: []complex128{7.09965484086874e-1},
 193                 },
 194                 {
 195                         h: blas64.General{
 196                                 Rows:   2,
 197                                 Cols:   2,
 198                                 Stride: 2,
 199                                 Data: []float64{
 200                                         0, -1,
 201                                         1, 0,
 202                                 },
 203                         },
 204                         ilo:    0,
 205                         ihi:    1,
 206                         iloz:   0,
 207                         ihiz:   1,
 208                         evWant: []complex128{1i, -1i},
 209                 },
 210                 {
 211                         h: blas64.General{
 212                                 Rows:   2,
 213                                 Cols:   2,
 214                                 Stride: 2,
 215                                 Data: []float64{
 216                                         6.25219991450918e-1, 8.17510791994361e-1,
 217                                         3.31218891622294e-1, 1.24103744878131e-1,
 218                                 },
 219                         },
 220                         ilo:    0,
 221                         ihi:    1,
 222                         iloz:   0,
 223                         ihiz:   1,
 224                         evWant: []complex128{9.52203547663447e-1, -2.02879811334398e-1},
 225                 },
 226                 {
 227                         h: blas64.General{
 228                                 Rows:   4,
 229                                 Cols:   4,
 230                                 Stride: 4,
 231                                 Data: []float64{
 232                                         1, 0, 0, 0,
 233                                         0, 6.25219991450918e-1, 8.17510791994361e-1, 0,
 234                                         0, 3.31218891622294e-1, 1.24103744878131e-1, 0,
 235                                         0, 0, 0, 1,
 236                                 },
 237                         },
 238                         ilo:    1,
 239                         ihi:    2,
 240                         iloz:   0,
 241                         ihiz:   3,
 242                         evWant: []complex128{9.52203547663447e-1, -2.02879811334398e-1},
 243                 },
 244                 {
 245                         h: blas64.General{
 246                                 Rows:   2,
 247                                 Cols:   2,
 248                                 Stride: 2,
 249                                 Data: []float64{
 250                                         -1.1219562276608, 6.85473513349362e-1,
 251                                         -8.19951061145131e-1, 1.93728523178888e-1,
 252                                 },
 253                         },
 254                         ilo:  0,
 255                         ihi:  1,
 256                         iloz: 0,
 257                         ihiz: 1,
 258                         evWant: []complex128{
 259                                 -4.64113852240958e-1 + 3.59580510817350e-1i,
 260                                 -4.64113852240958e-1 - 3.59580510817350e-1i,
 261                         },
 262                 },
 263                 {
 264                         h: blas64.General{
 265                                 Rows:   5,
 266                                 Cols:   5,
 267                                 Stride: 5,
 268                                 Data: []float64{
 269                                         9.57590178533658e-1, -5.10651295522708e-1, 9.24974510015869e-1, -1.30016306879522e-1, 2.92601986926954e-2,
 270                                         -1.08084756637964, 1.77529701001213, -1.36480197632509, 2.23196371219601e-1, 1.12912853063308e-1,
 271                                         0, -8.44075612174676e-1, 1.067867614486, -2.55782915176399e-1, -2.00598563137468e-1,
 272                                         0, 0, -5.67097237165410e-1, 2.07205057427341e-1, 6.54998340743380e-1,
 273                                         0, 0, 0, -1.89441413886041e-1, -4.18125416021786e-1,
 274                                 },
 275                         },
 276                         ilo:  0,
 277                         ihi:  4,
 278                         iloz: 0,
 279                         ihiz: 4,
 280                         evWant: []complex128{
 281                                 2.94393309555622,
 282                                 4.97029793606701e-1 + 3.63041654992384e-1i,
 283                                 4.97029793606701e-1 - 3.63041654992384e-1i,
 284                                 -1.74079119166145e-1 + 2.01570009462092e-1i,
 285                                 -1.74079119166145e-1 - 2.01570009462092e-1i,
 286                         },
 287                 },
 288         } {
 289                 test.wantt = true
 290                 test.wantz = true
 291                 testDlaqr04(t, impl, test, false, 1)
 292                 testDlaqr04(t, impl, test, true, 1)
 293         }
 294 }
 295
 296 func testDlaqr04(t *testing.T, impl Dlaqr04er, test dlaqr04Test, optwork bool, recur int) {
 297         const tol = 1e-14
 298
 299         h := cloneGeneral(test.h)
 300         n := h.Cols
 301         extra := h.Stride - h.Cols
 302         wantt := test.wantt
 303         wantz := test.wantz
 304         ilo := test.ilo
 305         ihi := test.ihi
 306         iloz := test.iloz
 307         ihiz := test.ihiz
 308
 309         var z, zCopy blas64.General
 310         if wantz {
 311                 z = eye(n, n+extra)
 312                 zCopy = cloneGeneral(z)
 313         }
 314
 315         wr := nanSlice(ihi + 1)
 316         wi := nanSlice(ihi + 1)
 317
 318         var work []float64
 319         if optwork {
 320                 work = nanSlice(1)
 321                 impl.Dlaqr04(wantt, wantz, n, ilo, ihi, nil, 0, nil, nil, iloz, ihiz, nil, 0, work, -1, recur)
 322                 work = nanSlice(int(work[0]))
 323         } else {
 324                 work = nanSlice(max(1, n))
 325         }
 326
 327         unconverged := impl.Dlaqr04(wantt, wantz, n, ilo, ihi, h.Data, h.Stride, wr, wi, iloz, ihiz, z.Data, z.Stride, work, len(work), recur)
 328
 329         prefix := fmt.Sprintf("Case wantt=%v, wantz=%v, n=%v, ilo=%v, ihi=%v, iloz=%v, ihiz=%v, extra=%v, opt=%v",
 330                 wantt, wantz, n, ilo, ihi, iloz, ihiz, extra, optwork)
 331
 332         if !generalOutsideAllNaN(h) {
 333                 t.Errorf("%v: out-of-range write to H\n%v", prefix, h.Data)
 334         }
 335         if !generalOutsideAllNaN(z) {
 336                 t.Errorf("%v: out-of-range write to Z\n%v", prefix, z.Data)
 337         }
 338
 339         start := ilo // Index of the first computed eigenvalue.
 340         if unconverged != 0 {
 341                 start = unconverged
 342                 if start == ihi+1 {
 343                         t.Logf("%v: no eigenvalue has converged", prefix)
 344                 }
 345         }
 346
 347         // Check that wr and wi have not been modified within [:start].
 348         if !isAllNaN(wr[:start]) {
 349                 t.Errorf("%v: unexpected modification of wr", prefix)
 350         }
 351         if !isAllNaN(wi[:start]) {
 352                 t.Errorf("%v: unexpected modification of wi", prefix)
 353         }
 354
 355         var hasReal bool
 356         for i := start; i <= ihi; {
 357                 if wi[i] == 0 { // Real eigenvalue.
 358                         hasReal = true
 359                         // Check that the eigenvalue corresponds to a 1×1 block
 360                         // on the diagonal of H.
 361                         if wantt && wr[i] != h.Data[i*h.Stride+i] {
 362                                 t.Errorf("%v: wr[%v] != H[%v,%v]", prefix, i, i, i)
 363                         }
 364                         i++
 365                         continue
 366                 }
 367
 368                 // Complex eigenvalue.
 369
 370                 // In the conjugate pair the real parts must be equal.
 371                 if wr[i] != wr[i+1] {
 372                         t.Errorf("%v: real part of conjugate pair not equal, i=%v", prefix, i)
 373                 }
 374                 // The first imaginary part must be positive.
 375                 if wi[i] < 0 {
 376                         t.Errorf("%v: wi[%v] not positive", prefix, i)
 377                 }
 378                 // The second imaginary part must be negative with the same
 379                 // magnitude.
 380                 if wi[i] != -wi[i+1] {
 381                         t.Errorf("%v: wi[%v] != wi[%v]", prefix, i, i+1)
 382                 }
 383                 if wantt {
 384                         // Check that wi[i] has the correct value.
 385                         if wr[i] != h.Data[i*h.Stride+i] {
 386                                 t.Errorf("%v: wr[%v] != H[%v,%v]", prefix, i, i, i)
 387                         }
 388                         if wr[i] != h.Data[(i+1)*h.Stride+i+1] {
 389                                 t.Errorf("%v: wr[%v] != H[%v,%v]", prefix, i, i+1, i+1)
 390                         }
 391                         im := math.Sqrt(math.Abs(h.Data[(i+1)*h.Stride+i])) * math.Sqrt(math.Abs(h.Data[i*h.Stride+i+1]))
 392                         if math.Abs(im-wi[i]) > tol {
 393                                 t.Errorf("%v: unexpected value of wi[%v]: want %v, got %v", prefix, i, im, wi[i])
 394                         }
 395                 }
 396                 i += 2
 397         }
 398         // If the number of found eigenvalues is odd, at least one must be real.
 399         if (ihi+1-start)%2 != 0 && !hasReal {
 400                 t.Errorf("%v: expected at least one real eigenvalue", prefix)
 401         }
 402
 403         // Compare found eigenvalues to the reference, if known.
 404         if test.evWant != nil {
 405                 for i := start; i <= ihi; i++ {
 406                         ev := complex(wr[i], wi[i])
 407                         found, _ := containsComplex(test.evWant, ev, tol)
 408                         if !found {
 409                                 t.Errorf("%v: unexpected eigenvalue %v", prefix, ev)
 410                         }
 411                 }
 412         }
 413
 414         if !wantz {
 415                 return
 416         }
 417
 418         // Z should contain the orthogonal matrix U.
 419         if !isOrthonormal(z) {
 420                 t.Errorf("%v: Z is not orthogonal", prefix)
 421         }
 422         // Z should have been modified only in the
 423         // [iloz:ihiz+1,ilo:ihi+1] block.
 424         for i := 0; i < n; i++ {
 425                 for j := 0; j < n; j++ {
 426                         if iloz <= i && i <= ihiz && ilo <= j && j <= ihi {
 427                                 continue
 428                         }
 429                         if z.Data[i*z.Stride+j] != zCopy.Data[i*zCopy.Stride+j] {
 430                                 t.Errorf("%v: Z modified outside of [iloz:ihiz+1,ilo:ihi+1] block", prefix)
 431                         }
 432                 }
 433         }
 434         if wantt {
 435                 // Zero out h under the subdiagonal because Dlaqr04 uses it as
 436                 // workspace.
 437                 for i := 2; i < n; i++ {
 438                         for j := 0; j < i-1; j++ {
 439                                 h.Data[i*h.Stride+j] = 0
 440                         }
 441                 }
 442                 hz := eye(n, n)
 443                 blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, test.h, z, 0, hz)
 444                 zhz := eye(n, n)
 445                 blas64.Gemm(blas.Trans, blas.NoTrans, 1, z, hz, 0, zhz)
 446                 if !equalApproxGeneral(zhz, h, 10*tol) {
 447                         t.Errorf("%v: Z^T*(initial H)*Z and (final H) are not equal", prefix)
 448                 }
 449         }
 450 }