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