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 b, c__, f, g, h__;
00024     static integer i__, j, l, m;
00025     static doublereal p, r__, s, t;
00026     static integer l1, ii;
00027     extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 
00028             *);
00029     static integer mml;
00030 
00031 
00032 
00033 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, */
00034 /*     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. */
00035 
00036 /*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC */
00037 /*     TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. */
00038 
00039 /*     ON INPUT */
00040 
00041 /*        N IS THE ORDER OF THE MATRIX. */
00042 
00043 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00044 
00045 /*        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE */
00046 /*          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY. 
00047 */
00048 
00049 /*      ON OUTPUT */
00050 
00051 /*        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
00052 /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
00053 /*          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
00054 /*          THE SMALLEST EIGENVALUES. */
00055 
00056 /*        E2 HAS BEEN DESTROYED. */
00057 
00058 /*        IERR IS SET TO */
00059 /*          ZERO       FOR NORMAL RETURN, */
00060 /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
00061 /*                     DETERMINED AFTER 30 ITERATIONS. */
00062 
00063 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00064 
00065 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00066 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00067 */
00068 
00069 /*     THIS VERSION DATED AUGUST 1983. */
00070 
00071 /*     ------------------------------------------------------------------ 
00072 */
00073 
00074     /* Parameter adjustments */
00075     --e2;
00076     --d__;
00077 
00078     /* Function Body */
00079     *ierr = 0;
00080     if (*n == 1) {
00081         goto L1001;
00082     }
00083 
00084     i__1 = *n;
00085     for (i__ = 2; i__ <= i__1; ++i__) {
00086 /* L100: */
00087         e2[i__ - 1] = e2[i__];
00088     }
00089 
00090     f = 0.;
00091     t = 0.;
00092     e2[*n] = 0.;
00093 
00094     i__1 = *n;
00095     for (l = 1; l <= i__1; ++l) {
00096         j = 0;
00097         h__ = (d__1 = d__[l], abs(d__1)) + sqrt(e2[l]);
00098         if (t > h__) {
00099             goto L105;
00100         }
00101         t = h__;
00102         b = epslon_(&t);
00103         c__ = b * b;
00104 /*     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ........
00105 .. */
00106 L105:
00107         i__2 = *n;
00108         for (m = l; m <= i__2; ++m) {
00109             if (e2[m] <= c__) {
00110                 goto L120;
00111             }
00112 /*     .......... E2(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         s = sqrt(e2[l]);
00129         g = d__[l];
00130         p = (d__[l1] - g) / (s * 2.);
00131         r__ = pythag_(&p, &c_b11);
00132         d__[l] = s / (p + d_sign(&r__, &p));
00133         h__ = g - d__[l];
00134 
00135         i__2 = *n;
00136         for (i__ = l1; i__ <= i__2; ++i__) {
00137 /* L140: */
00138             d__[i__] -= h__;
00139         }
00140 
00141         f += h__;
00142 /*     .......... RATIONAL QL TRANSFORMATION .......... */
00143         g = d__[m];
00144         if (g == 0.) {
00145             g = b;
00146         }
00147         h__ = g;
00148         s = 0.;
00149         mml = m - l;
00150 /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
00151         i__2 = mml;
00152         for (ii = 1; ii <= i__2; ++ii) {
00153             i__ = m - ii;
00154             p = g * h__;
00155             r__ = p + e2[i__];
00156             e2[i__ + 1] = s * r__;
00157             s = e2[i__] / r__;
00158             d__[i__ + 1] = h__ + s * (h__ + d__[i__]);
00159             g = d__[i__] - e2[i__] / g;
00160             if (g == 0.) {
00161                 g = b;
00162             }
00163             h__ = g * p / r__;
00164 /* L200: */
00165         }
00166 
00167         e2[l] = s * g;
00168         d__[l] = h__;
00169 /*     .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST ........
00170 .. */
00171         if (h__ == 0.) {
00172             goto L210;
00173         }
00174         if ((d__1 = e2[l], abs(d__1)) <= (d__2 = c__ / h__, abs(d__2))) {
00175             goto L210;
00176         }
00177         e2[l] = h__ * e2[l];
00178         if (e2[l] != 0.) {
00179             goto L130;
00180         }
00181 L210:
00182         p = d__[l] + f;
00183 /*     .......... ORDER EIGENVALUES .......... */
00184         if (l == 1) {
00185             goto L250;
00186         }
00187 /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
00188         i__2 = l;
00189         for (ii = 2; ii <= i__2; ++ii) {
00190             i__ = l + 2 - ii;
00191             if (p >= d__[i__ - 1]) {
00192                 goto L270;
00193             }
00194             d__[i__] = d__[i__ - 1];
00195 /* L230: */
00196         }
00197 
00198 L250:
00199         i__ = 1;
00200 L270:
00201         d__[i__] = p;
00202 /* L290: */
00203     }
00204 
00205     goto L1001;
00206 /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
00207 /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
00208 L1000:
00209     *ierr = l;
00210 L1001:
00211     return 0;
00212 } /* tqlrat_ */