Doxygen Source Code Documentation
eis_bisect.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | bisect_ (integer *n, doublereal *eps1, doublereal *d__, doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, integer *mm, integer *m, doublereal *w, integer *ind, integer *ierr, doublereal *rv4, doublereal *rv5) |
Variables | |
doublereal | c_b26 = 1. |
Function Documentation
|
Definition at line 12 of file eis_bisect.c. References abs, c_b26, epslon_(), ind, l, m1, m2, max, min, p, q, v, and x0.
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_ */ |
Variable Documentation
|
Definition at line 10 of file eis_bisect.c. Referenced by bisect_(). |