Doxygen Source Code Documentation
eis_tridib.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | tridib_ (integer *n, doublereal *eps1, doublereal *d__, doublereal *e, doublereal *e2, doublereal *lb, doublereal *ub, integer *m11, integer *m, doublereal *w, integer *ind, integer *ierr, doublereal *rv4, doublereal *rv5) |
Variables | |
doublereal | c_b33 = 1. |
Function Documentation
|
Definition at line 12 of file eis_tridib.c. References abs, c_b33, 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 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_ */ |
Variable Documentation
|
Definition at line 10 of file eis_tridib.c. Referenced by tridib_(). |