00011 {
00012     /* System generated locals */
00013     integer i__1, i__2;
00014     doublereal d__1, d__2, d__3;
00015 
00016     /* Local variables */
00017     static integer jdef;
00018     static doublereal f;
00019     static integer i__, j, k;
00020     static doublereal p, q, r__, s, delta;
00021     static integer k1, ii, jj;
00022     static doublereal ep, qp;
00023     extern doublereal epslon_(doublereal *);
00024     static doublereal err, tot;
00025 
00026 
00027 
00028 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE RATQR, */
00029 /*     NUM. MATH. 11, 264-272(1968) BY REINSCH AND BAUER. */
00030 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 257-265(1971). */
00031 
00032 /*     THIS SUBROUTINE FINDS THE ALGEBRAICALLY SMALLEST OR LARGEST */
00033 /*     EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE */
00034 /*     RATIONAL QR METHOD WITH NEWTON CORRECTIONS. */
00035 
00036 /*     ON INPUT */
00037 
00038 /*        N IS THE ORDER OF THE MATRIX. */
00039 
00040 /*        EPS1 IS A THEORETICAL ABSOLUTE ERROR TOLERANCE FOR THE */
00041 /*          COMPUTED EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE, */
00042 /*          OR INDEED SMALLER THAN ITS DEFAULT VALUE, IT IS RESET */
00043 /*          AT EACH ITERATION TO THE RESPECTIVE DEFAULT VALUE, */
00044 /*          NAMELY, THE PRODUCT OF THE RELATIVE MACHINE PRECISION */
00045 /*          AND THE MAGNITUDE OF THE CURRENT EIGENVALUE ITERATE. */
00046 /*          THE THEORETICAL ABSOLUTE ERROR IN THE K-TH EIGENVALUE */
00047 /*          IS USUALLY NOT GREATER THAN K TIMES EPS1. */
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 /*        M IS THE NUMBER OF EIGENVALUES TO BE FOUND. */
00058 
00059 /*        IDEF SHOULD BE SET TO 1 IF THE INPUT MATRIX IS KNOWN TO BE */
00060 /*          POSITIVE DEFINITE, TO -1 IF THE INPUT MATRIX IS KNOWN TO */
00061 /*          BE NEGATIVE DEFINITE, AND TO 0 OTHERWISE. */
00062 
00063 /*        TYPE SHOULD BE SET TO .TRUE. IF THE SMALLEST EIGENVALUES */
00064 /*          ARE TO BE FOUND, AND TO .FALSE. IF THE LARGEST EIGENVALUES */
00065 /*          ARE TO BE FOUND. */
00066 
00067 /*     ON OUTPUT */
00068 
00069 /*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
00070 /*          (LAST) DEFAULT VALUE. */
00071 
00072 /*        D AND E ARE UNALTERED (UNLESS W OVERWRITES D). */
00073 
00074 /*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
00075 /*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
00076 /*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
00077 /*          E2(1) IS SET TO 0.0D0 IF THE SMALLEST EIGENVALUES HAVE BEEN */
00078 /*          FOUND, AND TO 2.0D0 IF THE LARGEST EIGENVALUES HAVE BEEN */
00079 /*          FOUND.  E2 IS OTHERWISE UNALTERED (UNLESS OVERWRITTEN BY BD). 
00080 */
00081 
00082 /*        W CONTAINS THE M ALGEBRAICALLY SMALLEST EIGENVALUES IN */
00083 /*          ASCENDING ORDER, OR THE M LARGEST EIGENVALUES IN */
00084 /*          DESCENDING ORDER.  IF AN ERROR EXIT IS MADE BECAUSE OF */
00085 /*          AN INCORRECT SPECIFICATION OF IDEF, NO EIGENVALUES */
00086 /*          ARE FOUND.  IF THE NEWTON ITERATES FOR A PARTICULAR */
00087 /*          EIGENVALUE ARE NOT MONOTONE, THE BEST ESTIMATE OBTAINED */
00088 /*          IS RETURNED AND IERR IS SET.  W MAY COINCIDE WITH D. */
00089 
00090 /*        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
00091 /*          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
00092 /*          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
00093 /*          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. 
00094 */
00095 
00096 /*        BD CONTAINS REFINED BOUNDS FOR THE THEORETICAL ERRORS OF THE */
00097 /*          CORRESPONDING EIGENVALUES IN W.  THESE BOUNDS ARE USUALLY */
00098 /*          WITHIN THE TOLERANCE SPECIFIED BY EPS1.  BD MAY COINCIDE */
00099 /*          WITH E2. */
00100 
00101 /*        IERR IS SET TO */
00102 /*          ZERO       FOR NORMAL RETURN, */
00103 /*          6*N+1      IF  IDEF  IS SET TO 1 AND  TYPE  TO .TRUE. */
00104 /*                     WHEN THE MATRIX IS NOT POSITIVE DEFINITE, OR */
00105 /*                     IF  IDEF  IS SET TO -1 AND  TYPE  TO .FALSE. */
00106 /*                     WHEN THE MATRIX IS NOT NEGATIVE DEFINITE, */
00107 /*          5*N+K      IF SUCCESSIVE ITERATES TO THE K-TH EIGENVALUE */
00108 /*                     ARE NOT MONOTONE INCREASING, WHERE K REFERS */
00109 /*                     TO THE LAST SUCH OCCURRENCE. */
00110 
00111 /*     NOTE THAT SUBROUTINE TRIDIB IS GENERALLY FASTER AND MORE */
00112 /*     ACCURATE THAN RATQR IF THE EIGENVALUES ARE CLUSTERED. */
00113 
00114 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00115 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00116 */
00117 
00118 /*     THIS VERSION DATED AUGUST 1983. */
00119 
00120 /*     ------------------------------------------------------------------ 
00121 */
00122 
00123     /* Parameter adjustments */
00124     --bd;
00125     --ind;
00126     --w;
00127     --e2;
00128     --e;
00129     --d__;
00130 
00131     /* Function Body */
00132     *ierr = 0;
00133     jdef = *idef;
00134 /*     .......... COPY D ARRAY INTO W .......... */
00135     i__1 = *n;
00136     for (i__ = 1; i__ <= i__1; ++i__) {
00137 /* L20: */
00138         w[i__] = d__[i__];
00139     }
00140 
00141     if (*type__) {
00142         goto L40;
00143     }
00144     j = 1;
00145     goto L400;
00146 L40:
00147     err = 0.;
00148     s = 0.;
00149 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DEFINE */
00150 /*                INITIAL SHIFT FROM LOWER GERSCHGORIN BOUND. */
00151 /*                COPY E2 ARRAY INTO BD .......... */
00152     tot = w[1];
00153     q = 0.;
00154     j = 0;
00155 
00156     i__1 = *n;
00157     for (i__ = 1; i__ <= i__1; ++i__) {
00158         p = q;
00159         if (i__ == 1) {
00160             goto L60;
00161         }
00162         d__3 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00163                 ;
00164         if (p > epslon_(&d__3)) {
00165             goto L80;
00166         }
00167 L60:
00168         e2[i__] = 0.;
00169 L80:
00170         bd[i__] = e2[i__];
00171 /*     .......... COUNT ALSO IF ELEMENT OF E2 HAS UNDERFLOWED ........
00172 .. */
00173         if (e2[i__] == 0.) {
00174             ++j;
00175         }
00176         ind[i__] = j;
00177         q = 0.;
00178         if (i__ != *n) {
00179             q = (d__1 = e[i__ + 1], abs(d__1));
00180         }
00181 /* Computing MIN */
00182         d__1 = w[i__] - p - q;
00183         tot = min(d__1,tot);
00184 /* L100: */
00185     }
00186 
00187     if (jdef == 1 && tot < 0.) {
00188         goto L140;
00189     }
00190 
00191     i__1 = *n;
00192     for (i__ = 1; i__ <= i__1; ++i__) {
00193 /* L110: */
00194         w[i__] -= tot;
00195     }
00196 
00197     goto L160;
00198 L140:
00199     tot = 0.;
00200 
00201 L160:
00202     i__1 = *m;
00203     for (k = 1; k <= i__1; ++k) {
00204 /*     .......... NEXT QR TRANSFORMATION .......... */
00205 L180:
00206         tot += s;
00207         delta = w[*n] - s;
00208         i__ = *n;
00209         f = (d__1 = epslon_(&tot), abs(d__1));
00210         if (*eps1 < f) {
00211             *eps1 = f;
00212         }
00213         if (delta > *eps1) {
00214             goto L190;
00215         }
00216         if (delta < -(*eps1)) {
00217             goto L1000;
00218         }
00219         goto L300;
00220 /*     .......... REPLACE SMALL SUB-DIAGONAL SQUARES BY ZERO */
00221 /*                TO REDUCE THE INCIDENCE OF UNDERFLOWS .......... */
00222 L190:
00223         if (k == *n) {
00224             goto L210;
00225         }
00226         k1 = k + 1;
00227         i__2 = *n;
00228         for (j = k1; j <= i__2; ++j) {
00229             d__2 = w[j] + w[j - 1];
00230 /* Computing 2nd power */
00231             d__1 = epslon_(&d__2);
00232             if (bd[j] <= d__1 * d__1) {
00233                 bd[j] = 0.;
00234             }
00235 /* L200: */
00236         }
00237 
00238 L210:
00239         f = bd[*n] / delta;
00240         qp = delta + f;
00241         p = 1.;
00242         if (k == *n) {
00243             goto L260;
00244         }
00245         k1 = *n - k;
00246 /*     .......... FOR I=N-1 STEP -1 UNTIL K DO -- .......... */
00247         i__2 = k1;
00248         for (ii = 1; ii <= i__2; ++ii) {
00249             i__ = *n - ii;
00250             q = w[i__] - s - f;
00251             r__ = q / qp;
00252             p = p * r__ + 1.;
00253             ep = f * r__;
00254             w[i__ + 1] = qp + ep;
00255             delta = q - ep;
00256             if (delta > *eps1) {
00257                 goto L220;
00258             }
00259             if (delta < -(*eps1)) {
00260                 goto L1000;
00261             }
00262             goto L300;
00263 L220:
00264             f = bd[i__] / q;
00265             qp = delta + f;
00266             bd[i__ + 1] = qp * ep;
00267 /* L240: */
00268         }
00269 
00270 L260:
00271         w[k] = qp;
00272         s = qp / p;
00273         if (tot + s > tot) {
00274             goto L180;
00275         }
00276 /*     .......... SET ERROR -- IRREGULAR END OF ITERATION. */
00277 /*                DEFLATE MINIMUM DIAGONAL ELEMENT .......... */
00278         *ierr = *n * 5 + k;
00279         s = 0.;
00280         delta = qp;
00281 
00282         i__2 = *n;
00283         for (j = k; j <= i__2; ++j) {
00284             if (w[j] > delta) {
00285                 goto L280;
00286             }
00287             i__ = j;
00288             delta = w[j];
00289 L280:
00290             ;
00291         }
00292 /*     .......... CONVERGENCE .......... */
00293 L300:
00294         if (i__ < *n) {
00295             bd[i__ + 1] = bd[i__] * f / qp;
00296         }
00297         ii = ind[i__];
00298         if (i__ == k) {
00299             goto L340;
00300         }
00301         k1 = i__ - k;
00302 /*     .......... FOR J=I-1 STEP -1 UNTIL K DO -- .......... */
00303         i__2 = k1;
00304         for (jj = 1; jj <= i__2; ++jj) {
00305             j = i__ - jj;
00306             w[j + 1] = w[j] - s;
00307             bd[j + 1] = bd[j];
00308             ind[j + 1] = ind[j];
00309 /* L320: */
00310         }
00311 
00312 L340:
00313         w[k] = tot;
00314         err += abs(delta);
00315         bd[k] = err;
00316         ind[k] = ii;
00317 /* L360: */
00318     }
00319 
00320     if (*type__) {
00321         goto L1001;
00322     }
00323     f = bd[1];
00324     e2[1] = 2.;
00325     bd[1] = f;
00326     j = 2;
00327 /*     .......... NEGATE ELEMENTS OF W FOR LARGEST VALUES .......... */
00328 L400:
00329     i__1 = *n;
00330     for (i__ = 1; i__ <= i__1; ++i__) {
00331 /* L500: */
00332         w[i__] = -w[i__];
00333     }
00334 
00335     jdef = -jdef;
00336     switch (j) {
00337         case 1:  goto L40;
00338         case 2:  goto L1001;
00339     }
00340 /*     .......... SET ERROR -- IDEF SPECIFIED INCORRECTLY .......... */
00341 L1000:
00342     *ierr = *n * 6 + 1;
00343 L1001:
00344     return 0;
00345 } /* ratqr_ */