001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.jama; 022 023/** Eigenvalues and eigenvectors of a real matrix. 024<P> 025 If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is 026 diagonal and the eigenvector matrix V is orthogonal. 027 I.e. A = V.times(D.times(V.transpose())) and 028 V.times(V.transpose()) equals the identity matrix. 029<P> 030 If A is not symmetric, then the eigenvalue matrix D is block diagonal 031 with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, 032 lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The 033 columns of V represent the eigenvectors in the sense that A*V = V*D, 034 i.e. A.times(V) equals V.times(D). The matrix V may be badly 035 conditioned, or even singular, so the validity of the equation 036 A = V*D*inverse(V) depends upon V.cond(). 037 */ 038 039public class EigenvalueDecomposition implements java.io.Serializable { 040 041 /** 042 * 043 */ 044 private static final long serialVersionUID = 3557806515310435894L; 045 046 047/* ------------------------ 048 Class variables 049 * ------------------------ */ 050 051/** Row and column dimension (square matrix). 052 @serial matrix dimension. 053 */ 054 private int n; 055 056 /** Symmetry flag. 057 @serial internal symmetry flag. 058 */ 059 private boolean issymmetric; 060 061 /** Arrays for internal storage of eigenvalues. 062 @serial internal storage of eigenvalues. 063 */ 064 private double[] d, e; 065 066 /** Array for internal storage of eigenvectors. 067 @serial internal storage of eigenvectors. 068 */ 069 private double[][] V; 070 071 /** Array for internal storage of nonsymmetric Hessenberg form. 072 @serial internal storage of nonsymmetric Hessenberg form. 073 */ 074 private double[][] H; 075 076 /** Working storage for nonsymmetric algorithm. 077 @serial working storage for nonsymmetric algorithm. 078 */ 079 private double[] ort; 080 081/* ------------------------ 082 Private Methods 083 * ------------------------ */ 084 085 // Symmetric Householder reduction to tridiagonal form. 086 087 private void tred2 () { 088 089 // This is derived from the Algol procedures tred2 by 090 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 091 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 092 // Fortran subroutine in EISPACK. 093 094 for (int j = 0; j < n; j++) { 095 d[j] = V[n-1][j]; 096 } 097 098 // Householder reduction to tridiagonal form. 099 100 for (int i = n-1; i > 0; i--) { 101 102 // Scale to avoid under/overflow. 103 104 double scale = 0.0; 105 double h = 0.0; 106 for (int k = 0; k < i; k++) { 107 scale = scale + Math.abs(d[k]); 108 } 109 if (scale == 0.0) { 110 e[i] = d[i-1]; 111 for (int j = 0; j < i; j++) { 112 d[j] = V[i-1][j]; 113 V[i][j] = 0.0; 114 V[j][i] = 0.0; 115 } 116 } else { 117 118 // Generate Householder vector. 119 120 for (int k = 0; k < i; k++) { 121 d[k] /= scale; 122 h += d[k] * d[k]; 123 } 124 double f = d[i-1]; 125 double g = Math.sqrt(h); 126 if (f > 0) { 127 g = -g; 128 } 129 e[i] = scale * g; 130 h = h - f * g; 131 d[i-1] = f - g; 132 for (int j = 0; j < i; j++) { 133 e[j] = 0.0; 134 } 135 136 // Apply similarity transformation to remaining columns. 137 138 for (int j = 0; j < i; j++) { 139 f = d[j]; 140 V[j][i] = f; 141 g = e[j] + V[j][j] * f; 142 for (int k = j+1; k <= i-1; k++) { 143 g += V[k][j] * d[k]; 144 e[k] += V[k][j] * f; 145 } 146 e[j] = g; 147 } 148 f = 0.0; 149 for (int j = 0; j < i; j++) { 150 e[j] /= h; 151 f += e[j] * d[j]; 152 } 153 double hh = f / (h + h); 154 for (int j = 0; j < i; j++) { 155 e[j] -= hh * d[j]; 156 } 157 for (int j = 0; j < i; j++) { 158 f = d[j]; 159 g = e[j]; 160 for (int k = j; k <= i-1; k++) { 161 V[k][j] -= (f * e[k] + g * d[k]); 162 } 163 d[j] = V[i-1][j]; 164 V[i][j] = 0.0; 165 } 166 } 167 d[i] = h; 168 } 169 170 // Accumulate transformations. 171 172 for (int i = 0; i < n-1; i++) { 173 V[n-1][i] = V[i][i]; 174 V[i][i] = 1.0; 175 double h = d[i+1]; 176 if (h != 0.0) { 177 for (int k = 0; k <= i; k++) { 178 d[k] = V[k][i+1] / h; 179 } 180 for (int j = 0; j <= i; j++) { 181 double g = 0.0; 182 for (int k = 0; k <= i; k++) { 183 g += V[k][i+1] * V[k][j]; 184 } 185 for (int k = 0; k <= i; k++) { 186 V[k][j] -= g * d[k]; 187 } 188 } 189 } 190 for (int k = 0; k <= i; k++) { 191 V[k][i+1] = 0.0; 192 } 193 } 194 for (int j = 0; j < n; j++) { 195 d[j] = V[n-1][j]; 196 V[n-1][j] = 0.0; 197 } 198 V[n-1][n-1] = 1.0; 199 e[0] = 0.0; 200 } 201 202 // Symmetric tridiagonal QL algorithm. 203 204 private void tql2 () { 205 206 // This is derived from the Algol procedures tql2, by 207 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 208 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 209 // Fortran subroutine in EISPACK. 210 211 for (int i = 1; i < n; i++) { 212 e[i-1] = e[i]; 213 } 214 e[n-1] = 0.0; 215 216 double f = 0.0; 217 double tst1 = 0.0; 218 double eps = Math.pow(2.0,-52.0); 219 for (int l = 0; l < n; l++) { 220 221 // Find small subdiagonal element 222 223 tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l])); 224 int m = l; 225 while (m < n) { 226 if (Math.abs(e[m]) <= eps*tst1) { 227 break; 228 } 229 m++; 230 } 231 232 // If m == l, d[l] is an eigenvalue, 233 // otherwise, iterate. 234 235 if (m > l) { 236 int iter = 0; 237 do { 238 iter = iter + 1; // (Could check iteration count here.) 239 240 // Compute implicit shift 241 242 double g = d[l]; 243 double p = (d[l+1] - g) / (2.0 * e[l]); 244 double r = Maths.hypot(p,1.0); 245 if (p < 0) { 246 r = -r; 247 } 248 d[l] = e[l] / (p + r); 249 d[l+1] = e[l] * (p + r); 250 double dl1 = d[l+1]; 251 double h = g - d[l]; 252 for (int i = l+2; i < n; i++) { 253 d[i] -= h; 254 } 255 f = f + h; 256 257 // Implicit QL transformation. 258 259 p = d[m]; 260 double c = 1.0; 261 double c2 = c; 262 double c3 = c; 263 double el1 = e[l+1]; 264 double s = 0.0; 265 double s2 = 0.0; 266 for (int i = m-1; i >= l; i--) { 267 c3 = c2; 268 c2 = c; 269 s2 = s; 270 g = c * e[i]; 271 h = c * p; 272 r = Maths.hypot(p,e[i]); 273 e[i+1] = s * r; 274 s = e[i] / r; 275 c = p / r; 276 p = c * d[i] - s * g; 277 d[i+1] = h + s * (c * g + s * d[i]); 278 279 // Accumulate transformation. 280 281 for (int k = 0; k < n; k++) { 282 h = V[k][i+1]; 283 V[k][i+1] = s * V[k][i] + c * h; 284 V[k][i] = c * V[k][i] - s * h; 285 } 286 } 287 p = -s * s2 * c3 * el1 * e[l] / dl1; 288 e[l] = s * p; 289 d[l] = c * p; 290 291 // Check for convergence. 292 293 } while (Math.abs(e[l]) > eps*tst1); 294 } 295 d[l] = d[l] + f; 296 e[l] = 0.0; 297 } 298 299 // Sort eigenvalues and corresponding vectors. 300 301 for (int i = 0; i < n-1; i++) { 302 int k = i; 303 double p = d[i]; 304 for (int j = i+1; j < n; j++) { 305 if (d[j] < p) { 306 k = j; 307 p = d[j]; 308 } 309 } 310 if (k != i) { 311 d[k] = d[i]; 312 d[i] = p; 313 for (int j = 0; j < n; j++) { 314 p = V[j][i]; 315 V[j][i] = V[j][k]; 316 V[j][k] = p; 317 } 318 } 319 } 320 } 321 322 // Nonsymmetric reduction to Hessenberg form. 323 324 private void orthes () { 325 326 // This is derived from the Algol procedures orthes and ortran, 327 // by Martin and Wilkinson, Handbook for Auto. Comp., 328 // Vol.ii-Linear Algebra, and the corresponding 329 // Fortran subroutines in EISPACK. 330 331 int low = 0; 332 int high = n-1; 333 334 for (int m = low+1; m <= high-1; m++) { 335 336 // Scale column. 337 338 double scale = 0.0; 339 for (int i = m; i <= high; i++) { 340 scale = scale + Math.abs(H[i][m-1]); 341 } 342 if (scale != 0.0) { 343 344 // Compute Householder transformation. 345 346 double h = 0.0; 347 for (int i = high; i >= m; i--) { 348 ort[i] = H[i][m-1]/scale; 349 h += ort[i] * ort[i]; 350 } 351 double g = Math.sqrt(h); 352 if (ort[m] > 0) { 353 g = -g; 354 } 355 h = h - ort[m] * g; 356 ort[m] = ort[m] - g; 357 358 // Apply Householder similarity transformation 359 // H = (I-u*u'/h)*H*(I-u*u')/h) 360 361 for (int j = m; j < n; j++) { 362 double f = 0.0; 363 for (int i = high; i >= m; i--) { 364 f += ort[i]*H[i][j]; 365 } 366 f = f/h; 367 for (int i = m; i <= high; i++) { 368 H[i][j] -= f*ort[i]; 369 } 370 } 371 372 for (int i = 0; i <= high; i++) { 373 double f = 0.0; 374 for (int j = high; j >= m; j--) { 375 f += ort[j]*H[i][j]; 376 } 377 f = f/h; 378 for (int j = m; j <= high; j++) { 379 H[i][j] -= f*ort[j]; 380 } 381 } 382 ort[m] = scale*ort[m]; 383 H[m][m-1] = scale*g; 384 } 385 } 386 387 // Accumulate transformations (Algol's ortran). 388 389 for (int i = 0; i < n; i++) { 390 for (int j = 0; j < n; j++) { 391 V[i][j] = (i == j ? 1.0 : 0.0); 392 } 393 } 394 395 for (int m = high-1; m >= low+1; m--) { 396 if (H[m][m-1] != 0.0) { 397 for (int i = m+1; i <= high; i++) { 398 ort[i] = H[i][m-1]; 399 } 400 for (int j = m; j <= high; j++) { 401 double g = 0.0; 402 for (int i = m; i <= high; i++) { 403 g += ort[i] * V[i][j]; 404 } 405 // Double division avoids possible underflow 406 g = (g / ort[m]) / H[m][m-1]; 407 for (int i = m; i <= high; i++) { 408 V[i][j] += g * ort[i]; 409 } 410 } 411 } 412 } 413 } 414 415 416 // Complex scalar division. 417 418 private transient double cdivr, cdivi; 419 private void cdiv(double xr, double xi, double yr, double yi) { 420 double r,d; 421 if (Math.abs(yr) > Math.abs(yi)) { 422 r = yi/yr; 423 d = yr + r*yi; 424 cdivr = (xr + r*xi)/d; 425 cdivi = (xi - r*xr)/d; 426 } else { 427 r = yr/yi; 428 d = yi + r*yr; 429 cdivr = (r*xr + xi)/d; 430 cdivi = (r*xi - xr)/d; 431 } 432 } 433 434 435 // Nonsymmetric reduction from Hessenberg to real Schur form. 436 437 private void hqr2 () { 438 439 // This is derived from the Algol procedure hqr2, 440 // by Martin and Wilkinson, Handbook for Auto. Comp., 441 // Vol.ii-Linear Algebra, and the corresponding 442 // Fortran subroutine in EISPACK. 443 444 // Initialize 445 446 int nn = this.n; 447 int n = nn-1; 448 int low = 0; 449 int high = nn-1; 450 double eps = Math.pow(2.0,-52.0); 451 double exshift = 0.0; 452 double p=0,q=0,r=0,s=0,z=0,t,w,x,y; 453 454 // Store roots isolated by balanc and compute matrix norm 455 456 double norm = 0.0; 457 for (int i = 0; i < nn; i++) { 458 if (i < low || i > high) { 459 d[i] = H[i][i]; 460 e[i] = 0.0; 461 } 462 for (int j = Math.max(i-1,0); j < nn; j++) { 463 norm = norm + Math.abs(H[i][j]); 464 } 465 } 466 467 // Outer loop over eigenvalue index 468 469 int iter = 0; 470 while (n >= low) { 471 472 // Look for single small sub-diagonal element 473 474 int l = n; 475 while (l > low) { 476 s = Math.abs(H[l-1][l-1]) + Math.abs(H[l][l]); 477 if (s == 0.0) { 478 s = norm; 479 } 480 if (Math.abs(H[l][l-1]) < eps * s) { 481 break; 482 } 483 l--; 484 } 485 486 // Check for convergence 487 // One root found 488 489 if (l == n) { 490 H[n][n] = H[n][n] + exshift; 491 d[n] = H[n][n]; 492 e[n] = 0.0; 493 n--; 494 iter = 0; 495 496 // Two roots found 497 498 } else if (l == n-1) { 499 w = H[n][n-1] * H[n-1][n]; 500 p = (H[n-1][n-1] - H[n][n]) / 2.0; 501 q = p * p + w; 502 z = Math.sqrt(Math.abs(q)); 503 H[n][n] = H[n][n] + exshift; 504 H[n-1][n-1] = H[n-1][n-1] + exshift; 505 x = H[n][n]; 506 507 // Real pair 508 509 if (q >= 0) { 510 if (p >= 0) { 511 z = p + z; 512 } else { 513 z = p - z; 514 } 515 d[n-1] = x + z; 516 d[n] = d[n-1]; 517 if (z != 0.0) { 518 d[n] = x - w / z; 519 } 520 e[n-1] = 0.0; 521 e[n] = 0.0; 522 x = H[n][n-1]; 523 s = Math.abs(x) + Math.abs(z); 524 p = x / s; 525 q = z / s; 526 r = Math.sqrt(p * p+q * q); 527 p = p / r; 528 q = q / r; 529 530 // Row modification 531 532 for (int j = n-1; j < nn; j++) { 533 z = H[n-1][j]; 534 H[n-1][j] = q * z + p * H[n][j]; 535 H[n][j] = q * H[n][j] - p * z; 536 } 537 538 // Column modification 539 540 for (int i = 0; i <= n; i++) { 541 z = H[i][n-1]; 542 H[i][n-1] = q * z + p * H[i][n]; 543 H[i][n] = q * H[i][n] - p * z; 544 } 545 546 // Accumulate transformations 547 548 for (int i = low; i <= high; i++) { 549 z = V[i][n-1]; 550 V[i][n-1] = q * z + p * V[i][n]; 551 V[i][n] = q * V[i][n] - p * z; 552 } 553 554 // Complex pair 555 556 } else { 557 d[n-1] = x + p; 558 d[n] = x + p; 559 e[n-1] = z; 560 e[n] = -z; 561 } 562 n = n - 2; 563 iter = 0; 564 565 // No convergence yet 566 567 } else { 568 569 // Form shift 570 571 x = H[n][n]; 572 y = 0.0; 573 w = 0.0; 574 if (l < n) { 575 y = H[n-1][n-1]; 576 w = H[n][n-1] * H[n-1][n]; 577 } 578 579 // Wilkinson's original ad hoc shift 580 581 if (iter == 10) { 582 exshift += x; 583 for (int i = low; i <= n; i++) { 584 H[i][i] -= x; 585 } 586 s = Math.abs(H[n][n-1]) + Math.abs(H[n-1][n-2]); 587 x = y = 0.75 * s; 588 w = -0.4375 * s * s; 589 } 590 591 // MATLAB's new ad hoc shift 592 593 if (iter == 30) { 594 s = (y - x) / 2.0; 595 s = s * s + w; 596 if (s > 0) { 597 s = Math.sqrt(s); 598 if (y < x) { 599 s = -s; 600 } 601 s = x - w / ((y - x) / 2.0 + s); 602 for (int i = low; i <= n; i++) { 603 H[i][i] -= s; 604 } 605 exshift += s; 606 x = y = w = 0.964; 607 } 608 } 609 610 iter = iter + 1; // (Could check iteration count here.) 611 612 // Look for two consecutive small sub-diagonal elements 613 614 int m = n-2; 615 while (m >= l) { 616 z = H[m][m]; 617 r = x - z; 618 s = y - z; 619 p = (r * s - w) / H[m+1][m] + H[m][m+1]; 620 q = H[m+1][m+1] - z - r - s; 621 r = H[m+2][m+1]; 622 s = Math.abs(p) + Math.abs(q) + Math.abs(r); 623 p = p / s; 624 q = q / s; 625 r = r / s; 626 if (m == l) { 627 break; 628 } 629 if (Math.abs(H[m][m-1]) * (Math.abs(q) + Math.abs(r)) < 630 eps * (Math.abs(p) * (Math.abs(H[m-1][m-1]) + Math.abs(z) + 631 Math.abs(H[m+1][m+1])))) { 632 break; 633 } 634 m--; 635 } 636 637 for (int i = m+2; i <= n; i++) { 638 H[i][i-2] = 0.0; 639 if (i > m+2) { 640 H[i][i-3] = 0.0; 641 } 642 } 643 644 // Double QR step involving rows l:n and columns m:n 645 646 647 for (int k = m; k <= n-1; k++) { 648 boolean notlast = (k != n-1); 649 if (k != m) { 650 p = H[k][k-1]; 651 q = H[k+1][k-1]; 652 r = (notlast ? H[k+2][k-1] : 0.0); 653 x = Math.abs(p) + Math.abs(q) + Math.abs(r); 654 if (x == 0.0) { 655 continue; 656 } 657 p = p / x; 658 q = q / x; 659 r = r / x; 660 } 661 662 s = Math.sqrt(p * p + q * q + r * r); 663 if (p < 0) { 664 s = -s; 665 } 666 if (s != 0) { 667 if (k != m) { 668 H[k][k-1] = -s * x; 669 } else if (l != m) { 670 H[k][k-1] = -H[k][k-1]; 671 } 672 p = p + s; 673 x = p / s; 674 y = q / s; 675 z = r / s; 676 q = q / p; 677 r = r / p; 678 679 // Row modification 680 681 for (int j = k; j < nn; j++) { 682 p = H[k][j] + q * H[k+1][j]; 683 if (notlast) { 684 p = p + r * H[k+2][j]; 685 H[k+2][j] = H[k+2][j] - p * z; 686 } 687 H[k][j] = H[k][j] - p * x; 688 H[k+1][j] = H[k+1][j] - p * y; 689 } 690 691 // Column modification 692 693 for (int i = 0; i <= Math.min(n,k+3); i++) { 694 p = x * H[i][k] + y * H[i][k+1]; 695 if (notlast) { 696 p = p + z * H[i][k+2]; 697 H[i][k+2] = H[i][k+2] - p * r; 698 } 699 H[i][k] = H[i][k] - p; 700 H[i][k+1] = H[i][k+1] - p * q; 701 } 702 703 // Accumulate transformations 704 705 for (int i = low; i <= high; i++) { 706 p = x * V[i][k] + y * V[i][k+1]; 707 if (notlast) { 708 p = p + z * V[i][k+2]; 709 V[i][k+2] = V[i][k+2] - p * r; 710 } 711 V[i][k] = V[i][k] - p; 712 V[i][k+1] = V[i][k+1] - p * q; 713 } 714 } // (s != 0) 715 } // k loop 716 } // check convergence 717 } // while (n >= low) 718 719 // Backsubstitute to find vectors of upper triangular form 720 721 if (norm == 0.0) { 722 return; 723 } 724 725 for (n = nn-1; n >= 0; n--) { 726 p = d[n]; 727 q = e[n]; 728 729 // Real vector 730 731 if (q == 0) { 732 int l = n; 733 H[n][n] = 1.0; 734 for (int i = n-1; i >= 0; i--) { 735 w = H[i][i] - p; 736 r = 0.0; 737 for (int j = l; j <= n; j++) { 738 r = r + H[i][j] * H[j][n]; 739 } 740 if (e[i] < 0.0) { 741 z = w; 742 s = r; 743 } else { 744 l = i; 745 if (e[i] == 0.0) { 746 if (w != 0.0) { 747 H[i][n] = -r / w; 748 } else { 749 H[i][n] = -r / (eps * norm); 750 } 751 752 // Solve real equations 753 754 } else { 755 x = H[i][i+1]; 756 y = H[i+1][i]; 757 q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; 758 t = (x * s - z * r) / q; 759 H[i][n] = t; 760 if (Math.abs(x) > Math.abs(z)) { 761 H[i+1][n] = (-r - w * t) / x; 762 } else { 763 H[i+1][n] = (-s - y * t) / z; 764 } 765 } 766 767 // Overflow control 768 769 t = Math.abs(H[i][n]); 770 if ((eps * t) * t > 1) { 771 for (int j = i; j <= n; j++) { 772 H[j][n] = H[j][n] / t; 773 } 774 } 775 } 776 } 777 778 // Complex vector 779 780 } else if (q < 0) { 781 int l = n-1; 782 783 // Last vector component imaginary so matrix is triangular 784 785 if (Math.abs(H[n][n-1]) > Math.abs(H[n-1][n])) { 786 H[n-1][n-1] = q / H[n][n-1]; 787 H[n-1][n] = -(H[n][n] - p) / H[n][n-1]; 788 } else { 789 cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q); 790 H[n-1][n-1] = cdivr; 791 H[n-1][n] = cdivi; 792 } 793 H[n][n-1] = 0.0; 794 H[n][n] = 1.0; 795 for (int i = n-2; i >= 0; i--) { 796 double ra,sa,vr,vi; 797 ra = 0.0; 798 sa = 0.0; 799 for (int j = l; j <= n; j++) { 800 ra = ra + H[i][j] * H[j][n-1]; 801 sa = sa + H[i][j] * H[j][n]; 802 } 803 w = H[i][i] - p; 804 805 if (e[i] < 0.0) { 806 z = w; 807 r = ra; 808 s = sa; 809 } else { 810 l = i; 811 if (e[i] == 0) { 812 cdiv(-ra,-sa,w,q); 813 H[i][n-1] = cdivr; 814 H[i][n] = cdivi; 815 } else { 816 817 // Solve complex equations 818 819 x = H[i][i+1]; 820 y = H[i+1][i]; 821 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; 822 vi = (d[i] - p) * 2.0 * q; 823 if (vr == 0.0 && vi == 0.0) { 824 vr = eps * norm * (Math.abs(w) + Math.abs(q) + 825 Math.abs(x) + Math.abs(y) + Math.abs(z)); 826 } 827 cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi); 828 H[i][n-1] = cdivr; 829 H[i][n] = cdivi; 830 if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) { 831 H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x; 832 H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x; 833 } else { 834 cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q); 835 H[i+1][n-1] = cdivr; 836 H[i+1][n] = cdivi; 837 } 838 } 839 840 // Overflow control 841 842 t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n])); 843 if ((eps * t) * t > 1) { 844 for (int j = i; j <= n; j++) { 845 H[j][n-1] = H[j][n-1] / t; 846 H[j][n] = H[j][n] / t; 847 } 848 } 849 } 850 } 851 } 852 } 853 854 // Vectors of isolated roots 855 856 for (int i = 0; i < nn; i++) { 857 if (i < low || i > high) { 858 for (int j = i; j < nn; j++) { 859 V[i][j] = H[i][j]; 860 } 861 } 862 } 863 864 // Back transformation to get eigenvectors of original matrix 865 866 for (int j = nn-1; j >= low; j--) { 867 for (int i = low; i <= high; i++) { 868 z = 0.0; 869 for (int k = low; k <= Math.min(j,high); k++) { 870 z = z + V[i][k] * H[k][j]; 871 } 872 V[i][j] = z; 873 } 874 } 875 } 876 877 878/* ------------------------ 879 Constructor 880 * ------------------------ */ 881 882 /** Check for symmetry, then construct the eigenvalue decomposition 883 Structure to access D and V. 884 @param Arg Square matrix 885 */ 886 887 public EigenvalueDecomposition (Matrix Arg) { 888 double[][] A = Arg.getArray(); 889 n = Arg.getColumnDimension(); 890 V = new double[n][n]; 891 d = new double[n]; 892 e = new double[n]; 893 894 issymmetric = true; 895 for (int j = 0; (j < n) && issymmetric; j++) { 896 for (int i = 0; (i < n) && issymmetric; i++) { 897 issymmetric = (A[i][j] == A[j][i]); 898 } 899 } 900 901 if (issymmetric) { 902 for (int i = 0; i < n; i++) { 903 for (int j = 0; j < n; j++) { 904 V[i][j] = A[i][j]; 905 } 906 } 907 908 // Tridiagonalize. 909 tred2(); 910 911 // Diagonalize. 912 tql2(); 913 914 } else { 915 H = new double[n][n]; 916 ort = new double[n]; 917 918 for (int j = 0; j < n; j++) { 919 for (int i = 0; i < n; i++) { 920 H[i][j] = A[i][j]; 921 } 922 } 923 924 // Reduce to Hessenberg form. 925 orthes(); 926 927 // Reduce Hessenberg to real Schur form. 928 hqr2(); 929 } 930 } 931 932/* ------------------------ 933 Public Methods 934 * ------------------------ */ 935 936 /** Return the eigenvector matrix 937 @return V 938 */ 939 940 public Matrix getV () { 941 return new Matrix(V,n,n); 942 } 943 944 /** Return the real parts of the eigenvalues 945 @return real(diag(D)) 946 */ 947 948 public double[] getRealEigenvalues () { 949 return d; 950 } 951 952 /** Return the imaginary parts of the eigenvalues 953 @return imag(diag(D)) 954 */ 955 956 public double[] getImagEigenvalues () { 957 return e; 958 } 959 960 /** Return the block diagonal eigenvalue matrix 961 @return D 962 */ 963 964 public Matrix getD () { 965 Matrix X = new Matrix(n,n); 966 double[][] D = X.getArray(); 967 for (int i = 0; i < n; i++) { 968 for (int j = 0; j < n; j++) { 969 D[i][j] = 0.0; 970 } 971 D[i][i] = d[i]; 972 if (e[i] > 0) { 973 D[i][i+1] = e[i]; 974 } else if (e[i] < 0) { 975 D[i][i-1] = e[i]; 976 } 977 } 978 return X; 979 } 980 981}