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