00001 /* imtql1.f -- translated by f2c (version 19961017).
00002    You must link the resulting object file with the libraries:
00003         -lf2c -lm   (in that order)
00004 */
00005 
00006 #include "f2c.h"
00007 
00008 /* Table of constant values */
00009 
00010 static doublereal c_b10 = 1.;
00011 
00012 /* Subroutine */ int imtql1_(integer *n, doublereal *d__, doublereal *e, 
00013         integer *ierr)
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;
00024     static integer i__, j, 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 IMTQL1, */
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 OF A SYMMETRIC */
00039 /*     TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. */
00040 
00041 /*     ON INPUT */
00042 
00043 /*        N IS THE ORDER OF THE MATRIX. */
00044 
00045 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00046 
00047 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00048 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00049 
00050 /*      ON OUTPUT */
00051 
00052 /*        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
00053 /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
00054 /*          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
00055 /*          THE SMALLEST EIGENVALUES. */
00056 
00057 /*        E HAS BEEN DESTROYED. */
00058 
00059 /*        IERR IS SET TO */
00060 /*          ZERO       FOR NORMAL RETURN, */
00061 /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
00062 /*                     DETERMINED AFTER 30 ITERATIONS. */
00063 
00064 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00065 
00066 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00067 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00068 */
00069 
00070 /*     THIS VERSION DATED AUGUST 1983. */
00071 
00072 /*     ------------------------------------------------------------------ 
00073 */
00074 
00075     /* Parameter adjustments */
00076     --e;
00077     --d__;
00078 
00079     /* Function Body */
00080     *ierr = 0;
00081     if (*n == 1) {
00082         goto L1001;
00083     }
00084 
00085     i__1 = *n;
00086     for (i__ = 2; i__ <= i__1; ++i__) {
00087 /* L100: */
00088         e[i__ - 1] = e[i__];
00089     }
00090 
00091     e[*n] = 0.;
00092 
00093     i__1 = *n;
00094     for (l = 1; l <= i__1; ++l) {
00095         j = 0;
00096 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
00097 L105:
00098         i__2 = *n;
00099         for (m = l; m <= i__2; ++m) {
00100             if (m == *n) {
00101                 goto L120;
00102             }
00103             tst1 = (d__1 = d__[m], abs(d__1)) + (d__2 = d__[m + 1], abs(d__2))
00104                     ;
00105             tst2 = tst1 + (d__1 = e[m], abs(d__1));
00106             if (tst2 == tst1) {
00107                 goto L120;
00108             }
00109 /* L110: */
00110         }
00111 
00112 L120:
00113         p = d__[l];
00114         if (m == l) {
00115             goto L215;
00116         }
00117         if (j == 30) {
00118             goto L1000;
00119         }
00120         ++j;
00121 /*     .......... FORM SHIFT .......... */
00122         g = (d__[l + 1] - p) / (e[l] * 2.);
00123         r__ = pythag_(&g, &c_b10);
00124         g = d__[m] - p + e[l] / (g + d_sign(&r__, &g));
00125         s = 1.;
00126         c__ = 1.;
00127         p = 0.;
00128         mml = m - l;
00129 /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
00130         i__2 = mml;
00131         for (ii = 1; ii <= i__2; ++ii) {
00132             i__ = m - ii;
00133             f = s * e[i__];
00134             b = c__ * e[i__];
00135             r__ = pythag_(&f, &g);
00136             e[i__ + 1] = r__;
00137             if (r__ == 0.) {
00138                 goto L210;
00139             }
00140             s = f / r__;
00141             c__ = g / r__;
00142             g = d__[i__ + 1] - p;
00143             r__ = (d__[i__] - g) * s + c__ * 2. * b;
00144             p = s * r__;
00145             d__[i__ + 1] = g + p;
00146             g = c__ * r__ - b;
00147 /* L200: */
00148         }
00149 
00150         d__[l] -= p;
00151         e[l] = g;
00152         e[m] = 0.;
00153         goto L105;
00154 /*     .......... RECOVER FROM UNDERFLOW .......... */
00155 L210:
00156         d__[i__ + 1] -= p;
00157         e[m] = 0.;
00158         goto L105;
00159 /*     .......... ORDER EIGENVALUES .......... */
00160 L215:
00161         if (l == 1) {
00162             goto L250;
00163         }
00164 /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
00165         i__2 = l;
00166         for (ii = 2; ii <= i__2; ++ii) {
00167             i__ = l + 2 - ii;
00168             if (p >= d__[i__ - 1]) {
00169                 goto L270;
00170             }
00171             d__[i__] = d__[i__ - 1];
00172 /* L230: */
00173         }
00174 
00175 L250:
00176         i__ = 1;
00177 L270:
00178         d__[i__] = p;
00179 /* L290: */
00180     }
00181 
00182     goto L1001;
00183 /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
00184 /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
00185 L1000:
00186     *ierr = l;
00187 L1001:
00188     return 0;
00189 } /* imtql1_ */
00190