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