Doxygen Source Code Documentation
eis_hqr2.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | hqr2_ (integer *nm, integer *n, integer *low, integer *igh, doublereal *h__, doublereal *wr, doublereal *wi, doublereal *z__, integer *ierr) |
Variables | |
doublereal | c_b49 = 0. |
Function Documentation
|
Definition at line 12 of file eis_hqr2.c. References abs, c_b49, cdiv_(), d_sign(), l, max, min, p, and q. Referenced by rg_().
00015 { 00016 /* System generated locals */ 00017 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3; 00018 doublereal d__1, d__2, d__3, d__4; 00019 00020 /* Builtin functions */ 00021 double sqrt(doublereal), d_sign(doublereal *, doublereal *); 00022 00023 /* Local variables */ 00024 extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal * 00025 , doublereal *, doublereal *, doublereal *); 00026 static doublereal norm; 00027 static integer i__, j, k, l, m; 00028 static doublereal p, q, r__, s, t, w, x, y; 00029 static integer na, ii, en, jj; 00030 static doublereal ra, sa; 00031 static integer ll, mm, nn; 00032 static doublereal vi, vr, zz; 00033 static logical notlas; 00034 static integer mp2, itn, its, enm2; 00035 static doublereal tst1, tst2; 00036 00037 00038 00039 /* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2, */ 00040 /* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */ 00041 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */ 00042 00043 /* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */ 00044 /* OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE */ 00045 /* EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND */ 00046 /* IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE */ 00047 /* BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM */ 00048 /* AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS. */ 00049 00050 /* ON INPUT */ 00051 00052 /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ 00053 /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ 00054 /* DIMENSION STATEMENT. */ 00055 00056 /* N IS THE ORDER OF THE MATRIX. */ 00057 00058 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */ 00059 /* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */ 00060 /* SET LOW=1, IGH=N. */ 00061 00062 /* H CONTAINS THE UPPER HESSENBERG MATRIX. */ 00063 00064 /* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN */ 00065 /* AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE */ 00066 /* REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS */ 00067 /* OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE */ 00068 /* IDENTITY MATRIX. */ 00069 00070 /* ON OUTPUT */ 00071 00072 /* H HAS BEEN DESTROYED. */ 00073 00074 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00075 /* RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES */ 00076 /* ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS */ 00077 /* OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE */ 00078 /* HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN */ 00079 /* ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */ 00080 /* FOR INDICES IERR+1,...,N. */ 00081 00082 /* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */ 00083 /* IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z */ 00084 /* CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX 00085 */ 00086 /* WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH */ 00087 /* COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS */ 00088 /* EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN */ 00089 /* ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND. 00090 */ 00091 00092 /* IERR IS SET TO */ 00093 /* ZERO FOR NORMAL RETURN, */ 00094 /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */ 00095 /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */ 00096 00097 /* CALLS CDIV FOR COMPLEX DIVISION. */ 00098 00099 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00100 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00101 */ 00102 00103 /* THIS VERSION DATED AUGUST 1983. */ 00104 00105 /* ------------------------------------------------------------------ 00106 */ 00107 00108 /* Parameter adjustments */ 00109 z_dim1 = *nm; 00110 z_offset = z_dim1 + 1; 00111 z__ -= z_offset; 00112 --wi; 00113 --wr; 00114 h_dim1 = *nm; 00115 h_offset = h_dim1 + 1; 00116 h__ -= h_offset; 00117 00118 /* Function Body */ 00119 *ierr = 0; 00120 norm = 0.; 00121 k = 1; 00122 /* .......... STORE ROOTS ISOLATED BY BALANC */ 00123 /* AND COMPUTE MATRIX NORM .......... */ 00124 i__1 = *n; 00125 for (i__ = 1; i__ <= i__1; ++i__) { 00126 00127 i__2 = *n; 00128 for (j = k; j <= i__2; ++j) { 00129 /* L40: */ 00130 norm += (d__1 = h__[i__ + j * h_dim1], abs(d__1)); 00131 } 00132 00133 k = i__; 00134 if (i__ >= *low && i__ <= *igh) { 00135 goto L50; 00136 } 00137 wr[i__] = h__[i__ + i__ * h_dim1]; 00138 wi[i__] = 0.; 00139 L50: 00140 ; 00141 } 00142 00143 en = *igh; 00144 t = 0.; 00145 itn = *n * 30; 00146 /* .......... SEARCH FOR NEXT EIGENVALUES .......... */ 00147 L60: 00148 if (en < *low) { 00149 goto L340; 00150 } 00151 its = 0; 00152 na = en - 1; 00153 enm2 = na - 1; 00154 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */ 00155 /* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */ 00156 L70: 00157 i__1 = en; 00158 for (ll = *low; ll <= i__1; ++ll) { 00159 l = en + *low - ll; 00160 if (l == *low) { 00161 goto L100; 00162 } 00163 s = (d__1 = h__[l - 1 + (l - 1) * h_dim1], abs(d__1)) + (d__2 = h__[l 00164 + l * h_dim1], abs(d__2)); 00165 if (s == 0.) { 00166 s = norm; 00167 } 00168 tst1 = s; 00169 tst2 = tst1 + (d__1 = h__[l + (l - 1) * h_dim1], abs(d__1)); 00170 if (tst2 == tst1) { 00171 goto L100; 00172 } 00173 /* L80: */ 00174 } 00175 /* .......... FORM SHIFT .......... */ 00176 L100: 00177 x = h__[en + en * h_dim1]; 00178 if (l == en) { 00179 goto L270; 00180 } 00181 y = h__[na + na * h_dim1]; 00182 w = h__[en + na * h_dim1] * h__[na + en * h_dim1]; 00183 if (l == na) { 00184 goto L280; 00185 } 00186 if (itn == 0) { 00187 goto L1000; 00188 } 00189 if (its != 10 && its != 20) { 00190 goto L130; 00191 } 00192 /* .......... FORM EXCEPTIONAL SHIFT .......... */ 00193 t += x; 00194 00195 i__1 = en; 00196 for (i__ = *low; i__ <= i__1; ++i__) { 00197 /* L120: */ 00198 h__[i__ + i__ * h_dim1] -= x; 00199 } 00200 00201 s = (d__1 = h__[en + na * h_dim1], abs(d__1)) + (d__2 = h__[na + enm2 * 00202 h_dim1], abs(d__2)); 00203 x = s * .75; 00204 y = x; 00205 w = s * -.4375 * s; 00206 L130: 00207 ++its; 00208 --itn; 00209 /* .......... LOOK FOR TWO CONSECUTIVE SMALL */ 00210 /* SUB-DIAGONAL ELEMENTS. */ 00211 /* FOR M=EN-2 STEP -1 UNTIL L DO -- .......... */ 00212 i__1 = enm2; 00213 for (mm = l; mm <= i__1; ++mm) { 00214 m = enm2 + l - mm; 00215 zz = h__[m + m * h_dim1]; 00216 r__ = x - zz; 00217 s = y - zz; 00218 p = (r__ * s - w) / h__[m + 1 + m * h_dim1] + h__[m + (m + 1) * 00219 h_dim1]; 00220 q = h__[m + 1 + (m + 1) * h_dim1] - zz - r__ - s; 00221 r__ = h__[m + 2 + (m + 1) * h_dim1]; 00222 s = abs(p) + abs(q) + abs(r__); 00223 p /= s; 00224 q /= s; 00225 r__ /= s; 00226 if (m == l) { 00227 goto L150; 00228 } 00229 tst1 = abs(p) * ((d__1 = h__[m - 1 + (m - 1) * h_dim1], abs(d__1)) + 00230 abs(zz) + (d__2 = h__[m + 1 + (m + 1) * h_dim1], abs(d__2))); 00231 tst2 = tst1 + (d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(q) 00232 + abs(r__)); 00233 if (tst2 == tst1) { 00234 goto L150; 00235 } 00236 /* L140: */ 00237 } 00238 00239 L150: 00240 mp2 = m + 2; 00241 00242 i__1 = en; 00243 for (i__ = mp2; i__ <= i__1; ++i__) { 00244 h__[i__ + (i__ - 2) * h_dim1] = 0.; 00245 if (i__ == mp2) { 00246 goto L160; 00247 } 00248 h__[i__ + (i__ - 3) * h_dim1] = 0.; 00249 L160: 00250 ; 00251 } 00252 /* .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND */ 00253 /* COLUMNS M TO EN .......... */ 00254 i__1 = na; 00255 for (k = m; k <= i__1; ++k) { 00256 notlas = k != na; 00257 if (k == m) { 00258 goto L170; 00259 } 00260 p = h__[k + (k - 1) * h_dim1]; 00261 q = h__[k + 1 + (k - 1) * h_dim1]; 00262 r__ = 0.; 00263 if (notlas) { 00264 r__ = h__[k + 2 + (k - 1) * h_dim1]; 00265 } 00266 x = abs(p) + abs(q) + abs(r__); 00267 if (x == 0.) { 00268 goto L260; 00269 } 00270 p /= x; 00271 q /= x; 00272 r__ /= x; 00273 L170: 00274 d__1 = sqrt(p * p + q * q + r__ * r__); 00275 s = d_sign(&d__1, &p); 00276 if (k == m) { 00277 goto L180; 00278 } 00279 h__[k + (k - 1) * h_dim1] = -s * x; 00280 goto L190; 00281 L180: 00282 if (l != m) { 00283 h__[k + (k - 1) * h_dim1] = -h__[k + (k - 1) * h_dim1]; 00284 } 00285 L190: 00286 p += s; 00287 x = p / s; 00288 y = q / s; 00289 zz = r__ / s; 00290 q /= p; 00291 r__ /= p; 00292 if (notlas) { 00293 goto L225; 00294 } 00295 /* .......... ROW MODIFICATION .......... */ 00296 i__2 = *n; 00297 for (j = k; j <= i__2; ++j) { 00298 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1]; 00299 h__[k + j * h_dim1] -= p * x; 00300 h__[k + 1 + j * h_dim1] -= p * y; 00301 /* L200: */ 00302 } 00303 00304 /* Computing MIN */ 00305 i__2 = en, i__3 = k + 3; 00306 j = min(i__2,i__3); 00307 /* .......... COLUMN MODIFICATION .......... */ 00308 i__2 = j; 00309 for (i__ = 1; i__ <= i__2; ++i__) { 00310 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1]; 00311 h__[i__ + k * h_dim1] -= p; 00312 h__[i__ + (k + 1) * h_dim1] -= p * q; 00313 /* L210: */ 00314 } 00315 /* .......... ACCUMULATE TRANSFORMATIONS .......... */ 00316 i__2 = *igh; 00317 for (i__ = *low; i__ <= i__2; ++i__) { 00318 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1]; 00319 z__[i__ + k * z_dim1] -= p; 00320 z__[i__ + (k + 1) * z_dim1] -= p * q; 00321 /* L220: */ 00322 } 00323 goto L255; 00324 L225: 00325 /* .......... ROW MODIFICATION .......... */ 00326 i__2 = *n; 00327 for (j = k; j <= i__2; ++j) { 00328 p = h__[k + j * h_dim1] + q * h__[k + 1 + j * h_dim1] + r__ * h__[ 00329 k + 2 + j * h_dim1]; 00330 h__[k + j * h_dim1] -= p * x; 00331 h__[k + 1 + j * h_dim1] -= p * y; 00332 h__[k + 2 + j * h_dim1] -= p * zz; 00333 /* L230: */ 00334 } 00335 00336 /* Computing MIN */ 00337 i__2 = en, i__3 = k + 3; 00338 j = min(i__2,i__3); 00339 /* .......... COLUMN MODIFICATION .......... */ 00340 i__2 = j; 00341 for (i__ = 1; i__ <= i__2; ++i__) { 00342 p = x * h__[i__ + k * h_dim1] + y * h__[i__ + (k + 1) * h_dim1] + 00343 zz * h__[i__ + (k + 2) * h_dim1]; 00344 h__[i__ + k * h_dim1] -= p; 00345 h__[i__ + (k + 1) * h_dim1] -= p * q; 00346 h__[i__ + (k + 2) * h_dim1] -= p * r__; 00347 /* L240: */ 00348 } 00349 /* .......... ACCUMULATE TRANSFORMATIONS .......... */ 00350 i__2 = *igh; 00351 for (i__ = *low; i__ <= i__2; ++i__) { 00352 p = x * z__[i__ + k * z_dim1] + y * z__[i__ + (k + 1) * z_dim1] + 00353 zz * z__[i__ + (k + 2) * z_dim1]; 00354 z__[i__ + k * z_dim1] -= p; 00355 z__[i__ + (k + 1) * z_dim1] -= p * q; 00356 z__[i__ + (k + 2) * z_dim1] -= p * r__; 00357 /* L250: */ 00358 } 00359 L255: 00360 00361 L260: 00362 ; 00363 } 00364 00365 goto L70; 00366 /* .......... ONE ROOT FOUND .......... */ 00367 L270: 00368 h__[en + en * h_dim1] = x + t; 00369 wr[en] = h__[en + en * h_dim1]; 00370 wi[en] = 0.; 00371 en = na; 00372 goto L60; 00373 /* .......... TWO ROOTS FOUND .......... */ 00374 L280: 00375 p = (y - x) / 2.; 00376 q = p * p + w; 00377 zz = sqrt((abs(q))); 00378 h__[en + en * h_dim1] = x + t; 00379 x = h__[en + en * h_dim1]; 00380 h__[na + na * h_dim1] = y + t; 00381 if (q < 0.) { 00382 goto L320; 00383 } 00384 /* .......... REAL PAIR .......... */ 00385 zz = p + d_sign(&zz, &p); 00386 wr[na] = x + zz; 00387 wr[en] = wr[na]; 00388 if (zz != 0.) { 00389 wr[en] = x - w / zz; 00390 } 00391 wi[na] = 0.; 00392 wi[en] = 0.; 00393 x = h__[en + na * h_dim1]; 00394 s = abs(x) + abs(zz); 00395 p = x / s; 00396 q = zz / s; 00397 r__ = sqrt(p * p + q * q); 00398 p /= r__; 00399 q /= r__; 00400 /* .......... ROW MODIFICATION .......... */ 00401 i__1 = *n; 00402 for (j = na; j <= i__1; ++j) { 00403 zz = h__[na + j * h_dim1]; 00404 h__[na + j * h_dim1] = q * zz + p * h__[en + j * h_dim1]; 00405 h__[en + j * h_dim1] = q * h__[en + j * h_dim1] - p * zz; 00406 /* L290: */ 00407 } 00408 /* .......... COLUMN MODIFICATION .......... */ 00409 i__1 = en; 00410 for (i__ = 1; i__ <= i__1; ++i__) { 00411 zz = h__[i__ + na * h_dim1]; 00412 h__[i__ + na * h_dim1] = q * zz + p * h__[i__ + en * h_dim1]; 00413 h__[i__ + en * h_dim1] = q * h__[i__ + en * h_dim1] - p * zz; 00414 /* L300: */ 00415 } 00416 /* .......... ACCUMULATE TRANSFORMATIONS .......... */ 00417 i__1 = *igh; 00418 for (i__ = *low; i__ <= i__1; ++i__) { 00419 zz = z__[i__ + na * z_dim1]; 00420 z__[i__ + na * z_dim1] = q * zz + p * z__[i__ + en * z_dim1]; 00421 z__[i__ + en * z_dim1] = q * z__[i__ + en * z_dim1] - p * zz; 00422 /* L310: */ 00423 } 00424 00425 goto L330; 00426 /* .......... COMPLEX PAIR .......... */ 00427 L320: 00428 wr[na] = x + p; 00429 wr[en] = x + p; 00430 wi[na] = zz; 00431 wi[en] = -zz; 00432 L330: 00433 en = enm2; 00434 goto L60; 00435 /* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */ 00436 /* VECTORS OF UPPER TRIANGULAR FORM .......... */ 00437 L340: 00438 if (norm == 0.) { 00439 goto L1001; 00440 } 00441 /* .......... FOR EN=N STEP -1 UNTIL 1 DO -- .......... */ 00442 i__1 = *n; 00443 for (nn = 1; nn <= i__1; ++nn) { 00444 en = *n + 1 - nn; 00445 p = wr[en]; 00446 q = wi[en]; 00447 na = en - 1; 00448 if (q < 0.) { 00449 goto L710; 00450 } else if (q == 0) { 00451 goto L600; 00452 } else { 00453 goto L800; 00454 } 00455 /* .......... REAL VECTOR .......... */ 00456 L600: 00457 m = en; 00458 h__[en + en * h_dim1] = 1.; 00459 if (na == 0) { 00460 goto L800; 00461 } 00462 /* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */ 00463 i__2 = na; 00464 for (ii = 1; ii <= i__2; ++ii) { 00465 i__ = en - ii; 00466 w = h__[i__ + i__ * h_dim1] - p; 00467 r__ = 0.; 00468 00469 i__3 = en; 00470 for (j = m; j <= i__3; ++j) { 00471 /* L610: */ 00472 r__ += h__[i__ + j * h_dim1] * h__[j + en * h_dim1]; 00473 } 00474 00475 if (wi[i__] >= 0.) { 00476 goto L630; 00477 } 00478 zz = w; 00479 s = r__; 00480 goto L700; 00481 L630: 00482 m = i__; 00483 if (wi[i__] != 0.) { 00484 goto L640; 00485 } 00486 t = w; 00487 if (t != 0.) { 00488 goto L635; 00489 } 00490 tst1 = norm; 00491 t = tst1; 00492 L632: 00493 t *= .01; 00494 tst2 = norm + t; 00495 if (tst2 > tst1) { 00496 goto L632; 00497 } 00498 L635: 00499 h__[i__ + en * h_dim1] = -r__ / t; 00500 goto L680; 00501 /* .......... SOLVE REAL EQUATIONS .......... */ 00502 L640: 00503 x = h__[i__ + (i__ + 1) * h_dim1]; 00504 y = h__[i__ + 1 + i__ * h_dim1]; 00505 q = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__]; 00506 t = (x * s - zz * r__) / q; 00507 h__[i__ + en * h_dim1] = t; 00508 if (abs(x) <= abs(zz)) { 00509 goto L650; 00510 } 00511 h__[i__ + 1 + en * h_dim1] = (-r__ - w * t) / x; 00512 goto L680; 00513 L650: 00514 h__[i__ + 1 + en * h_dim1] = (-s - y * t) / zz; 00515 00516 /* .......... OVERFLOW CONTROL .......... */ 00517 L680: 00518 t = (d__1 = h__[i__ + en * h_dim1], abs(d__1)); 00519 if (t == 0.) { 00520 goto L700; 00521 } 00522 tst1 = t; 00523 tst2 = tst1 + 1. / tst1; 00524 if (tst2 > tst1) { 00525 goto L700; 00526 } 00527 i__3 = en; 00528 for (j = i__; j <= i__3; ++j) { 00529 h__[j + en * h_dim1] /= t; 00530 /* L690: */ 00531 } 00532 00533 L700: 00534 ; 00535 } 00536 /* .......... END REAL VECTOR .......... */ 00537 goto L800; 00538 /* .......... COMPLEX VECTOR .......... */ 00539 L710: 00540 m = na; 00541 /* .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT */ 00542 /* EIGENVECTOR MATRIX IS TRIANGULAR .......... */ 00543 if ((d__1 = h__[en + na * h_dim1], abs(d__1)) <= (d__2 = h__[na + en * 00544 h_dim1], abs(d__2))) { 00545 goto L720; 00546 } 00547 h__[na + na * h_dim1] = q / h__[en + na * h_dim1]; 00548 h__[na + en * h_dim1] = -(h__[en + en * h_dim1] - p) / h__[en + na * 00549 h_dim1]; 00550 goto L730; 00551 L720: 00552 d__1 = -h__[na + en * h_dim1]; 00553 d__2 = h__[na + na * h_dim1] - p; 00554 cdiv_(&c_b49, &d__1, &d__2, &q, &h__[na + na * h_dim1], &h__[na + en * 00555 h_dim1]); 00556 L730: 00557 h__[en + na * h_dim1] = 0.; 00558 h__[en + en * h_dim1] = 1.; 00559 enm2 = na - 1; 00560 if (enm2 == 0) { 00561 goto L800; 00562 } 00563 /* .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- .......... */ 00564 i__2 = enm2; 00565 for (ii = 1; ii <= i__2; ++ii) { 00566 i__ = na - ii; 00567 w = h__[i__ + i__ * h_dim1] - p; 00568 ra = 0.; 00569 sa = 0.; 00570 00571 i__3 = en; 00572 for (j = m; j <= i__3; ++j) { 00573 ra += h__[i__ + j * h_dim1] * h__[j + na * h_dim1]; 00574 sa += h__[i__ + j * h_dim1] * h__[j + en * h_dim1]; 00575 /* L760: */ 00576 } 00577 00578 if (wi[i__] >= 0.) { 00579 goto L770; 00580 } 00581 zz = w; 00582 r__ = ra; 00583 s = sa; 00584 goto L795; 00585 L770: 00586 m = i__; 00587 if (wi[i__] != 0.) { 00588 goto L780; 00589 } 00590 d__1 = -ra; 00591 d__2 = -sa; 00592 cdiv_(&d__1, &d__2, &w, &q, &h__[i__ + na * h_dim1], &h__[i__ + 00593 en * h_dim1]); 00594 goto L790; 00595 /* .......... SOLVE COMPLEX EQUATIONS .......... */ 00596 L780: 00597 x = h__[i__ + (i__ + 1) * h_dim1]; 00598 y = h__[i__ + 1 + i__ * h_dim1]; 00599 vr = (wr[i__] - p) * (wr[i__] - p) + wi[i__] * wi[i__] - q * q; 00600 vi = (wr[i__] - p) * 2. * q; 00601 if (vr != 0. || vi != 0.) { 00602 goto L784; 00603 } 00604 tst1 = norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(zz)); 00605 vr = tst1; 00606 L783: 00607 vr *= .01; 00608 tst2 = tst1 + vr; 00609 if (tst2 > tst1) { 00610 goto L783; 00611 } 00612 L784: 00613 d__1 = x * r__ - zz * ra + q * sa; 00614 d__2 = x * s - zz * sa - q * ra; 00615 cdiv_(&d__1, &d__2, &vr, &vi, &h__[i__ + na * h_dim1], &h__[i__ + 00616 en * h_dim1]); 00617 if (abs(x) <= abs(zz) + abs(q)) { 00618 goto L785; 00619 } 00620 h__[i__ + 1 + na * h_dim1] = (-ra - w * h__[i__ + na * h_dim1] + 00621 q * h__[i__ + en * h_dim1]) / x; 00622 h__[i__ + 1 + en * h_dim1] = (-sa - w * h__[i__ + en * h_dim1] - 00623 q * h__[i__ + na * h_dim1]) / x; 00624 goto L790; 00625 L785: 00626 d__1 = -r__ - y * h__[i__ + na * h_dim1]; 00627 d__2 = -s - y * h__[i__ + en * h_dim1]; 00628 cdiv_(&d__1, &d__2, &zz, &q, &h__[i__ + 1 + na * h_dim1], &h__[ 00629 i__ + 1 + en * h_dim1]); 00630 00631 /* .......... OVERFLOW CONTROL .......... */ 00632 L790: 00633 /* Computing MAX */ 00634 d__3 = (d__1 = h__[i__ + na * h_dim1], abs(d__1)), d__4 = (d__2 = 00635 h__[i__ + en * h_dim1], abs(d__2)); 00636 t = max(d__3,d__4); 00637 if (t == 0.) { 00638 goto L795; 00639 } 00640 tst1 = t; 00641 tst2 = tst1 + 1. / tst1; 00642 if (tst2 > tst1) { 00643 goto L795; 00644 } 00645 i__3 = en; 00646 for (j = i__; j <= i__3; ++j) { 00647 h__[j + na * h_dim1] /= t; 00648 h__[j + en * h_dim1] /= t; 00649 /* L792: */ 00650 } 00651 00652 L795: 00653 ; 00654 } 00655 /* .......... END COMPLEX VECTOR .......... */ 00656 L800: 00657 ; 00658 } 00659 /* .......... END BACK SUBSTITUTION. */ 00660 /* VECTORS OF ISOLATED ROOTS .......... */ 00661 i__1 = *n; 00662 for (i__ = 1; i__ <= i__1; ++i__) { 00663 if (i__ >= *low && i__ <= *igh) { 00664 goto L840; 00665 } 00666 00667 i__2 = *n; 00668 for (j = i__; j <= i__2; ++j) { 00669 /* L820: */ 00670 z__[i__ + j * z_dim1] = h__[i__ + j * h_dim1]; 00671 } 00672 00673 L840: 00674 ; 00675 } 00676 /* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */ 00677 /* VECTORS OF ORIGINAL FULL MATRIX. */ 00678 /* FOR J=N STEP -1 UNTIL LOW DO -- .......... */ 00679 i__1 = *n; 00680 for (jj = *low; jj <= i__1; ++jj) { 00681 j = *n + *low - jj; 00682 m = min(j,*igh); 00683 00684 i__2 = *igh; 00685 for (i__ = *low; i__ <= i__2; ++i__) { 00686 zz = 0.; 00687 00688 i__3 = m; 00689 for (k = *low; k <= i__3; ++k) { 00690 /* L860: */ 00691 zz += z__[i__ + k * z_dim1] * h__[k + j * h_dim1]; 00692 } 00693 00694 z__[i__ + j * z_dim1] = zz; 00695 /* L880: */ 00696 } 00697 } 00698 00699 goto L1001; 00700 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */ 00701 /* CONVERGED AFTER 30*N ITERATIONS .......... */ 00702 L1000: 00703 *ierr = en; 00704 L1001: 00705 return 0; 00706 } /* hqr2_ */ |
Variable Documentation
|
Definition at line 10 of file eis_hqr2.c. Referenced by hqr2_(). |