00001 /* tridib.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_b33 = 1.;
00011 
00012 /* Subroutine */ int tridib_(integer *n, doublereal *eps1, doublereal *d__, 
00013         doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, 
00014         integer *m11, integer *m, doublereal *w, integer *ind, integer *ierr, 
00015         doublereal *rv4, doublereal *rv5)
00016 {
00017     /* System generated locals */
00018     integer i__1, i__2;
00019     doublereal d__1, d__2, d__3;
00020 
00021     /* Local variables */
00022     static integer i__, j, k, l, p, q, r__, s;
00023     static doublereal u, v;
00024     static integer m1, m2;
00025     static doublereal t1, t2, x0, x1;
00026     static integer m22, ii;
00027     static doublereal xu;
00028     extern doublereal epslon_(doublereal *);
00029     static integer isturm, tag;
00030     static doublereal tst1, tst2;
00031 
00032 
00033 
00034 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BISECT, */
00035 /*     NUM. MATH. 9, 386-393(1967) BY BARTH, MARTIN, AND WILKINSON. */
00036 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 249-256(1971). */
00037 
00038 /*     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL */
00039 /*     SYMMETRIC MATRIX BETWEEN SPECIFIED BOUNDARY INDICES, */
00040 /*     USING BISECTION. */
00041 
00042 /*     ON INPUT */
00043 
00044 /*        N IS THE ORDER OF THE MATRIX. */
00045 
00046 /*        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED */
00047 /*          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE, */
00048 /*          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, */
00049 /*          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE */
00050 /*          PRECISION AND THE 1-NORM OF THE SUBMATRIX. */
00051 
00052 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00053 
00054 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00055 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00056 
00057 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
00058 /*          E2(1) IS ARBITRARY. */
00059 
00060 /*        M11 SPECIFIES THE LOWER BOUNDARY INDEX FOR THE DESIRED */
00061 /*          EIGENVALUES. */
00062 
00063 /*        M SPECIFIES THE NUMBER OF EIGENVALUES DESIRED.  THE UPPER */
00064 /*          BOUNDARY INDEX M22 IS THEN OBTAINED AS M22=M11+M-1. */
00065 
00066 /*     ON OUTPUT */
00067 
00068 /*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
00069 /*          (LAST) DEFAULT VALUE. */
00070 
00071 /*        D AND E ARE UNALTERED. */
00072 
00073 /*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
00074 /*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
00075 /*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
00076 /*          E2(1) IS ALSO SET TO ZERO. */
00077 
00078 /*        LB AND UB DEFINE AN INTERVAL CONTAINING EXACTLY THE DESIRED */
00079 /*          EIGENVALUES. */
00080 
00081 /*        W CONTAINS, IN ITS FIRST M POSITIONS, THE EIGENVALUES */
00082 /*          BETWEEN INDICES M11 AND M22 IN ASCENDING ORDER. */
00083 
00084 /*        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
00085 /*          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
00086 /*          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
00087 /*          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.. 
00088 */
00089 
00090 /*        IERR IS SET TO */
00091 /*          ZERO       FOR NORMAL RETURN, */
00092 /*          3*N+1      IF MULTIPLE EIGENVALUES AT INDEX M11 MAKE */
00093 /*                     UNIQUE SELECTION IMPOSSIBLE, */
00094 /*          3*N+2      IF MULTIPLE EIGENVALUES AT INDEX M22 MAKE */
00095 /*                     UNIQUE SELECTION IMPOSSIBLE. */
00096 
00097 /*        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. */
00098 
00099 /*     NOTE THAT SUBROUTINE TQL1, IMTQL1, OR TQLRAT IS GENERALLY FASTER */
00100 /*     THAN TRIDIB, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. */
00101 
00102 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00103 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00104 */
00105 
00106 /*     THIS VERSION DATED AUGUST 1983. */
00107 
00108 /*     ------------------------------------------------------------------ 
00109 */
00110 
00111     /* Parameter adjustments */
00112     --rv5;
00113     --rv4;
00114     --e2;
00115     --e;
00116     --d__;
00117     --ind;
00118     --w;
00119 
00120     /* Function Body */
00121     *ierr = 0;
00122     tag = 0;
00123     xu = d__[1];
00124     x0 = d__[1];
00125     u = 0.;
00126 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES AND DETERMINE AN */
00127 /*                INTERVAL CONTAINING ALL THE EIGENVALUES .......... */
00128     i__1 = *n;
00129     for (i__ = 1; i__ <= i__1; ++i__) {
00130         x1 = u;
00131         u = 0.;
00132         if (i__ != *n) {
00133             u = (d__1 = e[i__ + 1], abs(d__1));
00134         }
00135 /* Computing MIN */
00136         d__1 = d__[i__] - (x1 + u);
00137         xu = min(d__1,xu);
00138 /* Computing MAX */
00139         d__1 = d__[i__] + (x1 + u);
00140         x0 = max(d__1,x0);
00141         if (i__ == 1) {
00142             goto L20;
00143         }
00144         tst1 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00145                 ;
00146         tst2 = tst1 + (d__1 = e[i__], abs(d__1));
00147         if (tst2 > tst1) {
00148             goto L40;
00149         }
00150 L20:
00151         e2[i__] = 0.;
00152 L40:
00153         ;
00154     }
00155 
00156     x1 = (doublereal) (*n);
00157 /* Computing MAX */
00158     d__2 = abs(xu), d__3 = abs(x0);
00159     d__1 = max(d__2,d__3);
00160     x1 *= epslon_(&d__1);
00161     xu -= x1;
00162     t1 = xu;
00163     x0 += x1;
00164     t2 = x0;
00165 /*     .......... DETERMINE AN INTERVAL CONTAINING EXACTLY */
00166 /*                THE DESIRED EIGENVALUES .......... */
00167     p = 1;
00168     q = *n;
00169     m1 = *m11 - 1;
00170     if (m1 == 0) {
00171         goto L75;
00172     }
00173     isturm = 1;
00174 L50:
00175     v = x1;
00176     x1 = xu + (x0 - xu) * .5;
00177     if (x1 == v) {
00178         goto L980;
00179     }
00180     goto L320;
00181 L60:
00182     if ((i__1 = s - m1) < 0) {
00183         goto L65;
00184     } else if (i__1 == 0) {
00185         goto L73;
00186     } else {
00187         goto L70;
00188     }
00189 L65:
00190     xu = x1;
00191     goto L50;
00192 L70:
00193     x0 = x1;
00194     goto L50;
00195 L73:
00196     xu = x1;
00197     t1 = x1;
00198 L75:
00199     m22 = m1 + *m;
00200     if (m22 == *n) {
00201         goto L90;
00202     }
00203     x0 = t2;
00204     isturm = 2;
00205     goto L50;
00206 L80:
00207     if ((i__1 = s - m22) < 0) {
00208         goto L65;
00209     } else if (i__1 == 0) {
00210         goto L85;
00211     } else {
00212         goto L70;
00213     }
00214 L85:
00215     t2 = x1;
00216 L90:
00217     q = 0;
00218     r__ = 0;
00219 /*     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING */
00220 /*                INTERVAL BY THE GERSCHGORIN BOUNDS .......... */
00221 L100:
00222     if (r__ == *m) {
00223         goto L1001;
00224     }
00225     ++tag;
00226     p = q + 1;
00227     xu = d__[p];
00228     x0 = d__[p];
00229     u = 0.;
00230 
00231     i__1 = *n;
00232     for (q = p; q <= i__1; ++q) {
00233         x1 = u;
00234         u = 0.;
00235         v = 0.;
00236         if (q == *n) {
00237             goto L110;
00238         }
00239         u = (d__1 = e[q + 1], abs(d__1));
00240         v = e2[q + 1];
00241 L110:
00242 /* Computing MIN */
00243         d__1 = d__[q] - (x1 + u);
00244         xu = min(d__1,xu);
00245 /* Computing MAX */
00246         d__1 = d__[q] + (x1 + u);
00247         x0 = max(d__1,x0);
00248         if (v == 0.) {
00249             goto L140;
00250         }
00251 /* L120: */
00252     }
00253 
00254 L140:
00255 /* Computing MAX */
00256     d__2 = abs(xu), d__3 = abs(x0);
00257     d__1 = max(d__2,d__3);
00258     x1 = epslon_(&d__1);
00259     if (*eps1 <= 0.) {
00260         *eps1 = -x1;
00261     }
00262     if (p != q) {
00263         goto L180;
00264     }
00265 /*     .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... */
00266     if (t1 > d__[p] || d__[p] >= t2) {
00267         goto L940;
00268     }
00269     m1 = p;
00270     m2 = p;
00271     rv5[p] = d__[p];
00272     goto L900;
00273 L180:
00274     x1 *= q - p + 1;
00275 /* Computing MAX */
00276     d__1 = t1, d__2 = xu - x1;
00277     *lb = max(d__1,d__2);
00278 /* Computing MIN */
00279     d__1 = t2, d__2 = x0 + x1;
00280     *ub = min(d__1,d__2);
00281     x1 = *lb;
00282     isturm = 3;
00283     goto L320;
00284 L200:
00285     m1 = s + 1;
00286     x1 = *ub;
00287     isturm = 4;
00288     goto L320;
00289 L220:
00290     m2 = s;
00291     if (m1 > m2) {
00292         goto L940;
00293     }
00294 /*     .......... FIND ROOTS BY BISECTION .......... */
00295     x0 = *ub;
00296     isturm = 5;
00297 
00298     i__1 = m2;
00299     for (i__ = m1; i__ <= i__1; ++i__) {
00300         rv5[i__] = *ub;
00301         rv4[i__] = *lb;
00302 /* L240: */
00303     }
00304 /*     .......... LOOP FOR K-TH EIGENVALUE */
00305 /*                FOR K=M2 STEP -1 UNTIL M1 DO -- */
00306 /*                (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 
00307 */
00308     k = m2;
00309 L250:
00310     xu = *lb;
00311 /*     .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... */
00312     i__1 = k;
00313     for (ii = m1; ii <= i__1; ++ii) {
00314         i__ = m1 + k - ii;
00315         if (xu >= rv4[i__]) {
00316             goto L260;
00317         }
00318         xu = rv4[i__];
00319         goto L280;
00320 L260:
00321         ;
00322     }
00323 
00324 L280:
00325     if (x0 > rv5[k]) {
00326         x0 = rv5[k];
00327     }
00328 /*     .......... NEXT BISECTION STEP .......... */
00329 L300:
00330     x1 = (xu + x0) * .5;
00331     if (x0 - xu <= abs(*eps1)) {
00332         goto L420;
00333     }
00334     tst1 = (abs(xu) + abs(x0)) * 2.;
00335     tst2 = tst1 + (x0 - xu);
00336     if (tst2 == tst1) {
00337         goto L420;
00338     }
00339 /*     .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... */
00340 L320:
00341     s = p - 1;
00342     u = 1.;
00343 
00344     i__1 = q;
00345     for (i__ = p; i__ <= i__1; ++i__) {
00346         if (u != 0.) {
00347             goto L325;
00348         }
00349         v = (d__1 = e[i__], abs(d__1)) / epslon_(&c_b33);
00350         if (e2[i__] == 0.) {
00351             v = 0.;
00352         }
00353         goto L330;
00354 L325:
00355         v = e2[i__] / u;
00356 L330:
00357         u = d__[i__] - x1 - v;
00358         if (u < 0.) {
00359             ++s;
00360         }
00361 /* L340: */
00362     }
00363 
00364     switch (isturm) {
00365         case 1:  goto L60;
00366         case 2:  goto L80;
00367         case 3:  goto L200;
00368         case 4:  goto L220;
00369         case 5:  goto L360;
00370     }
00371 /*     .......... REFINE INTERVALS .......... */
00372 L360:
00373     if (s >= k) {
00374         goto L400;
00375     }
00376     xu = x1;
00377     if (s >= m1) {
00378         goto L380;
00379     }
00380     rv4[m1] = x1;
00381     goto L300;
00382 L380:
00383     rv4[s + 1] = x1;
00384     if (rv5[s] > x1) {
00385         rv5[s] = x1;
00386     }
00387     goto L300;
00388 L400:
00389     x0 = x1;
00390     goto L300;
00391 /*     .......... K-TH EIGENVALUE FOUND .......... */
00392 L420:
00393     rv5[k] = x1;
00394     --k;
00395     if (k >= m1) {
00396         goto L250;
00397     }
00398 /*     .......... ORDER EIGENVALUES TAGGED WITH THEIR */
00399 /*                SUBMATRIX ASSOCIATIONS .......... */
00400 L900:
00401     s = r__;
00402     r__ = r__ + m2 - m1 + 1;
00403     j = 1;
00404     k = m1;
00405 
00406     i__1 = r__;
00407     for (l = 1; l <= i__1; ++l) {
00408         if (j > s) {
00409             goto L910;
00410         }
00411         if (k > m2) {
00412             goto L940;
00413         }
00414         if (rv5[k] >= w[l]) {
00415             goto L915;
00416         }
00417 
00418         i__2 = s;
00419         for (ii = j; ii <= i__2; ++ii) {
00420             i__ = l + s - ii;
00421             w[i__ + 1] = w[i__];
00422             ind[i__ + 1] = ind[i__];
00423 /* L905: */
00424         }
00425 
00426 L910:
00427         w[l] = rv5[k];
00428         ind[l] = tag;
00429         ++k;
00430         goto L920;
00431 L915:
00432         ++j;
00433 L920:
00434         ;
00435     }
00436 
00437 L940:
00438     if (q < *n) {
00439         goto L100;
00440     }
00441     goto L1001;
00442 /*     .......... SET ERROR -- INTERVAL CANNOT BE FOUND CONTAINING */
00443 /*                EXACTLY THE DESIRED EIGENVALUES .......... */
00444 L980:
00445     *ierr = *n * 3 + isturm;
00446 L1001:
00447     *lb = t1;
00448     *ub = t2;
00449     return 0;
00450 } /* tridib_ */
00451