00014 {
00015     /* System generated locals */
00016     integer i__1, i__2;
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, 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 TQL1, */
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 OF A SYMMETRIC */
00042 /*     TRIDIAGONAL MATRIX BY THE QL METHOD. */
00043 
00044 /*     ON INPUT */
00045 
00046 /*        N IS THE ORDER OF THE MATRIX. */
00047 
00048 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00049 
00050 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00051 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00052 
00053 /*      ON OUTPUT */
00054 
00055 /*        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
00056 /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
00057 /*          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
00058 /*          THE SMALLEST EIGENVALUES. */
00059 
00060 /*        E HAS BEEN DESTROYED. */
00061 
00062 /*        IERR IS SET TO */
00063 /*          ZERO       FOR NORMAL RETURN, */
00064 /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
00065 /*                     DETERMINED AFTER 30 ITERATIONS. */
00066 
00067 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00068 
00069 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00070 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00071 */
00072 
00073 /*     THIS VERSION DATED AUGUST 1983. */
00074 
00075 /*     ------------------------------------------------------------------ 
00076 */
00077 
00078     /* Parameter adjustments */
00079     --e;
00080     --d__;
00081 
00082     /* Function Body */
00083     *ierr = 0;
00084     if (*n == 1) {
00085         goto L1001;
00086     }
00087 
00088     i__1 = *n;
00089     for (i__ = 2; i__ <= i__1; ++i__) {
00090 /* L100: */
00091         e[i__ - 1] = e[i__];
00092     }
00093 
00094     f = 0.;
00095     tst1 = 0.;
00096     e[*n] = 0.;
00097 
00098     i__1 = *n;
00099     for (l = 1; l <= i__1; ++l) {
00100         j = 0;
00101         h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
00102         if (tst1 < h__) {
00103             tst1 = h__;
00104         }
00105 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
00106         i__2 = *n;
00107         for (m = l; m <= i__2; ++m) {
00108             tst2 = tst1 + (d__1 = e[m], abs(d__1));
00109             if (tst2 == tst1) {
00110                 goto L120;
00111             }
00112 /*     .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */
00113 /*                THROUGH THE BOTTOM OF THE LOOP .......... */
00114 /* L110: */
00115         }
00116 
00117 L120:
00118         if (m == l) {
00119             goto L210;
00120         }
00121 L130:
00122         if (j == 30) {
00123             goto L1000;
00124         }
00125         ++j;
00126 /*     .......... FORM SHIFT .......... */
00127         l1 = l + 1;
00128         l2 = l1 + 1;
00129         g = d__[l];
00130         p = (d__[l1] - g) / (e[l] * 2.);
00131         r__ = pythag_(&p, &c_b10);
00132         d__[l] = e[l] / (p + d_sign(&r__, &p));
00133         d__[l1] = e[l] * (p + d_sign(&r__, &p));
00134         dl1 = d__[l1];
00135         h__ = g - d__[l];
00136         if (l2 > *n) {
00137             goto L145;
00138         }
00139 
00140         i__2 = *n;
00141         for (i__ = l2; i__ <= i__2; ++i__) {
00142 /* L140: */
00143             d__[i__] -= h__;
00144         }
00145 
00146 L145:
00147         f += h__;
00148 /*     .......... QL TRANSFORMATION .......... */
00149         p = d__[m];
00150         c__ = 1.;
00151         c2 = c__;
00152         el1 = e[l1];
00153         s = 0.;
00154         mml = m - l;
00155 /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
00156         i__2 = mml;
00157         for (ii = 1; ii <= i__2; ++ii) {
00158             c3 = c2;
00159             c2 = c__;
00160             s2 = s;
00161             i__ = m - ii;
00162             g = c__ * e[i__];
00163             h__ = c__ * p;
00164             r__ = pythag_(&p, &e[i__]);
00165             e[i__ + 1] = s * r__;
00166             s = e[i__] / r__;
00167             c__ = p / r__;
00168             p = c__ * d__[i__] - s * g;
00169             d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
00170 /* L200: */
00171         }
00172 
00173         p = -s * s2 * c3 * el1 * e[l] / dl1;
00174         e[l] = s * p;
00175         d__[l] = c__ * p;
00176         tst2 = tst1 + (d__1 = e[l], abs(d__1));
00177         if (tst2 > tst1) {
00178             goto L130;
00179         }
00180 L210:
00181         p = d__[l] + f;
00182 /*     .......... ORDER EIGENVALUES .......... */
00183         if (l == 1) {
00184             goto L250;
00185         }
00186 /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
00187         i__2 = l;
00188         for (ii = 2; ii <= i__2; ++ii) {
00189             i__ = l + 2 - ii;
00190             if (p >= d__[i__ - 1]) {
00191                 goto L270;
00192             }
00193             d__[i__] = d__[i__ - 1];
00194 /* L230: */
00195         }
00196 
00197 L250:
00198         i__ = 1;
00199 L270:
00200         d__[i__] = p;
00201 /* L290: */
00202     }
00203 
00204     goto L1001;
00205 /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
00206 /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
00207 L1000:
00208     *ierr = l;
00209 L1001:
00210     return 0;
00211 } /* tql1_ */