src/javajs/util/Eigen.java

   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
  63   public Eigen() {}
  64
  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
  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
  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
  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
 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
 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
 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
 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
 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
 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
 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 }