00001 /* imtqlv.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_b11 = 1.;
00011 
00012 /* Subroutine */ int imtqlv_(integer *n, doublereal *d__, doublereal *e, 
00013         doublereal *e2, doublereal *w, integer *ind, integer *ierr, 
00014         doublereal *rv1)
00015 {
00016     /* System generated locals */
00017     integer i__1, i__2;
00018     doublereal d__1, d__2;
00019 
00020     /* Builtin functions */
00021     double d_sign(doublereal *, doublereal *);
00022 
00023     /* Local variables */
00024     static doublereal b, c__, f, g;
00025     static integer i__, j, k, l, m;
00026     static doublereal p, r__, s;
00027     static integer ii;
00028     extern doublereal pythag_(doublereal *, doublereal *);
00029     static integer tag, mml;
00030     static doublereal tst1, tst2;
00031 
00032 
00033 
00034 /*     THIS SUBROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF 
00035 */
00036 /*     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND 
00037 */
00038 /*     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. */
00039 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). */
00040 
00041 /*     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL */
00042 /*     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM */
00043 /*     THEIR CORRESPONDING SUBMATRIX INDICES. */
00044 
00045 /*     ON INPUT */
00046 
00047 /*        N IS THE ORDER OF THE MATRIX. */
00048 
00049 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00050 
00051 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00052 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00053 
00054 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00055 /*          E2(1) IS ARBITRARY. */
00056 
00057 /*     ON OUTPUT */
00058 
00059 /*        D AND E ARE UNALTERED. */
00060 
00061 /*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
00062 /*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
00063 /*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
00064 /*          E2(1) IS ALSO SET TO ZERO. */
00065 
00066 /*        W CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN */
00067 /*          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND */
00068 /*          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE */
00069 /*          THE SMALLEST EIGENVALUES. */
00070 
00071 /*        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE */
00072 /*          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES */
00073 /*          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP, */
00074 /*          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. */
00075 
00076 /*        IERR IS SET TO */
00077 /*          ZERO       FOR NORMAL RETURN, */
00078 /*          J          IF THE J-TH EIGENVALUE HAS NOT BEEN */
00079 /*                     DETERMINED AFTER 30 ITERATIONS. */
00080 
00081 /*        RV1 IS A TEMPORARY STORAGE ARRAY. */
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     --rv1;
00096     --ind;
00097     --w;
00098     --e2;
00099     --e;
00100     --d__;
00101 
00102     /* Function Body */
00103     *ierr = 0;
00104     k = 0;
00105     tag = 0;
00106 
00107     i__1 = *n;
00108     for (i__ = 1; i__ <= i__1; ++i__) {
00109         w[i__] = d__[i__];
00110         if (i__ != 1) {
00111             rv1[i__ - 1] = e[i__];
00112         }
00113 /* L100: */
00114     }
00115 
00116     e2[1] = 0.;
00117     rv1[*n] = 0.;
00118 
00119     i__1 = *n;
00120     for (l = 1; l <= i__1; ++l) {
00121         j = 0;
00122 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */
00123 L105:
00124         i__2 = *n;
00125         for (m = l; m <= i__2; ++m) {
00126             if (m == *n) {
00127                 goto L120;
00128             }
00129             tst1 = (d__1 = w[m], abs(d__1)) + (d__2 = w[m + 1], abs(d__2));
00130             tst2 = tst1 + (d__1 = rv1[m], abs(d__1));
00131             if (tst2 == tst1) {
00132                 goto L120;
00133             }
00134 /*     .......... GUARD AGAINST UNDERFLOWED ELEMENT OF E2 ........
00135 .. */
00136             if (e2[m + 1] == 0.) {
00137                 goto L125;
00138             }
00139 /* L110: */
00140         }
00141 
00142 L120:
00143         if (m <= k) {
00144             goto L130;
00145         }
00146         if (m != *n) {
00147             e2[m + 1] = 0.;
00148         }
00149 L125:
00150         k = m;
00151         ++tag;
00152 L130:
00153         p = w[l];
00154         if (m == l) {
00155             goto L215;
00156         }
00157         if (j == 30) {
00158             goto L1000;
00159         }
00160         ++j;
00161 /*     .......... FORM SHIFT .......... */
00162         g = (w[l + 1] - p) / (rv1[l] * 2.);
00163         r__ = pythag_(&g, &c_b11);
00164         g = w[m] - p + rv1[l] / (g + d_sign(&r__, &g));
00165         s = 1.;
00166         c__ = 1.;
00167         p = 0.;
00168         mml = m - l;
00169 /*     .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */
00170         i__2 = mml;
00171         for (ii = 1; ii <= i__2; ++ii) {
00172             i__ = m - ii;
00173             f = s * rv1[i__];
00174             b = c__ * rv1[i__];
00175             r__ = pythag_(&f, &g);
00176             rv1[i__ + 1] = r__;
00177             if (r__ == 0.) {
00178                 goto L210;
00179             }
00180             s = f / r__;
00181             c__ = g / r__;
00182             g = w[i__ + 1] - p;
00183             r__ = (w[i__] - g) * s + c__ * 2. * b;
00184             p = s * r__;
00185             w[i__ + 1] = g + p;
00186             g = c__ * r__ - b;
00187 /* L200: */
00188         }
00189 
00190         w[l] -= p;
00191         rv1[l] = g;
00192         rv1[m] = 0.;
00193         goto L105;
00194 /*     .......... RECOVER FROM UNDERFLOW .......... */
00195 L210:
00196         w[i__ + 1] -= p;
00197         rv1[m] = 0.;
00198         goto L105;
00199 /*     .......... ORDER EIGENVALUES .......... */
00200 L215:
00201         if (l == 1) {
00202             goto L250;
00203         }
00204 /*     .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... */
00205         i__2 = l;
00206         for (ii = 2; ii <= i__2; ++ii) {
00207             i__ = l + 2 - ii;
00208             if (p >= w[i__ - 1]) {
00209                 goto L270;
00210             }
00211             w[i__] = w[i__ - 1];
00212             ind[i__] = ind[i__ - 1];
00213 /* L230: */
00214         }
00215 
00216 L250:
00217         i__ = 1;
00218 L270:
00219         w[i__] = p;
00220         ind[i__] = tag;
00221 /* L290: */
00222     }
00223 
00224     goto L1001;
00225 /*     .......... SET ERROR -- NO CONVERGENCE TO AN */
00226 /*                EIGENVALUE AFTER 30 ITERATIONS .......... */
00227 L1000:
00228     *ierr = l;
00229 L1001:
00230     return 0;
00231 } /* imtqlv_ */
00232