00001 /* tsturm.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_b26 = 1.;
00011 
00012 /* Subroutine */ int tsturm_(integer *nm, integer *n, doublereal *eps1, 
00013         doublereal *d__, doublereal *e, doublereal *e2, doublereal *lb, 
00014         doublereal *ub, integer *mm, integer *m, doublereal *w, doublereal *
00015         z__, integer *ierr, doublereal *rv1, doublereal *rv2, doublereal *rv3,
00016          doublereal *rv4, doublereal *rv5, doublereal *rv6)
00017 {
00018     /* System generated locals */
00019     integer z_dim1, z_offset, i__1, i__2, i__3;
00020     doublereal d__1, d__2, d__3, d__4;
00021 
00022     /* Builtin functions */
00023     double sqrt(doublereal);
00024 
00025     /* Local variables */
00026     static doublereal norm;
00027     static integer i__, j, k, p, q, r__, s;
00028     static doublereal u, v;
00029     static integer group, m1, m2;
00030     static doublereal t1, t2, x0, x1;
00031     static integer ii, jj, ip;
00032     static doublereal uk, xu;
00033     extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 
00034             *);
00035     static integer isturm, its;
00036     static doublereal eps2, eps3, eps4, tst1, tst2;
00037 
00038 
00039 
00040 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRISTURM */
00041 /*     BY PETERS AND WILKINSON. */
00042 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */
00043 
00044 /*     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL */
00045 /*     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL AND THEIR */
00046 /*     ASSOCIATED EIGENVECTORS, USING BISECTION AND INVERSE ITERATION. */
00047 
00048 /*     ON INPUT */
00049 
00050 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00051 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00052 /*          DIMENSION STATEMENT. */
00053 
00054 /*        N IS THE ORDER OF THE MATRIX. */
00055 
00056 /*        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED */
00057 /*          EIGENVALUES.  IT SHOULD BE CHOSEN COMMENSURATE WITH */
00058 /*          RELATIVE PERTURBATIONS IN THE MATRIX ELEMENTS OF THE */
00059 /*          ORDER OF THE RELATIVE MACHINE PRECISION.  IF THE */
00060 /*          INPUT EPS1 IS NON-POSITIVE, IT IS RESET FOR EACH */
00061 /*          SUBMATRIX TO A DEFAULT VALUE, NAMELY, MINUS THE */
00062 /*          PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE */
00063 /*          1-NORM OF THE SUBMATRIX. */
00064 
00065 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00066 
00067 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00068 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00069 
00070 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00071 /*          E2(1) IS ARBITRARY. */
00072 
00073 /*        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. */
00074 /*          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. */
00075 
00076 /*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
00077 /*          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN */
00078 /*          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, */
00079 /*          AN ERROR RETURN IS MADE WITH NO VALUES OR VECTORS FOUND. */
00080 
00081 /*     ON OUTPUT */
00082 
00083 /*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
00084 /*          (LAST) DEFAULT VALUE. */
00085 
00086 /*        D AND E ARE UNALTERED. */
00087 
00088 /*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
00089 /*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
00090 /*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
00091 /*          E2(1) IS ALSO SET TO ZERO. */
00092 
00093 /*        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). */
00094 
00095 /*        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER IF THE MATRIX */
00096 /*          DOES NOT SPLIT.  IF THE MATRIX SPLITS, THE EIGENVALUES ARE */
00097 /*          IN ASCENDING ORDER FOR EACH SUBMATRIX.  IF A VECTOR ERROR */
00098 /*          EXIT IS MADE, W CONTAINS THOSE VALUES ALREADY FOUND. */
00099 
00100 /*        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. */
00101 /*          IF AN ERROR EXIT IS MADE, Z CONTAINS THOSE VECTORS */
00102 /*          ALREADY FOUND. */
00103 
00104 /*        IERR IS SET TO */
00105 /*          ZERO       FOR NORMAL RETURN, */
00106 /*          3*N+1      IF M EXCEEDS MM. */
00107 /*          4*N+R      IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */
00108 /*                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. */
00109 
00110 /*        RV1, RV2, RV3, RV4, RV5, AND RV6 ARE TEMPORARY STORAGE ARRAYS. 
00111 */
00112 
00113 /*     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM */
00114 /*     APPEARS IN TSTURM IN-LINE. */
00115 
00116 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00117 
00118 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00119 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00120 */
00121 
00122 /*     THIS VERSION DATED AUGUST 1983. */
00123 
00124 /*     ------------------------------------------------------------------ 
00125 */
00126 
00127     /* Parameter adjustments */
00128     --rv6;
00129     --rv5;
00130     --rv4;
00131     --rv3;
00132     --rv2;
00133     --rv1;
00134     --e2;
00135     --e;
00136     --d__;
00137     z_dim1 = *nm;
00138     z_offset = z_dim1 + 1;
00139     z__ -= z_offset;
00140     --w;
00141 
00142     /* Function Body */
00143     *ierr = 0;
00144     t1 = *lb;
00145     t2 = *ub;
00146 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... */
00147     i__1 = *n;
00148     for (i__ = 1; i__ <= i__1; ++i__) {
00149         if (i__ == 1) {
00150             goto L20;
00151         }
00152         tst1 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00153                 ;
00154         tst2 = tst1 + (d__1 = e[i__], abs(d__1));
00155         if (tst2 > tst1) {
00156             goto L40;
00157         }
00158 L20:
00159         e2[i__] = 0.;
00160 L40:
00161         ;
00162     }
00163 /*     .......... DETERMINE THE NUMBER OF EIGENVALUES */
00164 /*                IN THE INTERVAL .......... */
00165     p = 1;
00166     q = *n;
00167     x1 = *ub;
00168     isturm = 1;
00169     goto L320;
00170 L60:
00171     *m = s;
00172     x1 = *lb;
00173     isturm = 2;
00174     goto L320;
00175 L80:
00176     *m -= s;
00177     if (*m > *mm) {
00178         goto L980;
00179     }
00180     q = 0;
00181     r__ = 0;
00182 /*     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING */
00183 /*                INTERVAL BY THE GERSCHGORIN BOUNDS .......... */
00184 L100:
00185     if (r__ == *m) {
00186         goto L1001;
00187     }
00188     p = q + 1;
00189     xu = d__[p];
00190     x0 = d__[p];
00191     u = 0.;
00192 
00193     i__1 = *n;
00194     for (q = p; q <= i__1; ++q) {
00195         x1 = u;
00196         u = 0.;
00197         v = 0.;
00198         if (q == *n) {
00199             goto L110;
00200         }
00201         u = (d__1 = e[q + 1], abs(d__1));
00202         v = e2[q + 1];
00203 L110:
00204 /* Computing MIN */
00205         d__1 = d__[q] - (x1 + u);
00206         xu = min(d__1,xu);
00207 /* Computing MAX */
00208         d__1 = d__[q] + (x1 + u);
00209         x0 = max(d__1,x0);
00210         if (v == 0.) {
00211             goto L140;
00212         }
00213 /* L120: */
00214     }
00215 
00216 L140:
00217 /* Computing MAX */
00218     d__2 = abs(xu), d__3 = abs(x0);
00219     d__1 = max(d__2,d__3);
00220     x1 = epslon_(&d__1);
00221     if (*eps1 <= 0.) {
00222         *eps1 = -x1;
00223     }
00224     if (p != q) {
00225         goto L180;
00226     }
00227 /*     .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... */
00228     if (t1 > d__[p] || d__[p] >= t2) {
00229         goto L940;
00230     }
00231     ++r__;
00232 
00233     i__1 = *n;
00234     for (i__ = 1; i__ <= i__1; ++i__) {
00235 /* L160: */
00236         z__[i__ + r__ * z_dim1] = 0.;
00237     }
00238 
00239     w[r__] = d__[p];
00240     z__[p + r__ * z_dim1] = 1.;
00241     goto L940;
00242 L180:
00243     u = (doublereal) (q - p + 1);
00244     x1 = u * x1;
00245 /* Computing MAX */
00246     d__1 = t1, d__2 = xu - x1;
00247     *lb = max(d__1,d__2);
00248 /* Computing MIN */
00249     d__1 = t2, d__2 = x0 + x1;
00250     *ub = min(d__1,d__2);
00251     x1 = *lb;
00252     isturm = 3;
00253     goto L320;
00254 L200:
00255     m1 = s + 1;
00256     x1 = *ub;
00257     isturm = 4;
00258     goto L320;
00259 L220:
00260     m2 = s;
00261     if (m1 > m2) {
00262         goto L940;
00263     }
00264 /*     .......... FIND ROOTS BY BISECTION .......... */
00265     x0 = *ub;
00266     isturm = 5;
00267 
00268     i__1 = m2;
00269     for (i__ = m1; i__ <= i__1; ++i__) {
00270         rv5[i__] = *ub;
00271         rv4[i__] = *lb;
00272 /* L240: */
00273     }
00274 /*     .......... LOOP FOR K-TH EIGENVALUE */
00275 /*                FOR K=M2 STEP -1 UNTIL M1 DO -- */
00276 /*                (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 
00277 */
00278     k = m2;
00279 L250:
00280     xu = *lb;
00281 /*     .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... */
00282     i__1 = k;
00283     for (ii = m1; ii <= i__1; ++ii) {
00284         i__ = m1 + k - ii;
00285         if (xu >= rv4[i__]) {
00286             goto L260;
00287         }
00288         xu = rv4[i__];
00289         goto L280;
00290 L260:
00291         ;
00292     }
00293 
00294 L280:
00295     if (x0 > rv5[k]) {
00296         x0 = rv5[k];
00297     }
00298 /*     .......... NEXT BISECTION STEP .......... */
00299 L300:
00300     x1 = (xu + x0) * .5;
00301     if (x0 - xu <= abs(*eps1)) {
00302         goto L420;
00303     }
00304     tst1 = (abs(xu) + abs(x0)) * 2.;
00305     tst2 = tst1 + (x0 - xu);
00306     if (tst2 == tst1) {
00307         goto L420;
00308     }
00309 /*     .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... */
00310 L320:
00311     s = p - 1;
00312     u = 1.;
00313 
00314     i__1 = q;
00315     for (i__ = p; i__ <= i__1; ++i__) {
00316         if (u != 0.) {
00317             goto L325;
00318         }
00319         v = (d__1 = e[i__], abs(d__1)) / epslon_(&c_b26);
00320         if (e2[i__] == 0.) {
00321             v = 0.;
00322         }
00323         goto L330;
00324 L325:
00325         v = e2[i__] / u;
00326 L330:
00327         u = d__[i__] - x1 - v;
00328         if (u < 0.) {
00329             ++s;
00330         }
00331 /* L340: */
00332     }
00333 
00334     switch (isturm) {
00335         case 1:  goto L60;
00336         case 2:  goto L80;
00337         case 3:  goto L200;
00338         case 4:  goto L220;
00339         case 5:  goto L360;
00340     }
00341 /*     .......... REFINE INTERVALS .......... */
00342 L360:
00343     if (s >= k) {
00344         goto L400;
00345     }
00346     xu = x1;
00347     if (s >= m1) {
00348         goto L380;
00349     }
00350     rv4[m1] = x1;
00351     goto L300;
00352 L380:
00353     rv4[s + 1] = x1;
00354     if (rv5[s] > x1) {
00355         rv5[s] = x1;
00356     }
00357     goto L300;
00358 L400:
00359     x0 = x1;
00360     goto L300;
00361 /*     .......... K-TH EIGENVALUE FOUND .......... */
00362 L420:
00363     rv5[k] = x1;
00364     --k;
00365     if (k >= m1) {
00366         goto L250;
00367     }
00368 /*     .......... FIND VECTORS BY INVERSE ITERATION .......... */
00369     norm = (d__1 = d__[p], abs(d__1));
00370     ip = p + 1;
00371 
00372     i__1 = q;
00373     for (i__ = ip; i__ <= i__1; ++i__) {
00374 /* L500: */
00375 /* Computing MAX */
00376         d__3 = norm, d__4 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[i__], 
00377                 abs(d__2));
00378         norm = max(d__3,d__4);
00379     }
00380 /*     .......... EPS2 IS THE CRITERION FOR GROUPING, */
00381 /*                EPS3 REPLACES ZERO PIVOTS AND EQUAL */
00382 /*                ROOTS ARE MODIFIED BY EPS3, */
00383 /*                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... */
00384     eps2 = norm * .001;
00385     eps3 = epslon_(&norm);
00386     uk = (doublereal) (q - p + 1);
00387     eps4 = uk * eps3;
00388     uk = eps4 / sqrt(uk);
00389     group = 0;
00390     s = p;
00391 
00392     i__1 = m2;
00393     for (k = m1; k <= i__1; ++k) {
00394         ++r__;
00395         its = 1;
00396         w[r__] = rv5[k];
00397         x1 = rv5[k];
00398 /*     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */
00399         if (k == m1) {
00400             goto L520;
00401         }
00402         if (x1 - x0 >= eps2) {
00403             group = -1;
00404         }
00405         ++group;
00406         if (x1 <= x0) {
00407             x1 = x0 + eps3;
00408         }
00409 /*     .......... ELIMINATION WITH INTERCHANGES AND */
00410 /*                INITIALIZATION OF VECTOR .......... */
00411 L520:
00412         v = 0.;
00413 
00414         i__2 = q;
00415         for (i__ = p; i__ <= i__2; ++i__) {
00416             rv6[i__] = uk;
00417             if (i__ == p) {
00418                 goto L560;
00419             }
00420             if ((d__1 = e[i__], abs(d__1)) < abs(u)) {
00421                 goto L540;
00422             }
00423             xu = u / e[i__];
00424             rv4[i__] = xu;
00425             rv1[i__ - 1] = e[i__];
00426             rv2[i__ - 1] = d__[i__] - x1;
00427             rv3[i__ - 1] = 0.;
00428             if (i__ != q) {
00429                 rv3[i__ - 1] = e[i__ + 1];
00430             }
00431             u = v - xu * rv2[i__ - 1];
00432             v = -xu * rv3[i__ - 1];
00433             goto L580;
00434 L540:
00435             xu = e[i__] / u;
00436             rv4[i__] = xu;
00437             rv1[i__ - 1] = u;
00438             rv2[i__ - 1] = v;
00439             rv3[i__ - 1] = 0.;
00440 L560:
00441             u = d__[i__] - x1 - xu * v;
00442             if (i__ != q) {
00443                 v = e[i__ + 1];
00444             }
00445 L580:
00446             ;
00447         }
00448 
00449         if (u == 0.) {
00450             u = eps3;
00451         }
00452         rv1[q] = u;
00453         rv2[q] = 0.;
00454         rv3[q] = 0.;
00455 /*     .......... BACK SUBSTITUTION */
00456 /*                FOR I=Q STEP -1 UNTIL P DO -- .......... */
00457 L600:
00458         i__2 = q;
00459         for (ii = p; ii <= i__2; ++ii) {
00460             i__ = p + q - ii;
00461             rv6[i__] = (rv6[i__] - u * rv2[i__] - v * rv3[i__]) / rv1[i__];
00462             v = u;
00463             u = rv6[i__];
00464 /* L620: */
00465         }
00466 /*     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */
00467 /*                MEMBERS OF GROUP .......... */
00468         if (group == 0) {
00469             goto L700;
00470         }
00471 
00472         i__2 = group;
00473         for (jj = 1; jj <= i__2; ++jj) {
00474             j = r__ - group - 1 + jj;
00475             xu = 0.;
00476 
00477             i__3 = q;
00478             for (i__ = p; i__ <= i__3; ++i__) {
00479 /* L640: */
00480                 xu += rv6[i__] * z__[i__ + j * z_dim1];
00481             }
00482 
00483             i__3 = q;
00484             for (i__ = p; i__ <= i__3; ++i__) {
00485 /* L660: */
00486                 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00487             }
00488 
00489 /* L680: */
00490         }
00491 
00492 L700:
00493         norm = 0.;
00494 
00495         i__2 = q;
00496         for (i__ = p; i__ <= i__2; ++i__) {
00497 /* L720: */
00498             norm += (d__1 = rv6[i__], abs(d__1));
00499         }
00500 
00501         if (norm >= 1.) {
00502             goto L840;
00503         }
00504 /*     .......... FORWARD SUBSTITUTION .......... */
00505         if (its == 5) {
00506             goto L960;
00507         }
00508         if (norm != 0.) {
00509             goto L740;
00510         }
00511         rv6[s] = eps4;
00512         ++s;
00513         if (s > q) {
00514             s = p;
00515         }
00516         goto L780;
00517 L740:
00518         xu = eps4 / norm;
00519 
00520         i__2 = q;
00521         for (i__ = p; i__ <= i__2; ++i__) {
00522 /* L760: */
00523             rv6[i__] *= xu;
00524         }
00525 /*     .......... ELIMINATION OPERATIONS ON NEXT VECTOR */
00526 /*                ITERATE .......... */
00527 L780:
00528         i__2 = q;
00529         for (i__ = ip; i__ <= i__2; ++i__) {
00530             u = rv6[i__];
00531 /*     .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE */
00532 /*                WAS PERFORMED EARLIER IN THE */
00533 /*                TRIANGULARIZATION PROCESS .......... */
00534             if (rv1[i__ - 1] != e[i__]) {
00535                 goto L800;
00536             }
00537             u = rv6[i__ - 1];
00538             rv6[i__ - 1] = rv6[i__];
00539 L800:
00540             rv6[i__] = u - rv4[i__] * rv6[i__ - 1];
00541 /* L820: */
00542         }
00543 
00544         ++its;
00545         goto L600;
00546 /*     .......... NORMALIZE SO THAT SUM OF SQUARES IS */
00547 /*                1 AND EXPAND TO FULL ORDER .......... */
00548 L840:
00549         u = 0.;
00550 
00551         i__2 = q;
00552         for (i__ = p; i__ <= i__2; ++i__) {
00553 /* L860: */
00554             u = pythag_(&u, &rv6[i__]);
00555         }
00556 
00557         xu = 1. / u;
00558 
00559         i__2 = *n;
00560         for (i__ = 1; i__ <= i__2; ++i__) {
00561 /* L880: */
00562             z__[i__ + r__ * z_dim1] = 0.;
00563         }
00564 
00565         i__2 = q;
00566         for (i__ = p; i__ <= i__2; ++i__) {
00567 /* L900: */
00568             z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00569         }
00570 
00571         x0 = x1;
00572 /* L920: */
00573     }
00574 
00575 L940:
00576     if (q < *n) {
00577         goto L100;
00578     }
00579     goto L1001;
00580 /*     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */
00581 L960:
00582     *ierr = (*n << 2) + r__;
00583     goto L1001;
00584 /*     .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF */
00585 /*                EIGENVALUES IN INTERVAL .......... */
00586 L980:
00587     *ierr = *n * 3 + 1;
00588 L1001:
00589     *lb = t1;
00590     *ub = t2;
00591     return 0;
00592 } /* tsturm_ */
00593