1 /* $RCSfile$ 2 * $Author: egonw $ 3 * $Date: 2005-11-10 09:52:44f -0600 (Thu, 10 Nov 2005) $ 4 * $Revision: 4255 $ 5 * 6 * Copyright (C) 2003-2005 Miguel, Jmol Development, www.jmol.org 7 * 8 * Contact: jmol-developers@lists.sf.net 9 * 10 * This library is free software; you can redistribute it and/or 11 * modify it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 23 */ 24 25 package javajs.util; 26 27 import javajs.api.EigenInterface; 28 29 30 /** 31 * Eigenvalues and eigenvectors of a real matrix. 32 * See javajs.api.EigenInterface() as well. 33 * 34 * adapted by Bob Hanson from http://math.nist.gov/javanumerics/jama/ (public 35 * domain); adding quaternion superimposition capability; removing 36 * nonsymmetric reduction to Hessenberg form, which we do not need in Jmol. 37 * 38 * Output is as a set of double[n] columns, but for the EigenInterface 39 * we return them as V3[3] and float[3] (or double[3]) values. 40 * 41 * Eigenvalues and eigenvectors are sorted from smallest to largest eigenvalue. 42 * 43 * <P> 44 * If A is symmetric, then A = V*D*V' where the eigenvalue matrix D is diagonal 45 * and the eigenvector matrix V is orthogonal. I.e. A = 46 * V.times(D.times(V.transpose())) and V.times(V.transpose()) equals the 47 * identity matrix. 48 * <P> 49 * If A is not symmetric, then the eigenvalue matrix D is block diagonal with 50 * the real eigenvalues in 1-by-1 blocks and any complex eigenvalues, lambda + 51 * i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda]. The columns of V represent 52 * the eigenvectors in the sense that A*V = V*D, i.e. A.times(V) equals 53 * V.times(D). The matrix V may be badly conditioned, or even singular, so the 54 * validity of the equation A = V*D*inverse(V) depends upon V.cond(). 55 **/ 56 57 public class Eigen implements EigenInterface { 58 59 /* ------------------------ 60 Public Methods 61 * ------------------------ */ 62 Eigen()63 public Eigen() {} 64 set(int n)65 public Eigen set(int n) { 66 this.n = n; 67 V = new double[n][n]; 68 d = new double[n]; 69 e = new double[n]; 70 return this; 71 } 72 73 @Override setM(double[][] m)74 public Eigen setM(double[][] m) { 75 set(m.length); 76 calc(m); 77 return this; 78 } 79 80 /** 81 * return values sorted from smallest to largest value. 82 */ 83 @Override getEigenvalues()84 public double[] getEigenvalues() { 85 return d; 86 } 87 88 /** 89 * Specifically for 3x3 systems, returns eigenVectors as V3[3] 90 * and values as float[3]; sorted from smallest to largest value. 91 * 92 * @param eigenVectors returned vectors 93 * @param eigenValues returned values 94 * 95 */ 96 @Override fillFloatArrays(V3[] eigenVectors, float[] eigenValues)97 public void fillFloatArrays(V3[] eigenVectors, float[] eigenValues) { 98 for (int i = 0; i < 3; i++) { 99 if (eigenVectors != null) { 100 if (eigenVectors[i] == null) 101 eigenVectors[i] = new V3(); 102 eigenVectors[i].set((float) V[0][i], (float) V[1][i], (float) V[2][i]); 103 } 104 if (eigenValues != null) 105 eigenValues[i] = (float) d[i]; 106 } 107 } 108 109 /** 110 * Transpose V and turn into floats; sorted from smallest to largest value. 111 * 112 * @return ROWS of eigenvectors f[0], f[1], f[2], etc. 113 */ 114 @Override getEigenvectorsFloatTransposed()115 public float[][] getEigenvectorsFloatTransposed() { 116 float[][] f = new float[n][n]; 117 for (int i = n; --i >= 0;) 118 for (int j = n; --j >= 0;) 119 f[j][i] = (float) V[i][j]; 120 return f; 121 } 122 123 124 /** 125 * Check for symmetry, then construct the eigenvalue decomposition 126 * 127 * @param A 128 * Square matrix 129 */ 130 calc(double[][] A)131 public void calc(double[][] A) { 132 133 /* Jmol only has need of symmetric solutions 134 * 135 issymmetric = true; 136 137 for (int j = 0; (j < n) & issymmetric; j++) { 138 for (int i = 0; (i < n) & issymmetric; i++) { 139 issymmetric = (A[i][j] == A[j][i]); 140 } 141 } 142 143 if (issymmetric) { 144 */ 145 for (int i = 0; i < n; i++) { 146 for (int j = 0; j < n; j++) { 147 V[i][j] = A[i][j]; 148 } 149 } 150 151 // Tridiagonalize. 152 tred2(); 153 154 // Diagonalize. 155 tql2(); 156 /* 157 } else { 158 H = new double[n][n]; 159 ort = new double[n]; 160 161 for (int j = 0; j < n; j++) { 162 for (int i = 0; i < n; i++) { 163 H[i][j] = A[i][j]; 164 } 165 } 166 167 // Reduce to Hessenberg form. 168 orthes(); 169 170 // Reduce Hessenberg to real Schur form. 171 hqr2(); 172 } 173 */ 174 175 } 176 177 /** 178 * Return the real parts of the eigenvalues 179 * 180 * @return real(diag(D)) 181 */ 182 getRealEigenvalues()183 public double[] getRealEigenvalues() { 184 return d; 185 } 186 187 /** 188 * Return the imaginary parts of the eigenvalues 189 * 190 * @return imag(diag(D)) 191 */ 192 getImagEigenvalues()193 public double[] getImagEigenvalues() { 194 return e; 195 } 196 197 /* ------------------------ 198 Class variables 199 * ------------------------ */ 200 201 /** 202 * Row and column dimension (square matrix). 203 * 204 * @serial matrix dimension. 205 */ 206 private int n = 3; 207 208 /** 209 * Symmetry flag. 210 * 211 * @serial internal symmetry flag. 212 */ 213 //private boolean issymmetric = true; 214 215 /** 216 * Arrays for internal storage of eigenvalues. 217 * 218 * @serial internal storage of eigenvalues. 219 */ 220 private double[] d, e; 221 222 /** 223 * Array for internal storage of eigenvectors. 224 * 225 * @serial internal storage of eigenvectors. 226 */ 227 private double[][] V; 228 229 /** 230 * Array for internal storage of nonsymmetric Hessenberg form. 231 * 232 * @serial internal storage of nonsymmetric Hessenberg form. 233 */ 234 //private double[][] H; 235 236 /** 237 * Working storage for nonsymmetric algorithm. 238 * 239 * @serial working storage for nonsymmetric algorithm. 240 */ 241 //private double[] ort; 242 243 /* ------------------------ 244 Private Methods 245 * ------------------------ */ 246 247 // Symmetric Householder reduction to tridiagonal form. 248 tred2()249 private void tred2() { 250 251 // This is derived from the Algol procedures tred2 by 252 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 253 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 254 // Fortran subroutine in EISPACK. 255 256 for (int j = 0; j < n; j++) { 257 d[j] = V[n - 1][j]; 258 } 259 260 // Householder reduction to tridiagonal form. 261 262 for (int i = n - 1; i > 0; i--) { 263 264 // Scale to avoid under/overflow. 265 266 double scale = 0.0; 267 double h = 0.0; 268 for (int k = 0; k < i; k++) { 269 scale = scale + Math.abs(d[k]); 270 } 271 if (scale == 0.0) { 272 e[i] = d[i - 1]; 273 for (int j = 0; j < i; j++) { 274 d[j] = V[i - 1][j]; 275 V[i][j] = 0.0; 276 V[j][i] = 0.0; 277 } 278 } else { 279 280 // Generate Householder vector. 281 282 for (int k = 0; k < i; k++) { 283 d[k] /= scale; 284 h += d[k] * d[k]; 285 } 286 double f = d[i - 1]; 287 double g = Math.sqrt(h); 288 if (f > 0) { 289 g = -g; 290 } 291 e[i] = scale * g; 292 h = h - f * g; 293 d[i - 1] = f - g; 294 for (int j = 0; j < i; j++) { 295 e[j] = 0.0; 296 } 297 298 // Apply similarity transformation to remaining columns. 299 300 for (int j = 0; j < i; j++) { 301 f = d[j]; 302 V[j][i] = f; 303 g = e[j] + V[j][j] * f; 304 for (int k = j + 1; k <= i - 1; k++) { 305 g += V[k][j] * d[k]; 306 e[k] += V[k][j] * f; 307 } 308 e[j] = g; 309 } 310 f = 0.0; 311 for (int j = 0; j < i; j++) { 312 e[j] /= h; 313 f += e[j] * d[j]; 314 } 315 double hh = f / (h + h); 316 for (int j = 0; j < i; j++) { 317 e[j] -= hh * d[j]; 318 } 319 for (int j = 0; j < i; j++) { 320 f = d[j]; 321 g = e[j]; 322 for (int k = j; k <= i - 1; k++) { 323 V[k][j] -= (f * e[k] + g * d[k]); 324 } 325 d[j] = V[i - 1][j]; 326 V[i][j] = 0.0; 327 } 328 } 329 d[i] = h; 330 } 331 332 // Accumulate transformations. 333 334 for (int i = 0; i < n - 1; i++) { 335 V[n - 1][i] = V[i][i]; 336 V[i][i] = 1.0; 337 double h = d[i + 1]; 338 if (h != 0.0) { 339 for (int k = 0; k <= i; k++) { 340 d[k] = V[k][i + 1] / h; 341 } 342 for (int j = 0; j <= i; j++) { 343 double g = 0.0; 344 for (int k = 0; k <= i; k++) { 345 g += V[k][i + 1] * V[k][j]; 346 } 347 for (int k = 0; k <= i; k++) { 348 V[k][j] -= g * d[k]; 349 } 350 } 351 } 352 for (int k = 0; k <= i; k++) { 353 V[k][i + 1] = 0.0; 354 } 355 } 356 for (int j = 0; j < n; j++) { 357 d[j] = V[n - 1][j]; 358 V[n - 1][j] = 0.0; 359 } 360 V[n - 1][n - 1] = 1.0; 361 e[0] = 0.0; 362 } 363 364 // Symmetric tridiagonal QL algorithm. 365 tql2()366 private void tql2() { 367 368 // This is derived from the Algol procedures tql2, by 369 // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for 370 // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 371 // Fortran subroutine in EISPACK. 372 373 for (int i = 1; i < n; i++) { 374 e[i - 1] = e[i]; 375 } 376 e[n - 1] = 0.0; 377 378 double f = 0.0; 379 double tst1 = 0.0; 380 double eps = Math.pow(2.0, -52.0); 381 for (int l = 0; l < n; l++) { 382 383 // Find small subdiagonal element 384 385 tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l])); 386 int m = l; 387 while (m < n) { 388 if (Math.abs(e[m]) <= eps * tst1) { 389 break; 390 } 391 m++; 392 } 393 394 // If m == l, d[l] is an eigenvalue, 395 // otherwise, iterate. 396 397 if (m > l) { 398 int iter = 0; 399 do { 400 iter = iter + 1; // (Could check iteration count here.) 401 402 // Compute implicit shift 403 404 double g = d[l]; 405 double p = (d[l + 1] - g) / (2.0 * e[l]); 406 double r = hypot(p, 1.0); 407 if (p < 0) { 408 r = -r; 409 } 410 d[l] = e[l] / (p + r); 411 d[l + 1] = e[l] * (p + r); 412 double dl1 = d[l + 1]; 413 double h = g - d[l]; 414 for (int i = l + 2; i < n; i++) { 415 d[i] -= h; 416 } 417 f = f + h; 418 419 // Implicit QL transformation. 420 421 p = d[m]; 422 double c = 1.0; 423 double c2 = c; 424 double c3 = c; 425 double el1 = e[l + 1]; 426 double s = 0.0; 427 double s2 = 0.0; 428 for (int i = m - 1; i >= l; i--) { 429 c3 = c2; 430 c2 = c; 431 s2 = s; 432 g = c * e[i]; 433 h = c * p; 434 r = hypot(p, e[i]); 435 e[i + 1] = s * r; 436 s = e[i] / r; 437 c = p / r; 438 p = c * d[i] - s * g; 439 d[i + 1] = h + s * (c * g + s * d[i]); 440 441 // Accumulate transformation. 442 443 for (int k = 0; k < n; k++) { 444 h = V[k][i + 1]; 445 V[k][i + 1] = s * V[k][i] + c * h; 446 V[k][i] = c * V[k][i] - s * h; 447 } 448 } 449 p = -s * s2 * c3 * el1 * e[l] / dl1; 450 e[l] = s * p; 451 d[l] = c * p; 452 453 // Check for convergence. 454 455 } while (Math.abs(e[l]) > eps * tst1); 456 } 457 d[l] = d[l] + f; 458 e[l] = 0.0; 459 } 460 461 // Sort eigenvalues and corresponding vectors. 462 463 for (int i = 0; i < n - 1; i++) { 464 int k = i; 465 double p = d[i]; 466 for (int j = i + 1; j < n; j++) { 467 if (d[j] < p) { 468 k = j; 469 p = d[j]; 470 } 471 } 472 if (k != i) { 473 d[k] = d[i]; 474 d[i] = p; 475 for (int j = 0; j < n; j++) { 476 p = V[j][i]; 477 V[j][i] = V[j][k]; 478 V[j][k] = p; 479 } 480 } 481 } 482 } 483 hypot(double a, double b)484 private static double hypot(double a, double b) { 485 486 // sqrt(a^2 + b^2) without under/overflow. 487 488 double r; 489 if (Math.abs(a) > Math.abs(b)) { 490 r = b / a; 491 r = Math.abs(a) * Math.sqrt(1 + r * r); 492 } else if (b != 0) { 493 r = a / b; 494 r = Math.abs(b) * Math.sqrt(1 + r * r); 495 } else { 496 r = 0.0; 497 } 498 return r; 499 } 500 501 // Nonsymmetric reduction to Hessenberg form. 502 503 /* 504 private void orthes() { 505 506 // This is derived from the Algol procedures orthes and ortran, 507 // by Martin and Wilkinson, Handbook for Auto. Comp., 508 // Vol.ii-Linear Algebra, and the corresponding 509 // Fortran subroutines in EISPACK. 510 511 int low = 0; 512 int high = n - 1; 513 514 for (int m = low + 1; m <= high - 1; m++) { 515 516 // Scale column. 517 518 double scale = 0.0; 519 for (int i = m; i <= high; i++) { 520 scale = scale + Math.abs(H[i][m - 1]); 521 } 522 if (scale != 0.0) { 523 524 // Compute Householder transformation. 525 526 double h = 0.0; 527 for (int i = high; i >= m; i--) { 528 ort[i] = H[i][m - 1] / scale; 529 h += ort[i] * ort[i]; 530 } 531 double g = Math.sqrt(h); 532 if (ort[m] > 0) { 533 g = -g; 534 } 535 h = h - ort[m] * g; 536 ort[m] = ort[m] - g; 537 538 // Apply Householder similarity transformation 539 // H = (I-u*u'/h)*H*(I-u*u')/h) 540 541 for (int j = m; j < n; j++) { 542 double f = 0.0; 543 for (int i = high; i >= m; i--) { 544 f += ort[i] * H[i][j]; 545 } 546 f = f / h; 547 for (int i = m; i <= high; i++) { 548 H[i][j] -= f * ort[i]; 549 } 550 } 551 552 for (int i = 0; i <= high; i++) { 553 double f = 0.0; 554 for (int j = high; j >= m; j--) { 555 f += ort[j] * H[i][j]; 556 } 557 f = f / h; 558 for (int j = m; j <= high; j++) { 559 H[i][j] -= f * ort[j]; 560 } 561 } 562 ort[m] = scale * ort[m]; 563 H[m][m - 1] = scale * g; 564 } 565 } 566 567 // Accumulate transformations (Algol's ortran). 568 569 for (int i = 0; i < n; i++) { 570 for (int j = 0; j < n; j++) { 571 V[i][j] = (i == j ? 1.0 : 0.0); 572 } 573 } 574 575 for (int m = high - 1; m >= low + 1; m--) { 576 if (H[m][m - 1] != 0.0) { 577 for (int i = m + 1; i <= high; i++) { 578 ort[i] = H[i][m - 1]; 579 } 580 for (int j = m; j <= high; j++) { 581 double g = 0.0; 582 for (int i = m; i <= high; i++) { 583 g += ort[i] * V[i][j]; 584 } 585 // Double division avoids possible underflow 586 g = (g / ort[m]) / H[m][m - 1]; 587 for (int i = m; i <= high; i++) { 588 V[i][j] += g * ort[i]; 589 } 590 } 591 } 592 } 593 } 594 595 // Complex scalar division. 596 597 private transient double cdivr, cdivi; 598 599 private void cdiv(double xr, double xi, double yr, double yi) { 600 double r, d; 601 if (Math.abs(yr) > Math.abs(yi)) { 602 r = yi / yr; 603 d = yr + r * yi; 604 cdivr = (xr + r * xi) / d; 605 cdivi = (xi - r * xr) / d; 606 } else { 607 r = yr / yi; 608 d = yi + r * yr; 609 cdivr = (r * xr + xi) / d; 610 cdivi = (r * xi - xr) / d; 611 } 612 } 613 614 // Nonsymmetric reduction from Hessenberg to real Schur form. 615 616 private void hqr2() { 617 618 // This is derived from the Algol procedure hqr2, 619 // by Martin and Wilkinson, Handbook for Auto. Comp., 620 // Vol.ii-Linear Algebra, and the corresponding 621 // Fortran subroutine in EISPACK. 622 623 // Initialize 624 625 int nn = this.n; 626 int n = nn - 1; 627 int low = 0; 628 int high = nn - 1; 629 double eps = Math.pow(2.0, -52.0); 630 double exshift = 0.0; 631 double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y; 632 633 // Store roots isolated by balanc and compute matrix norm 634 635 double norm = 0.0; 636 for (int i = 0; i < nn; i++) { 637 if (i < low || i > high) { 638 d[i] = H[i][i]; 639 e[i] = 0.0; 640 } 641 for (int j = Math.max(i - 1, 0); j < nn; j++) { 642 norm = norm + Math.abs(H[i][j]); 643 } 644 } 645 646 // Outer loop over eigenvalue index 647 648 int iter = 0; 649 while (n >= low) { 650 651 // Look for single small sub-diagonal element 652 653 int l = n; 654 while (l > low) { 655 s = Math.abs(H[l - 1][l - 1]) + Math.abs(H[l][l]); 656 if (s == 0.0) { 657 s = norm; 658 } 659 if (Math.abs(H[l][l - 1]) < eps * s) { 660 break; 661 } 662 l--; 663 } 664 665 // Check for convergence 666 // One root found 667 668 if (l == n) { 669 H[n][n] = H[n][n] + exshift; 670 d[n] = H[n][n]; 671 e[n] = 0.0; 672 n--; 673 iter = 0; 674 675 // Two roots found 676 677 } else if (l == n - 1) { 678 w = H[n][n - 1] * H[n - 1][n]; 679 p = (H[n - 1][n - 1] - H[n][n]) / 2.0; 680 q = p * p + w; 681 z = Math.sqrt(Math.abs(q)); 682 H[n][n] = H[n][n] + exshift; 683 H[n - 1][n - 1] = H[n - 1][n - 1] + exshift; 684 x = H[n][n]; 685 686 // Real pair 687 688 if (q >= 0) { 689 if (p >= 0) { 690 z = p + z; 691 } else { 692 z = p - z; 693 } 694 d[n - 1] = x + z; 695 d[n] = d[n - 1]; 696 if (z != 0.0) { 697 d[n] = x - w / z; 698 } 699 e[n - 1] = 0.0; 700 e[n] = 0.0; 701 x = H[n][n - 1]; 702 s = Math.abs(x) + Math.abs(z); 703 p = x / s; 704 q = z / s; 705 r = Math.sqrt(p * p + q * q); 706 p = p / r; 707 q = q / r; 708 709 // Row modification 710 711 for (int j = n - 1; j < nn; j++) { 712 z = H[n - 1][j]; 713 H[n - 1][j] = q * z + p * H[n][j]; 714 H[n][j] = q * H[n][j] - p * z; 715 } 716 717 // Column modification 718 719 for (int i = 0; i <= n; i++) { 720 z = H[i][n - 1]; 721 H[i][n - 1] = q * z + p * H[i][n]; 722 H[i][n] = q * H[i][n] - p * z; 723 } 724 725 // Accumulate transformations 726 727 for (int i = low; i <= high; i++) { 728 z = V[i][n - 1]; 729 V[i][n - 1] = q * z + p * V[i][n]; 730 V[i][n] = q * V[i][n] - p * z; 731 } 732 733 // Complex pair 734 735 } else { 736 d[n - 1] = x + p; 737 d[n] = x + p; 738 e[n - 1] = z; 739 e[n] = -z; 740 } 741 n = n - 2; 742 iter = 0; 743 744 // No convergence yet 745 746 } else { 747 748 // Form shift 749 750 x = H[n][n]; 751 y = 0.0; 752 w = 0.0; 753 if (l < n) { 754 y = H[n - 1][n - 1]; 755 w = H[n][n - 1] * H[n - 1][n]; 756 } 757 758 // Wilkinson's original ad hoc shift 759 760 if (iter == 10) { 761 exshift += x; 762 for (int i = low; i <= n; i++) { 763 H[i][i] -= x; 764 } 765 s = Math.abs(H[n][n - 1]) + Math.abs(H[n - 1][n - 2]); 766 x = y = 0.75 * s; 767 w = -0.4375 * s * s; 768 } 769 770 // MATLAB's new ad hoc shift 771 772 if (iter == 30) { 773 s = (y - x) / 2.0; 774 s = s * s + w; 775 if (s > 0) { 776 s = Math.sqrt(s); 777 if (y < x) { 778 s = -s; 779 } 780 s = x - w / ((y - x) / 2.0 + s); 781 for (int i = low; i <= n; i++) { 782 H[i][i] -= s; 783 } 784 exshift += s; 785 x = y = w = 0.964; 786 } 787 } 788 789 iter = iter + 1; // (Could check iteration count here.) 790 791 // Look for two consecutive small sub-diagonal elements 792 793 int m = n - 2; 794 while (m >= l) { 795 z = H[m][m]; 796 r = x - z; 797 s = y - z; 798 p = (r * s - w) / H[m + 1][m] + H[m][m + 1]; 799 q = H[m + 1][m + 1] - z - r - s; 800 r = H[m + 2][m + 1]; 801 s = Math.abs(p) + Math.abs(q) + Math.abs(r); 802 p = p / s; 803 q = q / s; 804 r = r / s; 805 if (m == l) { 806 break; 807 } 808 if (Math.abs(H[m][m - 1]) * (Math.abs(q) + Math.abs(r)) < eps 809 * (Math.abs(p) * (Math.abs(H[m - 1][m - 1]) + Math.abs(z) + Math 810 .abs(H[m + 1][m + 1])))) { 811 break; 812 } 813 m--; 814 } 815 816 for (int i = m + 2; i <= n; i++) { 817 H[i][i - 2] = 0.0; 818 if (i > m + 2) { 819 H[i][i - 3] = 0.0; 820 } 821 } 822 823 // Double QR step involving rows l:n and columns m:n 824 825 for (int k = m; k <= n - 1; k++) { 826 boolean notlast = (k != n - 1); 827 if (k != m) { 828 p = H[k][k - 1]; 829 q = H[k + 1][k - 1]; 830 r = (notlast ? H[k + 2][k - 1] : 0.0); 831 x = Math.abs(p) + Math.abs(q) + Math.abs(r); 832 if (x != 0.0) { 833 p = p / x; 834 q = q / x; 835 r = r / x; 836 } 837 } 838 if (x == 0.0) { 839 break; 840 } 841 s = Math.sqrt(p * p + q * q + r * r); 842 if (p < 0) { 843 s = -s; 844 } 845 if (s != 0) { 846 if (k != m) { 847 H[k][k - 1] = -s * x; 848 } else if (l != m) { 849 H[k][k - 1] = -H[k][k - 1]; 850 } 851 p = p + s; 852 x = p / s; 853 y = q / s; 854 z = r / s; 855 q = q / p; 856 r = r / p; 857 858 // Row modification 859 860 for (int j = k; j < nn; j++) { 861 p = H[k][j] + q * H[k + 1][j]; 862 if (notlast) { 863 p = p + r * H[k + 2][j]; 864 H[k + 2][j] = H[k + 2][j] - p * z; 865 } 866 H[k][j] = H[k][j] - p * x; 867 H[k + 1][j] = H[k + 1][j] - p * y; 868 } 869 870 // Column modification 871 872 for (int i = 0; i <= Math.min(n, k + 3); i++) { 873 p = x * H[i][k] + y * H[i][k + 1]; 874 if (notlast) { 875 p = p + z * H[i][k + 2]; 876 H[i][k + 2] = H[i][k + 2] - p * r; 877 } 878 H[i][k] = H[i][k] - p; 879 H[i][k + 1] = H[i][k + 1] - p * q; 880 } 881 882 // Accumulate transformations 883 884 for (int i = low; i <= high; i++) { 885 p = x * V[i][k] + y * V[i][k + 1]; 886 if (notlast) { 887 p = p + z * V[i][k + 2]; 888 V[i][k + 2] = V[i][k + 2] - p * r; 889 } 890 V[i][k] = V[i][k] - p; 891 V[i][k + 1] = V[i][k + 1] - p * q; 892 } 893 } // (s != 0) 894 } // k loop 895 } // check convergence 896 } // while (n >= low) 897 898 // Backsubstitute to find vectors of upper triangular form 899 900 if (norm == 0.0) { 901 return; 902 } 903 904 for (n = nn - 1; n >= 0; n--) { 905 p = d[n]; 906 q = e[n]; 907 908 // Real vector 909 910 if (q == 0) { 911 int l = n; 912 H[n][n] = 1.0; 913 for (int i = n - 1; i >= 0; i--) { 914 w = H[i][i] - p; 915 r = 0.0; 916 for (int j = l; j <= n; j++) { 917 r = r + H[i][j] * H[j][n]; 918 } 919 if (e[i] < 0.0) { 920 z = w; 921 s = r; 922 } else { 923 l = i; 924 if (e[i] == 0.0) { 925 if (w != 0.0) { 926 H[i][n] = -r / w; 927 } else { 928 H[i][n] = -r / (eps * norm); 929 } 930 931 // Solve real equations 932 933 } else { 934 x = H[i][i + 1]; 935 y = H[i + 1][i]; 936 q = (d[i] - p) * (d[i] - p) + e[i] * e[i]; 937 t = (x * s - z * r) / q; 938 H[i][n] = t; 939 if (Math.abs(x) > Math.abs(z)) { 940 H[i + 1][n] = (-r - w * t) / x; 941 } else { 942 H[i + 1][n] = (-s - y * t) / z; 943 } 944 } 945 946 // Overflow control 947 948 t = Math.abs(H[i][n]); 949 if ((eps * t) * t > 1) { 950 for (int j = i; j <= n; j++) { 951 H[j][n] = H[j][n] / t; 952 } 953 } 954 } 955 } 956 957 // Complex vector 958 959 } else if (q < 0) { 960 int l = n - 1; 961 962 // Last vector component imaginary so matrix is triangular 963 964 if (Math.abs(H[n][n - 1]) > Math.abs(H[n - 1][n])) { 965 H[n - 1][n - 1] = q / H[n][n - 1]; 966 H[n - 1][n] = -(H[n][n] - p) / H[n][n - 1]; 967 } else { 968 cdiv(0.0, -H[n - 1][n], H[n - 1][n - 1] - p, q); 969 H[n - 1][n - 1] = cdivr; 970 H[n - 1][n] = cdivi; 971 } 972 H[n][n - 1] = 0.0; 973 H[n][n] = 1.0; 974 for (int i = n - 2; i >= 0; i--) { 975 double ra, sa, vr, vi; 976 ra = 0.0; 977 sa = 0.0; 978 for (int j = l; j <= n; j++) { 979 ra = ra + H[i][j] * H[j][n - 1]; 980 sa = sa + H[i][j] * H[j][n]; 981 } 982 w = H[i][i] - p; 983 984 if (e[i] < 0.0) { 985 z = w; 986 r = ra; 987 s = sa; 988 } else { 989 l = i; 990 if (e[i] == 0) { 991 cdiv(-ra, -sa, w, q); 992 H[i][n - 1] = cdivr; 993 H[i][n] = cdivi; 994 } else { 995 996 // Solve complex equations 997 998 x = H[i][i + 1]; 999 y = H[i + 1][i]; 1000 vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q; 1001 vi = (d[i] - p) * 2.0 * q; 1002 if (vr == 0.0 & vi == 0.0) { 1003 vr = eps 1004 * norm 1005 * (Math.abs(w) + Math.abs(q) + Math.abs(x) + Math.abs(y) + Math 1006 .abs(z)); 1007 } 1008 cdiv(x * r - z * ra + q * sa, x * s - z * sa - q * ra, vr, vi); 1009 H[i][n - 1] = cdivr; 1010 H[i][n] = cdivi; 1011 if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) { 1012 H[i + 1][n - 1] = (-ra - w * H[i][n - 1] + q * H[i][n]) / x; 1013 H[i + 1][n] = (-sa - w * H[i][n] - q * H[i][n - 1]) / x; 1014 } else { 1015 cdiv(-r - y * H[i][n - 1], -s - y * H[i][n], z, q); 1016 H[i + 1][n - 1] = cdivr; 1017 H[i + 1][n] = cdivi; 1018 } 1019 } 1020 1021 // Overflow control 1022 1023 t = Math.max(Math.abs(H[i][n - 1]), Math.abs(H[i][n])); 1024 if ((eps * t) * t > 1) { 1025 for (int j = i; j <= n; j++) { 1026 H[j][n - 1] = H[j][n - 1] / t; 1027 H[j][n] = H[j][n] / t; 1028 } 1029 } 1030 } 1031 } 1032 } 1033 } 1034 1035 // Vectors of isolated roots 1036 1037 for (int i = 0; i < nn; i++) { 1038 if (i < low || i > high) { 1039 for (int j = i; j < nn; j++) { 1040 V[i][j] = H[i][j]; 1041 } 1042 } 1043 } 1044 1045 // Back transformation to get eigenvectors of original matrix 1046 1047 for (int j = nn - 1; j >= low; j--) { 1048 for (int i = low; i <= high; i++) { 1049 z = 0.0; 1050 for (int k = low; k <= Math.min(j, high); k++) { 1051 z = z + V[i][k] * H[k][j]; 1052 } 1053 V[i][j] = z; 1054 } 1055 } 1056 } 1057 */ 1058 1059 1060 } 1061