vendor/gonum.org/v1/gonum/lapack/testlapack/dgeev.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         "math/cmplx"
  11         "strconv"
  12         "testing"
  13
  14         "golang.org/x/exp/rand"
  15
  16         "gonum.org/v1/gonum/blas"
  17         "gonum.org/v1/gonum/blas/blas64"
  18         "gonum.org/v1/gonum/floats"
  19         "gonum.org/v1/gonum/lapack"
  20 )
  21
  22 type Dgeever interface {
  23         Dgeev(jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, n int, a []float64, lda int,
  24                 wr, wi []float64, vl []float64, ldvl int, vr []float64, ldvr int, work []float64, lwork int) int
  25 }
  26
  27 type dgeevTest struct {
  28         a      blas64.General
  29         evWant []complex128 // If nil, the eigenvalues are not known.
  30         valTol float64      // Tolerance for eigenvalue checks.
  31         vecTol float64      // Tolerance for eigenvector checks.
  32 }
  33
  34 func DgeevTest(t *testing.T, impl Dgeever) {
  35         rnd := rand.New(rand.NewSource(1))
  36
  37         for i, test := range []dgeevTest{
  38                 {
  39                         a:      A123{}.Matrix(),
  40                         evWant: A123{}.Eigenvalues(),
  41                 },
  42
  43                 dgeevTestForAntisymRandom(10, rnd),
  44                 dgeevTestForAntisymRandom(11, rnd),
  45                 dgeevTestForAntisymRandom(50, rnd),
  46                 dgeevTestForAntisymRandom(51, rnd),
  47                 dgeevTestForAntisymRandom(100, rnd),
  48                 dgeevTestForAntisymRandom(101, rnd),
  49
  50                 {
  51                         a:      Circulant(2).Matrix(),
  52                         evWant: Circulant(2).Eigenvalues(),
  53                 },
  54                 {
  55                         a:      Circulant(3).Matrix(),
  56                         evWant: Circulant(3).Eigenvalues(),
  57                 },
  58                 {
  59                         a:      Circulant(4).Matrix(),
  60                         evWant: Circulant(4).Eigenvalues(),
  61                 },
  62                 {
  63                         a:      Circulant(5).Matrix(),
  64                         evWant: Circulant(5).Eigenvalues(),
  65                 },
  66                 {
  67                         a:      Circulant(10).Matrix(),
  68                         evWant: Circulant(10).Eigenvalues(),
  69                 },
  70                 {
  71                         a:      Circulant(15).Matrix(),
  72                         evWant: Circulant(15).Eigenvalues(),
  73                         valTol: 1e-12,
  74                 },
  75                 {
  76                         a:      Circulant(30).Matrix(),
  77                         evWant: Circulant(30).Eigenvalues(),
  78                         valTol: 1e-11,
  79                         vecTol: 1e-12,
  80                 },
  81                 {
  82                         a:      Circulant(50).Matrix(),
  83                         evWant: Circulant(50).Eigenvalues(),
  84                         valTol: 1e-11,
  85                         vecTol: 1e-12,
  86                 },
  87                 {
  88                         a:      Circulant(101).Matrix(),
  89                         evWant: Circulant(101).Eigenvalues(),
  90                         valTol: 1e-10,
  91                         vecTol: 1e-11,
  92                 },
  93                 {
  94                         a:      Circulant(150).Matrix(),
  95                         evWant: Circulant(150).Eigenvalues(),
  96                         valTol: 1e-9,
  97                         vecTol: 1e-10,
  98                 },
  99
 100                 {
 101                         a:      Clement(2).Matrix(),
 102                         evWant: Clement(2).Eigenvalues(),
 103                 },
 104                 {
 105                         a:      Clement(3).Matrix(),
 106                         evWant: Clement(3).Eigenvalues(),
 107                 },
 108                 {
 109                         a:      Clement(4).Matrix(),
 110                         evWant: Clement(4).Eigenvalues(),
 111                 },
 112                 {
 113                         a:      Clement(5).Matrix(),
 114                         evWant: Clement(5).Eigenvalues(),
 115                 },
 116                 {
 117                         a:      Clement(10).Matrix(),
 118                         evWant: Clement(10).Eigenvalues(),
 119                 },
 120                 {
 121                         a:      Clement(15).Matrix(),
 122                         evWant: Clement(15).Eigenvalues(),
 123                 },
 124                 {
 125                         a:      Clement(30).Matrix(),
 126                         evWant: Clement(30).Eigenvalues(),
 127                         valTol: 1e-11,
 128                 },
 129                 {
 130                         a:      Clement(50).Matrix(),
 131                         evWant: Clement(50).Eigenvalues(),
 132                         valTol: 1e-7,
 133                         vecTol: 1e-11,
 134                 },
 135
 136                 {
 137                         a:      Creation(2).Matrix(),
 138                         evWant: Creation(2).Eigenvalues(),
 139                 },
 140                 {
 141                         a:      Creation(3).Matrix(),
 142                         evWant: Creation(3).Eigenvalues(),
 143                 },
 144                 {
 145                         a:      Creation(4).Matrix(),
 146                         evWant: Creation(4).Eigenvalues(),
 147                 },
 148                 {
 149                         a:      Creation(5).Matrix(),
 150                         evWant: Creation(5).Eigenvalues(),
 151                 },
 152                 {
 153                         a:      Creation(10).Matrix(),
 154                         evWant: Creation(10).Eigenvalues(),
 155                 },
 156                 {
 157                         a:      Creation(15).Matrix(),
 158                         evWant: Creation(15).Eigenvalues(),
 159                 },
 160                 {
 161                         a:      Creation(30).Matrix(),
 162                         evWant: Creation(30).Eigenvalues(),
 163                 },
 164                 {
 165                         a:      Creation(50).Matrix(),
 166                         evWant: Creation(50).Eigenvalues(),
 167                 },
 168                 {
 169                         a:      Creation(101).Matrix(),
 170                         evWant: Creation(101).Eigenvalues(),
 171                 },
 172                 {
 173                         a:      Creation(150).Matrix(),
 174                         evWant: Creation(150).Eigenvalues(),
 175                 },
 176
 177                 {
 178                         a:      Diagonal(0).Matrix(),
 179                         evWant: Diagonal(0).Eigenvalues(),
 180                 },
 181                 {
 182                         a:      Diagonal(10).Matrix(),
 183                         evWant: Diagonal(10).Eigenvalues(),
 184                 },
 185                 {
 186                         a:      Diagonal(50).Matrix(),
 187                         evWant: Diagonal(50).Eigenvalues(),
 188                 },
 189                 {
 190                         a:      Diagonal(151).Matrix(),
 191                         evWant: Diagonal(151).Eigenvalues(),
 192                 },
 193
 194                 {
 195                         a:      Downshift(2).Matrix(),
 196                         evWant: Downshift(2).Eigenvalues(),
 197                 },
 198                 {
 199                         a:      Downshift(3).Matrix(),
 200                         evWant: Downshift(3).Eigenvalues(),
 201                 },
 202                 {
 203                         a:      Downshift(4).Matrix(),
 204                         evWant: Downshift(4).Eigenvalues(),
 205                 },
 206                 {
 207                         a:      Downshift(5).Matrix(),
 208                         evWant: Downshift(5).Eigenvalues(),
 209                 },
 210                 {
 211                         a:      Downshift(10).Matrix(),
 212                         evWant: Downshift(10).Eigenvalues(),
 213                 },
 214                 {
 215                         a:      Downshift(15).Matrix(),
 216                         evWant: Downshift(15).Eigenvalues(),
 217                 },
 218                 {
 219                         a:      Downshift(30).Matrix(),
 220                         evWant: Downshift(30).Eigenvalues(),
 221                 },
 222                 {
 223                         a:      Downshift(50).Matrix(),
 224                         evWant: Downshift(50).Eigenvalues(),
 225                 },
 226                 {
 227                         a:      Downshift(101).Matrix(),
 228                         evWant: Downshift(101).Eigenvalues(),
 229                 },
 230                 {
 231                         a:      Downshift(150).Matrix(),
 232                         evWant: Downshift(150).Eigenvalues(),
 233                 },
 234
 235                 {
 236                         a:      Fibonacci(2).Matrix(),
 237                         evWant: Fibonacci(2).Eigenvalues(),
 238                 },
 239                 {
 240                         a:      Fibonacci(3).Matrix(),
 241                         evWant: Fibonacci(3).Eigenvalues(),
 242                 },
 243                 {
 244                         a:      Fibonacci(4).Matrix(),
 245                         evWant: Fibonacci(4).Eigenvalues(),
 246                 },
 247                 {
 248                         a:      Fibonacci(5).Matrix(),
 249                         evWant: Fibonacci(5).Eigenvalues(),
 250                 },
 251                 {
 252                         a:      Fibonacci(10).Matrix(),
 253                         evWant: Fibonacci(10).Eigenvalues(),
 254                 },
 255                 {
 256                         a:      Fibonacci(15).Matrix(),
 257                         evWant: Fibonacci(15).Eigenvalues(),
 258                 },
 259                 {
 260                         a:      Fibonacci(30).Matrix(),
 261                         evWant: Fibonacci(30).Eigenvalues(),
 262                 },
 263                 {
 264                         a:      Fibonacci(50).Matrix(),
 265                         evWant: Fibonacci(50).Eigenvalues(),
 266                 },
 267                 {
 268                         a:      Fibonacci(101).Matrix(),
 269                         evWant: Fibonacci(101).Eigenvalues(),
 270                 },
 271                 {
 272                         a:      Fibonacci(150).Matrix(),
 273                         evWant: Fibonacci(150).Eigenvalues(),
 274                 },
 275
 276                 {
 277                         a:      Gear(2).Matrix(),
 278                         evWant: Gear(2).Eigenvalues(),
 279                 },
 280                 {
 281                         a:      Gear(3).Matrix(),
 282                         evWant: Gear(3).Eigenvalues(),
 283                 },
 284                 {
 285                         a:      Gear(4).Matrix(),
 286                         evWant: Gear(4).Eigenvalues(),
 287                         valTol: 1e-7,
 288                 },
 289                 {
 290                         a:      Gear(5).Matrix(),
 291                         evWant: Gear(5).Eigenvalues(),
 292                 },
 293                 {
 294                         a:      Gear(10).Matrix(),
 295                         evWant: Gear(10).Eigenvalues(),
 296                         valTol: 1e-8,
 297                 },
 298                 {
 299                         a:      Gear(15).Matrix(),
 300                         evWant: Gear(15).Eigenvalues(),
 301                 },
 302                 {
 303                         a:      Gear(30).Matrix(),
 304                         evWant: Gear(30).Eigenvalues(),
 305                         valTol: 1e-8,
 306                 },
 307                 {
 308                         a:      Gear(50).Matrix(),
 309                         evWant: Gear(50).Eigenvalues(),
 310                         valTol: 1e-8,
 311                 },
 312                 {
 313                         a:      Gear(101).Matrix(),
 314                         evWant: Gear(101).Eigenvalues(),
 315                 },
 316                 {
 317                         a:      Gear(150).Matrix(),
 318                         evWant: Gear(150).Eigenvalues(),
 319                         valTol: 1e-8,
 320                 },
 321
 322                 {
 323                         a:      Grcar{N: 10, K: 3}.Matrix(),
 324                         evWant: Grcar{N: 10, K: 3}.Eigenvalues(),
 325                 },
 326                 {
 327                         a:      Grcar{N: 10, K: 7}.Matrix(),
 328                         evWant: Grcar{N: 10, K: 7}.Eigenvalues(),
 329                 },
 330                 {
 331                         a:      Grcar{N: 11, K: 7}.Matrix(),
 332                         evWant: Grcar{N: 11, K: 7}.Eigenvalues(),
 333                 },
 334                 {
 335                         a:      Grcar{N: 50, K: 3}.Matrix(),
 336                         evWant: Grcar{N: 50, K: 3}.Eigenvalues(),
 337                 },
 338                 {
 339                         a:      Grcar{N: 51, K: 3}.Matrix(),
 340                         evWant: Grcar{N: 51, K: 3}.Eigenvalues(),
 341                 },
 342                 {
 343                         a:      Grcar{N: 50, K: 10}.Matrix(),
 344                         evWant: Grcar{N: 50, K: 10}.Eigenvalues(),
 345                 },
 346                 {
 347                         a:      Grcar{N: 51, K: 10}.Matrix(),
 348                         evWant: Grcar{N: 51, K: 10}.Eigenvalues(),
 349                 },
 350                 {
 351                         a:      Grcar{N: 50, K: 30}.Matrix(),
 352                         evWant: Grcar{N: 50, K: 30}.Eigenvalues(),
 353                 },
 354                 {
 355                         a:      Grcar{N: 150, K: 2}.Matrix(),
 356                         evWant: Grcar{N: 150, K: 2}.Eigenvalues(),
 357                 },
 358                 {
 359                         a:      Grcar{N: 150, K: 148}.Matrix(),
 360                         evWant: Grcar{N: 150, K: 148}.Eigenvalues(),
 361                 },
 362
 363                 {
 364                         a:      Hanowa{N: 6, Alpha: 17}.Matrix(),
 365                         evWant: Hanowa{N: 6, Alpha: 17}.Eigenvalues(),
 366                 },
 367                 {
 368                         a:      Hanowa{N: 50, Alpha: -1}.Matrix(),
 369                         evWant: Hanowa{N: 50, Alpha: -1}.Eigenvalues(),
 370                 },
 371                 {
 372                         a:      Hanowa{N: 100, Alpha: -1}.Matrix(),
 373                         evWant: Hanowa{N: 100, Alpha: -1}.Eigenvalues(),
 374                 },
 375
 376                 {
 377                         a:      Lesp(2).Matrix(),
 378                         evWant: Lesp(2).Eigenvalues(),
 379                 },
 380                 {
 381                         a:      Lesp(3).Matrix(),
 382                         evWant: Lesp(3).Eigenvalues(),
 383                 },
 384                 {
 385                         a:      Lesp(4).Matrix(),
 386                         evWant: Lesp(4).Eigenvalues(),
 387                 },
 388                 {
 389                         a:      Lesp(5).Matrix(),
 390                         evWant: Lesp(5).Eigenvalues(),
 391                 },
 392                 {
 393                         a:      Lesp(10).Matrix(),
 394                         evWant: Lesp(10).Eigenvalues(),
 395                 },
 396                 {
 397                         a:      Lesp(15).Matrix(),
 398                         evWant: Lesp(15).Eigenvalues(),
 399                 },
 400                 {
 401                         a:      Lesp(30).Matrix(),
 402                         evWant: Lesp(30).Eigenvalues(),
 403                 },
 404                 {
 405                         a:      Lesp(50).Matrix(),
 406                         evWant: Lesp(50).Eigenvalues(),
 407                         valTol: 1e-12,
 408                         vecTol: 1e-12,
 409                 },
 410                 {
 411                         a:      Lesp(101).Matrix(),
 412                         evWant: Lesp(101).Eigenvalues(),
 413                         valTol: 1e-12,
 414                         vecTol: 1e-12,
 415                 },
 416                 {
 417                         a:      Lesp(150).Matrix(),
 418                         evWant: Lesp(150).Eigenvalues(),
 419                         valTol: 1e-12,
 420                         vecTol: 1e-12,
 421                 },
 422
 423                 {
 424                         a:      Rutis{}.Matrix(),
 425                         evWant: Rutis{}.Eigenvalues(),
 426                 },
 427
 428                 {
 429                         a:      Tris{N: 74, X: 1, Y: -2, Z: 1}.Matrix(),
 430                         evWant: Tris{N: 74, X: 1, Y: -2, Z: 1}.Eigenvalues(),
 431                 },
 432                 {
 433                         a:      Tris{N: 74, X: 1, Y: 2, Z: -3}.Matrix(),
 434                         evWant: Tris{N: 74, X: 1, Y: 2, Z: -3}.Eigenvalues(),
 435                 },
 436                 {
 437                         a:      Tris{N: 75, X: 1, Y: 2, Z: -3}.Matrix(),
 438                         evWant: Tris{N: 75, X: 1, Y: 2, Z: -3}.Eigenvalues(),
 439                 },
 440
 441                 {
 442                         a:      Wilk4{}.Matrix(),
 443                         evWant: Wilk4{}.Eigenvalues(),
 444                 },
 445                 {
 446                         a:      Wilk12{}.Matrix(),
 447                         evWant: Wilk12{}.Eigenvalues(),
 448                         valTol: 1e-7,
 449                 },
 450                 {
 451                         a:      Wilk20(0).Matrix(),
 452                         evWant: Wilk20(0).Eigenvalues(),
 453                 },
 454                 {
 455                         a:      Wilk20(1e-10).Matrix(),
 456                         evWant: Wilk20(1e-10).Eigenvalues(),
 457                         valTol: 1e-12,
 458                         vecTol: 1e-12,
 459                 },
 460
 461                 {
 462                         a:      Zero(1).Matrix(),
 463                         evWant: Zero(1).Eigenvalues(),
 464                 },
 465                 {
 466                         a:      Zero(10).Matrix(),
 467                         evWant: Zero(10).Eigenvalues(),
 468                 },
 469                 {
 470                         a:      Zero(50).Matrix(),
 471                         evWant: Zero(50).Eigenvalues(),
 472                 },
 473                 {
 474                         a:      Zero(100).Matrix(),
 475                         evWant: Zero(100).Eigenvalues(),
 476                 },
 477         } {
 478                 for _, jobvl := range []lapack.LeftEVJob{lapack.ComputeLeftEV, lapack.None} {
 479                         for _, jobvr := range []lapack.RightEVJob{lapack.ComputeRightEV, lapack.None} {
 480                                 for _, extra := range []int{0, 11} {
 481                                         for _, wl := range []worklen{minimumWork, mediumWork, optimumWork} {
 482                                                 testDgeev(t, impl, strconv.Itoa(i), test, jobvl, jobvr, extra, wl)
 483                                         }
 484                                 }
 485                         }
 486                 }
 487         }
 488
 489         for _, n := range []int{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 20, 50, 51, 100, 101} {
 490                 for _, jobvl := range []lapack.LeftEVJob{lapack.ComputeLeftEV, lapack.None} {
 491                         for _, jobvr := range []lapack.RightEVJob{lapack.ComputeRightEV, lapack.None} {
 492                                 for cas := 0; cas < 10; cas++ {
 493                                         // Create a block diagonal matrix with
 494                                         // random eigenvalues of random multiplicity.
 495                                         ev := make([]complex128, n)
 496                                         tmat := zeros(n, n, n)
 497                                         for i := 0; i < n; {
 498                                                 re := rnd.NormFloat64()
 499                                                 if i == n-1 || rnd.Float64() < 0.5 {
 500                                                         // Real eigenvalue.
 501                                                         nb := rnd.Intn(min(4, n-i)) + 1
 502                                                         for k := 0; k < nb; k++ {
 503                                                                 tmat.Data[i*tmat.Stride+i] = re
 504                                                                 ev[i] = complex(re, 0)
 505                                                                 i++
 506                                                         }
 507                                                         continue
 508                                                 }
 509                                                 // Complex eigenvalue.
 510                                                 im := rnd.NormFloat64()
 511                                                 nb := rnd.Intn(min(4, (n-i)/2)) + 1
 512                                                 for k := 0; k < nb; k++ {
 513                                                         // 2×2 block for the complex eigenvalue.
 514                                                         tmat.Data[i*tmat.Stride+i] = re
 515                                                         tmat.Data[(i+1)*tmat.Stride+i+1] = re
 516                                                         tmat.Data[(i+1)*tmat.Stride+i] = -im
 517                                                         tmat.Data[i*tmat.Stride+i+1] = im
 518                                                         ev[i] = complex(re, im)
 519                                                         ev[i+1] = complex(re, -im)
 520                                                         i += 2
 521                                                 }
 522                                         }
 523
 524                                         // Compute A = Q T Q^T where Q is an
 525                                         // orthogonal matrix.
 526                                         q := randomOrthogonal(n, rnd)
 527                                         tq := zeros(n, n, n)
 528                                         blas64.Gemm(blas.NoTrans, blas.Trans, 1, tmat, q, 0, tq)
 529                                         a := zeros(n, n, n)
 530                                         blas64.Gemm(blas.NoTrans, blas.NoTrans, 1, q, tq, 0, a)
 531
 532                                         test := dgeevTest{
 533                                                 a:      a,
 534                                                 evWant: ev,
 535                                                 valTol: 1e-12,
 536                                                 vecTol: 1e-7,
 537                                         }
 538                                         testDgeev(t, impl, "random", test, jobvl, jobvr, 0, optimumWork)
 539                                 }
 540                         }
 541                 }
 542         }
 543 }
 544
 545 func testDgeev(t *testing.T, impl Dgeever, tc string, test dgeevTest, jobvl lapack.LeftEVJob, jobvr lapack.RightEVJob, extra int, wl worklen) {
 546         const defaultTol = 1e-12
 547         valTol := test.valTol
 548         if valTol == 0 {
 549                 valTol = defaultTol
 550         }
 551         vecTol := test.vecTol
 552         if vecTol == 0 {
 553                 vecTol = defaultTol
 554         }
 555
 556         a := cloneGeneral(test.a)
 557         n := a.Rows
 558
 559         var vl blas64.General
 560         if jobvl == lapack.ComputeLeftEV {
 561                 vl = nanGeneral(n, n, n)
 562         }
 563
 564         var vr blas64.General
 565         if jobvr == lapack.ComputeRightEV {
 566                 vr = nanGeneral(n, n, n)
 567         }
 568
 569         wr := make([]float64, n)
 570         wi := make([]float64, n)
 571
 572         var lwork int
 573         switch wl {
 574         case minimumWork:
 575                 if jobvl == lapack.ComputeLeftEV || jobvr == lapack.ComputeRightEV {
 576                         lwork = max(1, 4*n)
 577                 } else {
 578                         lwork = max(1, 3*n)
 579                 }
 580         case mediumWork:
 581                 work := make([]float64, 1)
 582                 impl.Dgeev(jobvl, jobvr, n, nil, 1, nil, nil, nil, 1, nil, 1, work, -1)
 583                 if jobvl == lapack.ComputeLeftEV || jobvr == lapack.ComputeRightEV {
 584                         lwork = (int(work[0]) + 4*n) / 2
 585                 } else {
 586                         lwork = (int(work[0]) + 3*n) / 2
 587                 }
 588                 lwork = max(1, lwork)
 589         case optimumWork:
 590                 work := make([]float64, 1)
 591                 impl.Dgeev(jobvl, jobvr, n, nil, 1, nil, nil, nil, 1, nil, 1, work, -1)
 592                 lwork = int(work[0])
 593         }
 594         work := make([]float64, lwork)
 595
 596         first := impl.Dgeev(jobvl, jobvr, n, a.Data, a.Stride, wr, wi,
 597                 vl.Data, vl.Stride, vr.Data, vr.Stride, work, len(work))
 598
 599         prefix := fmt.Sprintf("Case #%v: n=%v, jobvl=%v, jobvr=%v, extra=%v, work=%v",
 600                 tc, n, jobvl, jobvr, extra, wl)
 601
 602         if !generalOutsideAllNaN(vl) {
 603                 t.Errorf("%v: out-of-range write to VL", prefix)
 604         }
 605         if !generalOutsideAllNaN(vr) {
 606                 t.Errorf("%v: out-of-range write to VR", prefix)
 607         }
 608
 609         if first > 0 {
 610                 t.Logf("%v: all eigenvalues haven't been computed, first=%v", prefix, first)
 611         }
 612
 613         // Check that conjugate pair eigevalues are ordered correctly.
 614         for i := first; i < n; {
 615                 if wi[i] == 0 {
 616                         i++
 617                         continue
 618                 }
 619                 if wr[i] != wr[i+1] {
 620                         t.Errorf("%v: real parts of %vth conjugate pair not equal", prefix, i)
 621                 }
 622                 if wi[i] < 0 || wi[i+1] > 0 {
 623                         t.Errorf("%v: unexpected ordering of %vth conjugate pair", prefix, i)
 624                 }
 625                 i += 2
 626         }
 627
 628         // Check the computed eigenvalues against provided known eigenvalues.
 629         if test.evWant != nil {
 630                 used := make([]bool, n)
 631                 for i := first; i < n; i++ {
 632                         evGot := complex(wr[i], wi[i])
 633                         idx := -1
 634                         for k, evWant := range test.evWant {
 635                                 if !used[k] && cmplx.Abs(evWant-evGot) < valTol {
 636                                         idx = k
 637                                         used[k] = true
 638                                         break
 639                                 }
 640                         }
 641                         if idx == -1 {
 642                                 t.Errorf("%v: unexpected eigenvalue %v", prefix, evGot)
 643                         }
 644                 }
 645         }
 646
 647         if first > 0 || (jobvl == lapack.None && jobvr == lapack.None) {
 648                 // No eigenvectors have been computed.
 649                 return
 650         }
 651
 652         // Check that the columns of VL and VR are eigenvectors that correspond
 653         // to the computed eigenvalues.
 654         for k := 0; k < n; {
 655                 if wi[k] == 0 {
 656                         if jobvl == lapack.ComputeLeftEV {
 657                                 ev := columnOf(vl, k)
 658                                 if !isLeftEigenvectorOf(test.a, ev, nil, complex(wr[k], 0), vecTol) {
 659                                         t.Errorf("%v: VL[:,%v] is not left real eigenvector",
 660                                                 prefix, k)
 661                                 }
 662
 663                                 norm := floats.Norm(ev, 2)
 664                                 if math.Abs(norm-1) >= defaultTol {
 665                                         t.Errorf("%v: norm of left real eigenvector %v not equal to 1: got %v",
 666                                                 prefix, k, norm)
 667                                 }
 668                         }
 669                         if jobvr == lapack.ComputeRightEV {
 670                                 ev := columnOf(vr, k)
 671                                 if !isRightEigenvectorOf(test.a, ev, nil, complex(wr[k], 0), vecTol) {
 672                                         t.Errorf("%v: VR[:,%v] is not right real eigenvector",
 673                                                 prefix, k)
 674                                 }
 675
 676                                 norm := floats.Norm(ev, 2)
 677                                 if math.Abs(norm-1) >= defaultTol {
 678                                         t.Errorf("%v: norm of right real eigenvector %v not equal to 1: got %v",
 679                                                 prefix, k, norm)
 680                                 }
 681                         }
 682                         k++
 683                 } else {
 684                         if jobvl == lapack.ComputeLeftEV {
 685                                 evre := columnOf(vl, k)
 686                                 evim := columnOf(vl, k+1)
 687                                 if !isLeftEigenvectorOf(test.a, evre, evim, complex(wr[k], wi[k]), vecTol) {
 688                                         t.Errorf("%v: VL[:,%v:%v] is not left complex eigenvector",
 689                                                 prefix, k, k+1)
 690                                 }
 691                                 floats.Scale(-1, evim)
 692                                 if !isLeftEigenvectorOf(test.a, evre, evim, complex(wr[k+1], wi[k+1]), vecTol) {
 693                                         t.Errorf("%v: VL[:,%v:%v] is not left complex eigenvector",
 694                                                 prefix, k, k+1)
 695                                 }
 696
 697                                 norm := math.Hypot(floats.Norm(evre, 2), floats.Norm(evim, 2))
 698                                 if math.Abs(norm-1) > defaultTol {
 699                                         t.Errorf("%v: norm of left complex eigenvector %v not equal to 1: got %v",
 700                                                 prefix, k, norm)
 701                                 }
 702                         }
 703                         if jobvr == lapack.ComputeRightEV {
 704                                 evre := columnOf(vr, k)
 705                                 evim := columnOf(vr, k+1)
 706                                 if !isRightEigenvectorOf(test.a, evre, evim, complex(wr[k], wi[k]), vecTol) {
 707                                         t.Errorf("%v: VR[:,%v:%v] is not right complex eigenvector",
 708                                                 prefix, k, k+1)
 709                                 }
 710                                 floats.Scale(-1, evim)
 711                                 if !isRightEigenvectorOf(test.a, evre, evim, complex(wr[k+1], wi[k+1]), vecTol) {
 712                                         t.Errorf("%v: VR[:,%v:%v] is not right complex eigenvector",
 713                                                 prefix, k, k+1)
 714                                 }
 715
 716                                 norm := math.Hypot(floats.Norm(evre, 2), floats.Norm(evim, 2))
 717                                 if math.Abs(norm-1) > defaultTol {
 718                                         t.Errorf("%v: norm of right complex eigenvector %v not equal to 1: got %v",
 719                                                 prefix, k, norm)
 720                                 }
 721                         }
 722                         // We don't test whether the largest component is real
 723                         // because checking it is flaky due to rounding errors.
 724
 725                         k += 2
 726                 }
 727         }
 728 }
 729
 730 func dgeevTestForAntisymRandom(n int, rnd *rand.Rand) dgeevTest {
 731         a := NewAntisymRandom(n, rnd)
 732         return dgeevTest{
 733                 a:      a.Matrix(),
 734                 evWant: a.Eigenvalues(),
 735         }
 736 }