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