Doxygen Source Code Documentation
eis_comqr2.c File Reference
#include "f2c.h"
Go to the source code of this file.
Functions | |
int | comqr2_ (integer *nm, integer *n, integer *low, integer *igh, doublereal *ortr, doublereal *orti, doublereal *hr, doublereal *hi, doublereal *wr, doublereal *wi, doublereal *zr, doublereal *zi, integer *ierr) |
Function Documentation
|
Definition at line 8 of file eis_comqr2.c. References abs, cdiv_(), csroot_(), l, min, and pythag_(). Referenced by cg_().
00012 { 00013 /* System generated locals */ 00014 integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset, 00015 zi_dim1, zi_offset, i__1, i__2, i__3; 00016 doublereal d__1, d__2, d__3, d__4; 00017 00018 /* Local variables */ 00019 static integer iend; 00020 extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal * 00021 , doublereal *, doublereal *, doublereal *); 00022 static doublereal norm; 00023 static integer i__, j, k, l, m, ii, en, jj, ll, nn; 00024 static doublereal si, ti, xi, yi, sr, tr, xr, yr; 00025 extern doublereal pythag_(doublereal *, doublereal *); 00026 extern /* Subroutine */ int csroot_(doublereal *, doublereal *, 00027 doublereal *, doublereal *); 00028 static integer ip1, lp1, itn, its; 00029 static doublereal zzi, zzr; 00030 static integer enm1; 00031 static doublereal tst1, tst2; 00032 00033 00034 00035 /* THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */ 00036 /* ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS */ 00037 /* AND WILKINSON. */ 00038 /* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */ 00039 /* THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */ 00040 /* (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */ 00041 00042 /* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */ 00043 /* OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR */ 00044 /* METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */ 00045 /* CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE */ 00046 /* THIS GENERAL MATRIX TO HESSENBERG FORM. */ 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 /* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */ 00057 /* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */ 00058 /* SET LOW=1, IGH=N. */ 00059 00060 /* ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */ 00061 /* FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED. */ 00062 /* ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS 00063 */ 00064 /* OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND */ 00065 /* ORTI(J) TO 0.0D0 FOR THESE ELEMENTS. */ 00066 00067 /* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00068 /* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */ 00069 /* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER */ 00070 /* INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE 00071 */ 00072 /* REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF */ 00073 /* THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE */ 00074 /* ARBITRARY. */ 00075 00076 /* ON OUTPUT */ 00077 00078 /* ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI */ 00079 /* HAVE BEEN DESTROYED. */ 00080 00081 /* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00082 /* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */ 00083 /* EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */ 00084 /* FOR INDICES IERR+1,...,N. */ 00085 00086 /* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */ 00087 /* RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS */ 00088 /* ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF */ 00089 /* THE EIGENVECTORS HAS BEEN FOUND. */ 00090 00091 /* IERR IS SET TO */ 00092 /* ZERO FOR NORMAL RETURN, */ 00093 /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */ 00094 /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */ 00095 00096 /* CALLS CDIV FOR COMPLEX DIVISION. */ 00097 /* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */ 00098 /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ 00099 00100 /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ 00101 /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 00102 */ 00103 00104 /* THIS VERSION DATED AUGUST 1983. */ 00105 00106 /* ------------------------------------------------------------------ 00107 */ 00108 00109 /* Parameter adjustments */ 00110 zi_dim1 = *nm; 00111 zi_offset = zi_dim1 + 1; 00112 zi -= zi_offset; 00113 zr_dim1 = *nm; 00114 zr_offset = zr_dim1 + 1; 00115 zr -= zr_offset; 00116 --wi; 00117 --wr; 00118 hi_dim1 = *nm; 00119 hi_offset = hi_dim1 + 1; 00120 hi -= hi_offset; 00121 hr_dim1 = *nm; 00122 hr_offset = hr_dim1 + 1; 00123 hr -= hr_offset; 00124 --orti; 00125 --ortr; 00126 00127 /* Function Body */ 00128 *ierr = 0; 00129 /* .......... INITIALIZE EIGENVECTOR MATRIX .......... */ 00130 i__1 = *n; 00131 for (j = 1; j <= i__1; ++j) { 00132 00133 i__2 = *n; 00134 for (i__ = 1; i__ <= i__2; ++i__) { 00135 zr[i__ + j * zr_dim1] = 0.; 00136 zi[i__ + j * zi_dim1] = 0.; 00137 /* L100: */ 00138 } 00139 zr[j + j * zr_dim1] = 1.; 00140 /* L101: */ 00141 } 00142 /* .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */ 00143 /* FROM THE INFORMATION LEFT BY CORTH .......... */ 00144 iend = *igh - *low - 1; 00145 if (iend < 0) { 00146 goto L180; 00147 } else if (iend == 0) { 00148 goto L150; 00149 } else { 00150 goto L105; 00151 } 00152 /* .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */ 00153 L105: 00154 i__1 = iend; 00155 for (ii = 1; ii <= i__1; ++ii) { 00156 i__ = *igh - ii; 00157 if (ortr[i__] == 0. && orti[i__] == 0.) { 00158 goto L140; 00159 } 00160 if (hr[i__ + (i__ - 1) * hr_dim1] == 0. && hi[i__ + (i__ - 1) * 00161 hi_dim1] == 0.) { 00162 goto L140; 00163 } 00164 /* .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ........ 00165 .. */ 00166 norm = hr[i__ + (i__ - 1) * hr_dim1] * ortr[i__] + hi[i__ + (i__ - 1) 00167 * hi_dim1] * orti[i__]; 00168 ip1 = i__ + 1; 00169 00170 i__2 = *igh; 00171 for (k = ip1; k <= i__2; ++k) { 00172 ortr[k] = hr[k + (i__ - 1) * hr_dim1]; 00173 orti[k] = hi[k + (i__ - 1) * hi_dim1]; 00174 /* L110: */ 00175 } 00176 00177 i__2 = *igh; 00178 for (j = i__; j <= i__2; ++j) { 00179 sr = 0.; 00180 si = 0.; 00181 00182 i__3 = *igh; 00183 for (k = i__; k <= i__3; ++k) { 00184 sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j * 00185 zi_dim1]; 00186 si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j * 00187 zr_dim1]; 00188 /* L115: */ 00189 } 00190 00191 sr /= norm; 00192 si /= norm; 00193 00194 i__3 = *igh; 00195 for (k = i__; k <= i__3; ++k) { 00196 zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si 00197 * orti[k]; 00198 zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si 00199 * ortr[k]; 00200 /* L120: */ 00201 } 00202 00203 /* L130: */ 00204 } 00205 00206 L140: 00207 ; 00208 } 00209 /* .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */ 00210 L150: 00211 l = *low + 1; 00212 00213 i__1 = *igh; 00214 for (i__ = l; i__ <= i__1; ++i__) { 00215 /* Computing MIN */ 00216 i__2 = i__ + 1; 00217 ll = min(i__2,*igh); 00218 if (hi[i__ + (i__ - 1) * hi_dim1] == 0.) { 00219 goto L170; 00220 } 00221 norm = pythag_(&hr[i__ + (i__ - 1) * hr_dim1], &hi[i__ + (i__ - 1) * 00222 hi_dim1]); 00223 yr = hr[i__ + (i__ - 1) * hr_dim1] / norm; 00224 yi = hi[i__ + (i__ - 1) * hi_dim1] / norm; 00225 hr[i__ + (i__ - 1) * hr_dim1] = norm; 00226 hi[i__ + (i__ - 1) * hi_dim1] = 0.; 00227 00228 i__2 = *n; 00229 for (j = i__; j <= i__2; ++j) { 00230 si = yr * hi[i__ + j * hi_dim1] - yi * hr[i__ + j * hr_dim1]; 00231 hr[i__ + j * hr_dim1] = yr * hr[i__ + j * hr_dim1] + yi * hi[i__ 00232 + j * hi_dim1]; 00233 hi[i__ + j * hi_dim1] = si; 00234 /* L155: */ 00235 } 00236 00237 i__2 = ll; 00238 for (j = 1; j <= i__2; ++j) { 00239 si = yr * hi[j + i__ * hi_dim1] + yi * hr[j + i__ * hr_dim1]; 00240 hr[j + i__ * hr_dim1] = yr * hr[j + i__ * hr_dim1] - yi * hi[j + 00241 i__ * hi_dim1]; 00242 hi[j + i__ * hi_dim1] = si; 00243 /* L160: */ 00244 } 00245 00246 i__2 = *igh; 00247 for (j = *low; j <= i__2; ++j) { 00248 si = yr * zi[j + i__ * zi_dim1] + yi * zr[j + i__ * zr_dim1]; 00249 zr[j + i__ * zr_dim1] = yr * zr[j + i__ * zr_dim1] - yi * zi[j + 00250 i__ * zi_dim1]; 00251 zi[j + i__ * zi_dim1] = si; 00252 /* L165: */ 00253 } 00254 00255 L170: 00256 ; 00257 } 00258 /* .......... STORE ROOTS ISOLATED BY CBAL .......... */ 00259 L180: 00260 i__1 = *n; 00261 for (i__ = 1; i__ <= i__1; ++i__) { 00262 if (i__ >= *low && i__ <= *igh) { 00263 goto L200; 00264 } 00265 wr[i__] = hr[i__ + i__ * hr_dim1]; 00266 wi[i__] = hi[i__ + i__ * hi_dim1]; 00267 L200: 00268 ; 00269 } 00270 00271 en = *igh; 00272 tr = 0.; 00273 ti = 0.; 00274 itn = *n * 30; 00275 /* .......... SEARCH FOR NEXT EIGENVALUE .......... */ 00276 L220: 00277 if (en < *low) { 00278 goto L680; 00279 } 00280 its = 0; 00281 enm1 = en - 1; 00282 /* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */ 00283 /* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */ 00284 L240: 00285 i__1 = en; 00286 for (ll = *low; ll <= i__1; ++ll) { 00287 l = en + *low - ll; 00288 if (l == *low) { 00289 goto L300; 00290 } 00291 tst1 = (d__1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d__1)) + (d__2 = hi[ 00292 l - 1 + (l - 1) * hi_dim1], abs(d__2)) + (d__3 = hr[l + l * 00293 hr_dim1], abs(d__3)) + (d__4 = hi[l + l * hi_dim1], abs(d__4)) 00294 ; 00295 tst2 = tst1 + (d__1 = hr[l + (l - 1) * hr_dim1], abs(d__1)); 00296 if (tst2 == tst1) { 00297 goto L300; 00298 } 00299 /* L260: */ 00300 } 00301 /* .......... FORM SHIFT .......... */ 00302 L300: 00303 if (l == en) { 00304 goto L660; 00305 } 00306 if (itn == 0) { 00307 goto L1000; 00308 } 00309 if (its == 10 || its == 20) { 00310 goto L320; 00311 } 00312 sr = hr[en + en * hr_dim1]; 00313 si = hi[en + en * hi_dim1]; 00314 xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1]; 00315 xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1]; 00316 if (xr == 0. && xi == 0.) { 00317 goto L340; 00318 } 00319 yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.; 00320 yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.; 00321 /* Computing 2nd power */ 00322 d__2 = yr; 00323 /* Computing 2nd power */ 00324 d__3 = yi; 00325 d__1 = d__2 * d__2 - d__3 * d__3 + xr; 00326 d__4 = yr * 2. * yi + xi; 00327 csroot_(&d__1, &d__4, &zzr, &zzi); 00328 if (yr * zzr + yi * zzi >= 0.) { 00329 goto L310; 00330 } 00331 zzr = -zzr; 00332 zzi = -zzi; 00333 L310: 00334 d__1 = yr + zzr; 00335 d__2 = yi + zzi; 00336 cdiv_(&xr, &xi, &d__1, &d__2, &xr, &xi); 00337 sr -= xr; 00338 si -= xi; 00339 goto L340; 00340 /* .......... FORM EXCEPTIONAL SHIFT .......... */ 00341 L320: 00342 sr = (d__1 = hr[en + enm1 * hr_dim1], abs(d__1)) + (d__2 = hr[enm1 + (en 00343 - 2) * hr_dim1], abs(d__2)); 00344 si = 0.; 00345 00346 L340: 00347 i__1 = en; 00348 for (i__ = *low; i__ <= i__1; ++i__) { 00349 hr[i__ + i__ * hr_dim1] -= sr; 00350 hi[i__ + i__ * hi_dim1] -= si; 00351 /* L360: */ 00352 } 00353 00354 tr += sr; 00355 ti += si; 00356 ++its; 00357 --itn; 00358 /* .......... REDUCE TO TRIANGLE (ROWS) .......... */ 00359 lp1 = l + 1; 00360 00361 i__1 = en; 00362 for (i__ = lp1; i__ <= i__1; ++i__) { 00363 sr = hr[i__ + (i__ - 1) * hr_dim1]; 00364 hr[i__ + (i__ - 1) * hr_dim1] = 0.; 00365 d__1 = pythag_(&hr[i__ - 1 + (i__ - 1) * hr_dim1], &hi[i__ - 1 + (i__ 00366 - 1) * hi_dim1]); 00367 norm = pythag_(&d__1, &sr); 00368 xr = hr[i__ - 1 + (i__ - 1) * hr_dim1] / norm; 00369 wr[i__ - 1] = xr; 00370 xi = hi[i__ - 1 + (i__ - 1) * hi_dim1] / norm; 00371 wi[i__ - 1] = xi; 00372 hr[i__ - 1 + (i__ - 1) * hr_dim1] = norm; 00373 hi[i__ - 1 + (i__ - 1) * hi_dim1] = 0.; 00374 hi[i__ + (i__ - 1) * hi_dim1] = sr / norm; 00375 00376 i__2 = *n; 00377 for (j = i__; j <= i__2; ++j) { 00378 yr = hr[i__ - 1 + j * hr_dim1]; 00379 yi = hi[i__ - 1 + j * hi_dim1]; 00380 zzr = hr[i__ + j * hr_dim1]; 00381 zzi = hi[i__ + j * hi_dim1]; 00382 hr[i__ - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i__ + (i__ - 1) 00383 * hi_dim1] * zzr; 00384 hi[i__ - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i__ + (i__ - 1) 00385 * hi_dim1] * zzi; 00386 hr[i__ + j * hr_dim1] = xr * zzr - xi * zzi - hi[i__ + (i__ - 1) * 00387 hi_dim1] * yr; 00388 hi[i__ + j * hi_dim1] = xr * zzi + xi * zzr - hi[i__ + (i__ - 1) * 00389 hi_dim1] * yi; 00390 /* L490: */ 00391 } 00392 00393 /* L500: */ 00394 } 00395 00396 si = hi[en + en * hi_dim1]; 00397 if (si == 0.) { 00398 goto L540; 00399 } 00400 norm = pythag_(&hr[en + en * hr_dim1], &si); 00401 sr = hr[en + en * hr_dim1] / norm; 00402 si /= norm; 00403 hr[en + en * hr_dim1] = norm; 00404 hi[en + en * hi_dim1] = 0.; 00405 if (en == *n) { 00406 goto L540; 00407 } 00408 ip1 = en + 1; 00409 00410 i__1 = *n; 00411 for (j = ip1; j <= i__1; ++j) { 00412 yr = hr[en + j * hr_dim1]; 00413 yi = hi[en + j * hi_dim1]; 00414 hr[en + j * hr_dim1] = sr * yr + si * yi; 00415 hi[en + j * hi_dim1] = sr * yi - si * yr; 00416 /* L520: */ 00417 } 00418 /* .......... INVERSE OPERATION (COLUMNS) .......... */ 00419 L540: 00420 i__1 = en; 00421 for (j = lp1; j <= i__1; ++j) { 00422 xr = wr[j - 1]; 00423 xi = wi[j - 1]; 00424 00425 i__2 = j; 00426 for (i__ = 1; i__ <= i__2; ++i__) { 00427 yr = hr[i__ + (j - 1) * hr_dim1]; 00428 yi = 0.; 00429 zzr = hr[i__ + j * hr_dim1]; 00430 zzi = hi[i__ + j * hi_dim1]; 00431 if (i__ == j) { 00432 goto L560; 00433 } 00434 yi = hi[i__ + (j - 1) * hi_dim1]; 00435 hi[i__ + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 00436 hi_dim1] * zzi; 00437 L560: 00438 hr[i__ + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 00439 hi_dim1] * zzr; 00440 hr[i__ + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 00441 hi_dim1] * yr; 00442 hi[i__ + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 00443 hi_dim1] * yi; 00444 /* L580: */ 00445 } 00446 00447 i__2 = *igh; 00448 for (i__ = *low; i__ <= i__2; ++i__) { 00449 yr = zr[i__ + (j - 1) * zr_dim1]; 00450 yi = zi[i__ + (j - 1) * zi_dim1]; 00451 zzr = zr[i__ + j * zr_dim1]; 00452 zzi = zi[i__ + j * zi_dim1]; 00453 zr[i__ + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) * 00454 hi_dim1] * zzr; 00455 zi[i__ + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) * 00456 hi_dim1] * zzi; 00457 zr[i__ + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) * 00458 hi_dim1] * yr; 00459 zi[i__ + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) * 00460 hi_dim1] * yi; 00461 /* L590: */ 00462 } 00463 00464 /* L600: */ 00465 } 00466 00467 if (si == 0.) { 00468 goto L240; 00469 } 00470 00471 i__1 = en; 00472 for (i__ = 1; i__ <= i__1; ++i__) { 00473 yr = hr[i__ + en * hr_dim1]; 00474 yi = hi[i__ + en * hi_dim1]; 00475 hr[i__ + en * hr_dim1] = sr * yr - si * yi; 00476 hi[i__ + en * hi_dim1] = sr * yi + si * yr; 00477 /* L630: */ 00478 } 00479 00480 i__1 = *igh; 00481 for (i__ = *low; i__ <= i__1; ++i__) { 00482 yr = zr[i__ + en * zr_dim1]; 00483 yi = zi[i__ + en * zi_dim1]; 00484 zr[i__ + en * zr_dim1] = sr * yr - si * yi; 00485 zi[i__ + en * zi_dim1] = sr * yi + si * yr; 00486 /* L640: */ 00487 } 00488 00489 goto L240; 00490 /* .......... A ROOT FOUND .......... */ 00491 L660: 00492 hr[en + en * hr_dim1] += tr; 00493 wr[en] = hr[en + en * hr_dim1]; 00494 hi[en + en * hi_dim1] += ti; 00495 wi[en] = hi[en + en * hi_dim1]; 00496 en = enm1; 00497 goto L220; 00498 /* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */ 00499 /* VECTORS OF UPPER TRIANGULAR FORM .......... */ 00500 L680: 00501 norm = 0.; 00502 00503 i__1 = *n; 00504 for (i__ = 1; i__ <= i__1; ++i__) { 00505 00506 i__2 = *n; 00507 for (j = i__; j <= i__2; ++j) { 00508 tr = (d__1 = hr[i__ + j * hr_dim1], abs(d__1)) + (d__2 = hi[i__ + 00509 j * hi_dim1], abs(d__2)); 00510 if (tr > norm) { 00511 norm = tr; 00512 } 00513 /* L720: */ 00514 } 00515 } 00516 00517 if (*n == 1 || norm == 0.) { 00518 goto L1001; 00519 } 00520 /* .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */ 00521 i__2 = *n; 00522 for (nn = 2; nn <= i__2; ++nn) { 00523 en = *n + 2 - nn; 00524 xr = wr[en]; 00525 xi = wi[en]; 00526 hr[en + en * hr_dim1] = 1.; 00527 hi[en + en * hi_dim1] = 0.; 00528 enm1 = en - 1; 00529 /* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */ 00530 i__1 = enm1; 00531 for (ii = 1; ii <= i__1; ++ii) { 00532 i__ = en - ii; 00533 zzr = 0.; 00534 zzi = 0.; 00535 ip1 = i__ + 1; 00536 00537 i__3 = en; 00538 for (j = ip1; j <= i__3; ++j) { 00539 zzr = zzr + hr[i__ + j * hr_dim1] * hr[j + en * hr_dim1] - hi[ 00540 i__ + j * hi_dim1] * hi[j + en * hi_dim1]; 00541 zzi = zzi + hr[i__ + j * hr_dim1] * hi[j + en * hi_dim1] + hi[ 00542 i__ + j * hi_dim1] * hr[j + en * hr_dim1]; 00543 /* L740: */ 00544 } 00545 00546 yr = xr - wr[i__]; 00547 yi = xi - wi[i__]; 00548 if (yr != 0. || yi != 0.) { 00549 goto L765; 00550 } 00551 tst1 = norm; 00552 yr = tst1; 00553 L760: 00554 yr *= .01; 00555 tst2 = norm + yr; 00556 if (tst2 > tst1) { 00557 goto L760; 00558 } 00559 L765: 00560 cdiv_(&zzr, &zzi, &yr, &yi, &hr[i__ + en * hr_dim1], &hi[i__ + en 00561 * hi_dim1]); 00562 /* .......... OVERFLOW CONTROL .......... */ 00563 tr = (d__1 = hr[i__ + en * hr_dim1], abs(d__1)) + (d__2 = hi[i__ 00564 + en * hi_dim1], abs(d__2)); 00565 if (tr == 0.) { 00566 goto L780; 00567 } 00568 tst1 = tr; 00569 tst2 = tst1 + 1. / tst1; 00570 if (tst2 > tst1) { 00571 goto L780; 00572 } 00573 i__3 = en; 00574 for (j = i__; j <= i__3; ++j) { 00575 hr[j + en * hr_dim1] /= tr; 00576 hi[j + en * hi_dim1] /= tr; 00577 /* L770: */ 00578 } 00579 00580 L780: 00581 ; 00582 } 00583 00584 /* L800: */ 00585 } 00586 /* .......... END BACKSUBSTITUTION .......... */ 00587 enm1 = *n - 1; 00588 /* .......... VECTORS OF ISOLATED ROOTS .......... */ 00589 i__2 = enm1; 00590 for (i__ = 1; i__ <= i__2; ++i__) { 00591 if (i__ >= *low && i__ <= *igh) { 00592 goto L840; 00593 } 00594 ip1 = i__ + 1; 00595 00596 i__1 = *n; 00597 for (j = ip1; j <= i__1; ++j) { 00598 zr[i__ + j * zr_dim1] = hr[i__ + j * hr_dim1]; 00599 zi[i__ + j * zi_dim1] = hi[i__ + j * hi_dim1]; 00600 /* L820: */ 00601 } 00602 00603 L840: 00604 ; 00605 } 00606 /* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */ 00607 /* VECTORS OF ORIGINAL FULL MATRIX. */ 00608 /* FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */ 00609 i__2 = enm1; 00610 for (jj = *low; jj <= i__2; ++jj) { 00611 j = *n + *low - jj; 00612 m = min(j,*igh); 00613 00614 i__1 = *igh; 00615 for (i__ = *low; i__ <= i__1; ++i__) { 00616 zzr = 0.; 00617 zzi = 0.; 00618 00619 i__3 = m; 00620 for (k = *low; k <= i__3; ++k) { 00621 zzr = zzr + zr[i__ + k * zr_dim1] * hr[k + j * hr_dim1] - zi[ 00622 i__ + k * zi_dim1] * hi[k + j * hi_dim1]; 00623 zzi = zzi + zr[i__ + k * zr_dim1] * hi[k + j * hi_dim1] + zi[ 00624 i__ + k * zi_dim1] * hr[k + j * hr_dim1]; 00625 /* L860: */ 00626 } 00627 00628 zr[i__ + j * zr_dim1] = zzr; 00629 zi[i__ + j * zi_dim1] = zzi; 00630 /* L880: */ 00631 } 00632 } 00633 00634 goto L1001; 00635 /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */ 00636 /* CONVERGED AFTER 30*N ITERATIONS .......... */ 00637 L1000: 00638 *ierr = en; 00639 L1001: 00640 return 0; 00641 } /* comqr2_ */ |