00014 {
00015     /* System generated locals */
00016     integer z_dim1, z_offset, i__1, i__2, i__3;
00017     doublereal d__1, d__2;
00018 
00019     /* Builtin functions */
00020     double d_sign(doublereal *, doublereal *);
00021 
00022     /* Local variables */
00023     static doublereal b, c__, f, g;
00024     static integer i__, j, k, l, m;
00025     static doublereal p, r__, s;
00026     static integer ii;
00027     extern doublereal pythag_(doublereal *, doublereal *);
00028     static integer mml;
00029     static doublereal tst1, tst2;
00030 
00031 
00032 
00033 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, */
00034 /*     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, */
00035 /*     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. */
00036 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). */
00037 
00038 /*     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
00039 /*     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. */
00040 /*     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO */
00041 /*     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS */
00042 /*     FULL MATRIX TO TRIDIAGONAL FORM. */
00043 
00044 /*     ON INPUT */
00045 
00046 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00047 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00048 /*          DIMENSION STATEMENT. */
00049 
00050 /*        N IS THE ORDER OF THE MATRIX. */
00051 
00052 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00053 
00054 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00055 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00056 
00057 /*        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */
00058 /*          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS */
00059 /*          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN */
00060 /*          THE IDENTITY MATRIX. */
00061 
00062 /*      ON OUTPUT */
00063 
00064 /*        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
00065 /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT */
00066 /*          UNORDERED FOR INDICES 1,2,...,IERR-1. */
00067 
00068 /*        E HAS BEEN DESTROYED. */
00069 
00070 /*        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC */
00071 /*          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE, */
00072 /*          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED */
00073 /*          EIGENVALUES. */
00074 
00075 /*        IERR IS SET TO */
00076 /*          ZERO       FOR NORMAL RETURN, */
00077 /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
00078 /*                     DETERMINED AFTER 30 ITERATIONS. */
00079 
00080 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00081 
00082 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00083 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00084 */
00085 
00086 /*     THIS VERSION DATED AUGUST 1983. */
00087 
00088 /*     ------------------------------------------------------------------ 
00089 */
00090 
00091     /* Parameter adjustments */
00092     z_dim1 = *nm;
00093     z_offset = z_dim1 + 1;
00094     z__ -= z_offset;
00095     --e;
00096     --d__;
00097 
00098     /* Function Body */
00099     *ierr = 0;
00100     if (*n == 1) {
00101         goto L1001;
00102     }
00103 
00104     i__1 = *n;
00105     for (i__ = 2; i__ <= i__1; ++i__) {
00106 /* L100: */
00107         e[i__ - 1] = e[i__];
00108     }
00109 
00110     e[*n] = 0.;
00111 
00112     i__1 = *n;
00113     for (l = 1; l <= i__1; ++l) {
00114         j = 0;
00115 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
00116 L105:
00117         i__2 = *n;
00118         for (m = l; m <= i__2; ++m) {
00119             if (m == *n) {
00120                 goto L120;
00121             }
00122             tst1 = (d__1 = d__[m], abs(d__1)) + (d__2 = d__[m + 1], abs(d__2))
00123                     ;
00124             tst2 = tst1 + (d__1 = e[m], abs(d__1));
00125             if (tst2 == tst1) {
00126                 goto L120;
00127             }
00128 /* L110: */
00129         }
00130 
00131 L120:
00132         p = d__[l];
00133         if (m == l) {
00134             goto L240;
00135         }
00136         if (j == 30) {
00137             goto L1000;
00138         }
00139         ++j;
00140 /*     .......... FORM SHIFT .......... */
00141         g = (d__[l + 1] - p) / (e[l] * 2.);
00142         r__ = pythag_(&g, &c_b9);
00143         g = d__[m] - p + e[l] / (g + d_sign(&r__, &g));
00144         s = 1.;
00145         c__ = 1.;
00146         p = 0.;
00147         mml = m - l;
00148 /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
00149         i__2 = mml;
00150         for (ii = 1; ii <= i__2; ++ii) {
00151             i__ = m - ii;
00152             f = s * e[i__];
00153             b = c__ * e[i__];
00154             r__ = pythag_(&f, &g);
00155             e[i__ + 1] = r__;
00156             if (r__ == 0.) {
00157                 goto L210;
00158             }
00159             s = f / r__;
00160             c__ = g / r__;
00161             g = d__[i__ + 1] - p;
00162             r__ = (d__[i__] - g) * s + c__ * 2. * b;
00163             p = s * r__;
00164             d__[i__ + 1] = g + p;
00165             g = c__ * r__ - b;
00166 /*     .......... FORM VECTOR .......... */
00167             i__3 = *n;
00168             for (k = 1; k <= i__3; ++k) {
00169                 f = z__[k + (i__ + 1) * z_dim1];
00170                 z__[k + (i__ + 1) * z_dim1] = s * z__[k + i__ * z_dim1] + c__ 
00171                         * f;
00172                 z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * f;
00173 /* L180: */
00174             }
00175 
00176 /* L200: */
00177         }
00178 
00179         d__[l] -= p;
00180         e[l] = g;
00181         e[m] = 0.;
00182         goto L105;
00183 /*     .......... RECOVER FROM UNDERFLOW .......... */
00184 L210:
00185         d__[i__ + 1] -= p;
00186         e[m] = 0.;
00187         goto L105;
00188 L240:
00189         ;
00190     }
00191 /*     .......... ORDER EIGENVALUES AND EIGENVECTORS .......... */
00192     i__1 = *n;
00193     for (ii = 2; ii <= i__1; ++ii) {
00194         i__ = ii - 1;
00195         k = i__;
00196         p = d__[i__];
00197 
00198         i__2 = *n;
00199         for (j = ii; j <= i__2; ++j) {
00200             if (d__[j] >= p) {
00201                 goto L260;
00202             }
00203             k = j;
00204             p = d__[j];
00205 L260:
00206             ;
00207         }
00208 
00209         if (k == i__) {
00210             goto L300;
00211         }
00212         d__[k] = d__[i__];
00213         d__[i__] = p;
00214 
00215         i__2 = *n;
00216         for (j = 1; j <= i__2; ++j) {
00217             p = z__[j + i__ * z_dim1];
00218             z__[j + i__ * z_dim1] = z__[j + k * z_dim1];
00219             z__[j + k * z_dim1] = p;
00220 /* L280: */
00221         }
00222 
00223 L300:
00224         ;
00225     }
00226 
00227     goto L1001;
00228 /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
00229 /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
00230 L1000:
00231     *ierr = l;
00232 L1001:
00233     return 0;
00234 } /* imtql2_ */