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