00001 /* bisect.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 bisect_(integer *n, doublereal *eps1, doublereal *d__, 
00013         doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, 
00014         integer *mm, 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 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 BISECTION TECHNIQUE */
00035 /*     IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. */
00036 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */
00037 
00038 /*     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL */
00039 /*     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL, */
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 /*        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. */
00061 /*          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. */
00062 
00063 /*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
00064 /*          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN */
00065 /*          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, */
00066 /*          AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND. */
00067 
00068 /*     ON OUTPUT */
00069 
00070 /*        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
00071 /*          (LAST) DEFAULT VALUE. */
00072 
00073 /*        D AND E ARE UNALTERED. */
00074 
00075 /*        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
00076 /*          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
00077 /*          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
00078 /*          E2(1) IS ALSO SET TO ZERO. */
00079 
00080 /*        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). */
00081 
00082 /*        W CONTAINS THE M EIGENVALUES 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 M EXCEEDS MM. */
00093 
00094 /*        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. */
00095 
00096 /*     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM */
00097 /*     APPEARS IN BISECT IN-LINE. */
00098 
00099 /*     NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN */
00100 /*     BISECT, 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     t1 = *lb;
00124     t2 = *ub;
00125 /*     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... */
00126     i__1 = *n;
00127     for (i__ = 1; i__ <= i__1; ++i__) {
00128         if (i__ == 1) {
00129             goto L20;
00130         }
00131         tst1 = (d__1 = d__[i__], abs(d__1)) + (d__2 = d__[i__ - 1], abs(d__2))
00132                 ;
00133         tst2 = tst1 + (d__1 = e[i__], abs(d__1));
00134         if (tst2 > tst1) {
00135             goto L40;
00136         }
00137 L20:
00138         e2[i__] = 0.;
00139 L40:
00140         ;
00141     }
00142 /*     .......... DETERMINE THE NUMBER OF EIGENVALUES */
00143 /*                IN THE INTERVAL .......... */
00144     p = 1;
00145     q = *n;
00146     x1 = *ub;
00147     isturm = 1;
00148     goto L320;
00149 L60:
00150     *m = s;
00151     x1 = *lb;
00152     isturm = 2;
00153     goto L320;
00154 L80:
00155     *m -= s;
00156     if (*m > *mm) {
00157         goto L980;
00158     }
00159     q = 0;
00160     r__ = 0;
00161 /*     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING */
00162 /*                INTERVAL BY THE GERSCHGORIN BOUNDS .......... */
00163 L100:
00164     if (r__ == *m) {
00165         goto L1001;
00166     }
00167     ++tag;
00168     p = q + 1;
00169     xu = d__[p];
00170     x0 = d__[p];
00171     u = 0.;
00172 
00173     i__1 = *n;
00174     for (q = p; q <= i__1; ++q) {
00175         x1 = u;
00176         u = 0.;
00177         v = 0.;
00178         if (q == *n) {
00179             goto L110;
00180         }
00181         u = (d__1 = e[q + 1], abs(d__1));
00182         v = e2[q + 1];
00183 L110:
00184 /* Computing MIN */
00185         d__1 = d__[q] - (x1 + u);
00186         xu = min(d__1,xu);
00187 /* Computing MAX */
00188         d__1 = d__[q] + (x1 + u);
00189         x0 = max(d__1,x0);
00190         if (v == 0.) {
00191             goto L140;
00192         }
00193 /* L120: */
00194     }
00195 
00196 L140:
00197 /* Computing MAX */
00198     d__2 = abs(xu), d__3 = abs(x0);
00199     d__1 = max(d__2,d__3);
00200     x1 = epslon_(&d__1);
00201     if (*eps1 <= 0.) {
00202         *eps1 = -x1;
00203     }
00204     if (p != q) {
00205         goto L180;
00206     }
00207 /*     .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... */
00208     if (t1 > d__[p] || d__[p] >= t2) {
00209         goto L940;
00210     }
00211     m1 = p;
00212     m2 = p;
00213     rv5[p] = d__[p];
00214     goto L900;
00215 L180:
00216     x1 *= q - p + 1;
00217 /* Computing MAX */
00218     d__1 = t1, d__2 = xu - x1;
00219     *lb = max(d__1,d__2);
00220 /* Computing MIN */
00221     d__1 = t2, d__2 = x0 + x1;
00222     *ub = min(d__1,d__2);
00223     x1 = *lb;
00224     isturm = 3;
00225     goto L320;
00226 L200:
00227     m1 = s + 1;
00228     x1 = *ub;
00229     isturm = 4;
00230     goto L320;
00231 L220:
00232     m2 = s;
00233     if (m1 > m2) {
00234         goto L940;
00235     }
00236 /*     .......... FIND ROOTS BY BISECTION .......... */
00237     x0 = *ub;
00238     isturm = 5;
00239 
00240     i__1 = m2;
00241     for (i__ = m1; i__ <= i__1; ++i__) {
00242         rv5[i__] = *ub;
00243         rv4[i__] = *lb;
00244 /* L240: */
00245     }
00246 /*     .......... LOOP FOR K-TH EIGENVALUE */
00247 /*                FOR K=M2 STEP -1 UNTIL M1 DO -- */
00248 /*                (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) .......... 
00249 */
00250     k = m2;
00251 L250:
00252     xu = *lb;
00253 /*     .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... */
00254     i__1 = k;
00255     for (ii = m1; ii <= i__1; ++ii) {
00256         i__ = m1 + k - ii;
00257         if (xu >= rv4[i__]) {
00258             goto L260;
00259         }
00260         xu = rv4[i__];
00261         goto L280;
00262 L260:
00263         ;
00264     }
00265 
00266 L280:
00267     if (x0 > rv5[k]) {
00268         x0 = rv5[k];
00269     }
00270 /*     .......... NEXT BISECTION STEP .......... */
00271 L300:
00272     x1 = (xu + x0) * .5;
00273     if (x0 - xu <= abs(*eps1)) {
00274         goto L420;
00275     }
00276     tst1 = (abs(xu) + abs(x0)) * 2.;
00277     tst2 = tst1 + (x0 - xu);
00278     if (tst2 == tst1) {
00279         goto L420;
00280     }
00281 /*     .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... */
00282 L320:
00283     s = p - 1;
00284     u = 1.;
00285 
00286     i__1 = q;
00287     for (i__ = p; i__ <= i__1; ++i__) {
00288         if (u != 0.) {
00289             goto L325;
00290         }
00291         v = (d__1 = e[i__], abs(d__1)) / epslon_(&c_b26);
00292         if (e2[i__] == 0.) {
00293             v = 0.;
00294         }
00295         goto L330;
00296 L325:
00297         v = e2[i__] / u;
00298 L330:
00299         u = d__[i__] - x1 - v;
00300         if (u < 0.) {
00301             ++s;
00302         }
00303 /* L340: */
00304     }
00305 
00306     switch (isturm) {
00307         case 1:  goto L60;
00308         case 2:  goto L80;
00309         case 3:  goto L200;
00310         case 4:  goto L220;
00311         case 5:  goto L360;
00312     }
00313 /*     .......... REFINE INTERVALS .......... */
00314 L360:
00315     if (s >= k) {
00316         goto L400;
00317     }
00318     xu = x1;
00319     if (s >= m1) {
00320         goto L380;
00321     }
00322     rv4[m1] = x1;
00323     goto L300;
00324 L380:
00325     rv4[s + 1] = x1;
00326     if (rv5[s] > x1) {
00327         rv5[s] = x1;
00328     }
00329     goto L300;
00330 L400:
00331     x0 = x1;
00332     goto L300;
00333 /*     .......... K-TH EIGENVALUE FOUND .......... */
00334 L420:
00335     rv5[k] = x1;
00336     --k;
00337     if (k >= m1) {
00338         goto L250;
00339     }
00340 /*     .......... ORDER EIGENVALUES TAGGED WITH THEIR */
00341 /*                SUBMATRIX ASSOCIATIONS .......... */
00342 L900:
00343     s = r__;
00344     r__ = r__ + m2 - m1 + 1;
00345     j = 1;
00346     k = m1;
00347 
00348     i__1 = r__;
00349     for (l = 1; l <= i__1; ++l) {
00350         if (j > s) {
00351             goto L910;
00352         }
00353         if (k > m2) {
00354             goto L940;
00355         }
00356         if (rv5[k] >= w[l]) {
00357             goto L915;
00358         }
00359 
00360         i__2 = s;
00361         for (ii = j; ii <= i__2; ++ii) {
00362             i__ = l + s - ii;
00363             w[i__ + 1] = w[i__];
00364             ind[i__ + 1] = ind[i__];
00365 /* L905: */
00366         }
00367 
00368 L910:
00369         w[l] = rv5[k];
00370         ind[l] = tag;
00371         ++k;
00372         goto L920;
00373 L915:
00374         ++j;
00375 L920:
00376         ;
00377     }
00378 
00379 L940:
00380     if (q < *n) {
00381         goto L100;
00382     }
00383     goto L1001;
00384 /*     .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF */
00385 /*                EIGENVALUES IN INTERVAL .......... */
00386 L980:
00387     *ierr = *n * 3 + 1;
00388 L1001:
00389     *lb = t1;
00390     *ub = t2;
00391     return 0;
00392 } /* bisect_ */
00393