src/javajs/util/Eigen.java

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