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_ */
|
