00013 {
00014     /* System generated locals */
00015     integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset, 
00016             zi_dim1, zi_offset, rm1_dim1, rm1_offset, rm2_dim1, rm2_offset, 
00017             i__1, i__2, i__3;
00018     doublereal d__1, d__2;
00019 
00020     /* Builtin functions */
00021     double sqrt(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, s;
00028     static doublereal x, y, normv;
00029     static integer ii;
00030     static doublereal ilambd;
00031     static integer mp, uk;
00032     static doublereal rlambd;
00033     extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 
00034             *);
00035     static integer km1, ip1;
00036     static doublereal growto, ukroot;
00037     static integer its;
00038     static doublereal eps3;
00039 
00040 
00041 
00042 /*     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT */
00043 /*     BY PETERS AND WILKINSON. */
00044 /*     HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971). */
00045 
00046 /*     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER */
00047 /*     HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, */
00048 /*     USING INVERSE ITERATION. */
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 /*        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
00059 /*          RESPECTIVELY, OF THE HESSENBERG MATRIX. */
00060 
00061 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
00062 /*          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE */
00063 /*          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  COMLR, */
00064 /*          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. */
00065 
00066 /*        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND.  THE */
00067 /*          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS */
00068 /*          SPECIFIED BY SETTING SELECT(J) TO .TRUE.. */
00069 
00070 /*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
00071 /*          EIGENVECTORS TO BE FOUND. */
00072 
00073 /*     ON OUTPUT */
00074 
00075 /*        AR, AI, WI, AND SELECT ARE UNALTERED. */
00076 
00077 /*        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED 
00078 */
00079 /*          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. */
00080 
00081 /*        M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND. */
00082 
00083 /*        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
00084 /*          OF THE EIGENVECTORS.  THE EIGENVECTORS ARE NORMALIZED */
00085 /*          SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. */
00086 /*          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. */
00087 
00088 /*        IERR IS SET TO */
00089 /*          ZERO       FOR NORMAL RETURN, */
00090 /*          -(2*N+1)   IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED, 
00091 */
00092 /*          -K         IF THE ITERATION CORRESPONDING TO THE K-TH */
00093 /*                     VALUE FAILS, */
00094 /*          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR. */
00095 
00096 /*        RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. */
00097 
00098 /*     THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE. */
00099 
00100 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00101 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00102 
00103 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00104 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00105 */
00106 
00107 /*     THIS VERSION DATED AUGUST 1983. */
00108 
00109 /*     ------------------------------------------------------------------ 
00110 */
00111 
00112     /* Parameter adjustments */
00113     --rv2;
00114     --rv1;
00115     rm2_dim1 = *n;
00116     rm2_offset = rm2_dim1 + 1;
00117     rm2 -= rm2_offset;
00118     rm1_dim1 = *n;
00119     rm1_offset = rm1_dim1 + 1;
00120     rm1 -= rm1_offset;
00121     --select;
00122     --wi;
00123     --wr;
00124     ai_dim1 = *nm;
00125     ai_offset = ai_dim1 + 1;
00126     ai -= ai_offset;
00127     ar_dim1 = *nm;
00128     ar_offset = ar_dim1 + 1;
00129     ar -= ar_offset;
00130     zi_dim1 = *nm;
00131     zi_offset = zi_dim1 + 1;
00132     zi -= zi_offset;
00133     zr_dim1 = *nm;
00134     zr_offset = zr_dim1 + 1;
00135     zr -= zr_offset;
00136 
00137     /* Function Body */
00138     *ierr = 0;
00139     uk = 0;
00140     s = 1;
00141 
00142     i__1 = *n;
00143     for (k = 1; k <= i__1; ++k) {
00144         if (! select[k]) {
00145             goto L980;
00146         }
00147         if (s > *mm) {
00148             goto L1000;
00149         }
00150         if (uk >= k) {
00151             goto L200;
00152         }
00153 /*     .......... CHECK FOR POSSIBLE SPLITTING .......... */
00154         i__2 = *n;
00155         for (uk = k; uk <= i__2; ++uk) {
00156             if (uk == *n) {
00157                 goto L140;
00158             }
00159             if (ar[uk + 1 + uk * ar_dim1] == 0. && ai[uk + 1 + uk * ai_dim1] 
00160                     == 0.) {
00161                 goto L140;
00162             }
00163 /* L120: */
00164         }
00165 /*     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK */
00166 /*                (HESSENBERG) MATRIX .......... */
00167 L140:
00168         norm = 0.;
00169         mp = 1;
00170 
00171         i__2 = uk;
00172         for (i__ = 1; i__ <= i__2; ++i__) {
00173             x = 0.;
00174 
00175             i__3 = uk;
00176             for (j = mp; j <= i__3; ++j) {
00177 /* L160: */
00178                 x += pythag_(&ar[i__ + j * ar_dim1], &ai[i__ + j * ai_dim1]);
00179             }
00180 
00181             if (x > norm) {
00182                 norm = x;
00183             }
00184             mp = i__;
00185 /* L180: */
00186         }
00187 /*     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION */
00188 /*                AND CLOSE ROOTS ARE MODIFIED BY EPS3 .......... */
00189         if (norm == 0.) {
00190             norm = 1.;
00191         }
00192         eps3 = epslon_(&norm);
00193 /*     .......... GROWTO IS THE CRITERION FOR GROWTH .......... */
00194         ukroot = (doublereal) uk;
00195         ukroot = sqrt(ukroot);
00196         growto = .1 / ukroot;
00197 L200:
00198         rlambd = wr[k];
00199         ilambd = wi[k];
00200         if (k == 1) {
00201             goto L280;
00202         }
00203         km1 = k - 1;
00204         goto L240;
00205 /*     .......... PERTURB EIGENVALUE IF IT IS CLOSE */
00206 /*                TO ANY PREVIOUS EIGENVALUE .......... */
00207 L220:
00208         rlambd += eps3;
00209 /*     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- .......... */
00210 L240:
00211         i__2 = km1;
00212         for (ii = 1; ii <= i__2; ++ii) {
00213             i__ = k - ii;
00214             if (select[i__] && (d__1 = wr[i__] - rlambd, abs(d__1)) < eps3 && 
00215                     (d__2 = wi[i__] - ilambd, abs(d__2)) < eps3) {
00216                 goto L220;
00217             }
00218 /* L260: */
00219         }
00220 
00221         wr[k] = rlambd;
00222 /*     .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I */
00223 /*                AND INITIAL COMPLEX VECTOR .......... */
00224 L280:
00225         mp = 1;
00226 
00227         i__2 = uk;
00228         for (i__ = 1; i__ <= i__2; ++i__) {
00229 
00230             i__3 = uk;
00231             for (j = mp; j <= i__3; ++j) {
00232                 rm1[i__ + j * rm1_dim1] = ar[i__ + j * ar_dim1];
00233                 rm2[i__ + j * rm2_dim1] = ai[i__ + j * ai_dim1];
00234 /* L300: */
00235             }
00236 
00237             rm1[i__ + i__ * rm1_dim1] -= rlambd;
00238             rm2[i__ + i__ * rm2_dim1] -= ilambd;
00239             mp = i__;
00240             rv1[i__] = eps3;
00241 /* L320: */
00242         }
00243 /*     .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES, */
00244 /*                REPLACING ZERO PIVOTS BY EPS3 .......... */
00245         if (uk == 1) {
00246             goto L420;
00247         }
00248 
00249         i__2 = uk;
00250         for (i__ = 2; i__ <= i__2; ++i__) {
00251             mp = i__ - 1;
00252             if (pythag_(&rm1[i__ + mp * rm1_dim1], &rm2[i__ + mp * rm2_dim1]) 
00253                     <= pythag_(&rm1[mp + mp * rm1_dim1], &rm2[mp + mp * 
00254                     rm2_dim1])) {
00255                 goto L360;
00256             }
00257 
00258             i__3 = uk;
00259             for (j = mp; j <= i__3; ++j) {
00260                 y = rm1[i__ + j * rm1_dim1];
00261                 rm1[i__ + j * rm1_dim1] = rm1[mp + j * rm1_dim1];
00262                 rm1[mp + j * rm1_dim1] = y;
00263                 y = rm2[i__ + j * rm2_dim1];
00264                 rm2[i__ + j * rm2_dim1] = rm2[mp + j * rm2_dim1];
00265                 rm2[mp + j * rm2_dim1] = y;
00266 /* L340: */
00267             }
00268 
00269 L360:
00270             if (rm1[mp + mp * rm1_dim1] == 0. && rm2[mp + mp * rm2_dim1] == 
00271                     0.) {
00272                 rm1[mp + mp * rm1_dim1] = eps3;
00273             }
00274             cdiv_(&rm1[i__ + mp * rm1_dim1], &rm2[i__ + mp * rm2_dim1], &rm1[
00275                     mp + mp * rm1_dim1], &rm2[mp + mp * rm2_dim1], &x, &y);
00276             if (x == 0. && y == 0.) {
00277                 goto L400;
00278             }
00279 
00280             i__3 = uk;
00281             for (j = i__; j <= i__3; ++j) {
00282                 rm1[i__ + j * rm1_dim1] = rm1[i__ + j * rm1_dim1] - x * rm1[
00283                         mp + j * rm1_dim1] + y * rm2[mp + j * rm2_dim1];
00284                 rm2[i__ + j * rm2_dim1] = rm2[i__ + j * rm2_dim1] - x * rm2[
00285                         mp + j * rm2_dim1] - y * rm1[mp + j * rm1_dim1];
00286 /* L380: */
00287             }
00288 
00289 L400:
00290             ;
00291         }
00292 
00293 L420:
00294         if (rm1[uk + uk * rm1_dim1] == 0. && rm2[uk + uk * rm2_dim1] == 0.) {
00295             rm1[uk + uk * rm1_dim1] = eps3;
00296         }
00297         its = 0;
00298 /*     .......... BACK SUBSTITUTION */
00299 /*                FOR I=UK STEP -1 UNTIL 1 DO -- .......... */
00300 L660:
00301         i__2 = uk;
00302         for (ii = 1; ii <= i__2; ++ii) {
00303             i__ = uk + 1 - ii;
00304             x = rv1[i__];
00305             y = 0.;
00306             if (i__ == uk) {
00307                 goto L700;
00308             }
00309             ip1 = i__ + 1;
00310 
00311             i__3 = uk;
00312             for (j = ip1; j <= i__3; ++j) {
00313                 x = x - rm1[i__ + j * rm1_dim1] * rv1[j] + rm2[i__ + j * 
00314                         rm2_dim1] * rv2[j];
00315                 y = y - rm1[i__ + j * rm1_dim1] * rv2[j] - rm2[i__ + j * 
00316                         rm2_dim1] * rv1[j];
00317 /* L680: */
00318             }
00319 
00320 L700:
00321             cdiv_(&x, &y, &rm1[i__ + i__ * rm1_dim1], &rm2[i__ + i__ * 
00322                     rm2_dim1], &rv1[i__], &rv2[i__]);
00323 /* L720: */
00324         }
00325 /*     .......... ACCEPTANCE TEST FOR EIGENVECTOR */
00326 /*                AND NORMALIZATION .......... */
00327         ++its;
00328         norm = 0.;
00329         normv = 0.;
00330 
00331         i__2 = uk;
00332         for (i__ = 1; i__ <= i__2; ++i__) {
00333             x = pythag_(&rv1[i__], &rv2[i__]);
00334             if (normv >= x) {
00335                 goto L760;
00336             }
00337             normv = x;
00338             j = i__;
00339 L760:
00340             norm += x;
00341 /* L780: */
00342         }
00343 
00344         if (norm < growto) {
00345             goto L840;
00346         }
00347 /*     .......... ACCEPT VECTOR .......... */
00348         x = rv1[j];
00349         y = rv2[j];
00350 
00351         i__2 = uk;
00352         for (i__ = 1; i__ <= i__2; ++i__) {
00353             cdiv_(&rv1[i__], &rv2[i__], &x, &y, &zr[i__ + s * zr_dim1], &zi[
00354                     i__ + s * zi_dim1]);
00355 /* L820: */
00356         }
00357 
00358         if (uk == *n) {
00359             goto L940;
00360         }
00361         j = uk + 1;
00362         goto L900;
00363 /*     .......... IN-LINE PROCEDURE FOR CHOOSING */
00364 /*                A NEW STARTING VECTOR .......... */
00365 L840:
00366         if (its >= uk) {
00367             goto L880;
00368         }
00369         x = ukroot;
00370         y = eps3 / (x + 1.);
00371         rv1[1] = eps3;
00372 
00373         i__2 = uk;
00374         for (i__ = 2; i__ <= i__2; ++i__) {
00375 /* L860: */
00376             rv1[i__] = y;
00377         }
00378 
00379         j = uk - its + 1;
00380         rv1[j] -= eps3 * x;
00381         goto L660;
00382 /*     .......... SET ERROR -- UNACCEPTED EIGENVECTOR .......... */
00383 L880:
00384         j = 1;
00385         *ierr = -k;
00386 /*     .......... SET REMAINING VECTOR COMPONENTS TO ZERO .......... 
00387 */
00388 L900:
00389         i__2 = *n;
00390         for (i__ = j; i__ <= i__2; ++i__) {
00391             zr[i__ + s * zr_dim1] = 0.;
00392             zi[i__ + s * zi_dim1] = 0.;
00393 /* L920: */
00394         }
00395 
00396 L940:
00397         ++s;
00398 L980:
00399         ;
00400     }
00401 
00402     goto L1001;
00403 /*     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR */
00404 /*                SPACE REQUIRED .......... */
00405 L1000:
00406     if (*ierr != 0) {
00407         *ierr -= *n;
00408     }
00409     if (*ierr == 0) {
00410         *ierr = -((*n << 1) + 1);
00411     }
00412 L1001:
00413     *m = s - 1;
00414     return 0;
00415 } /* cinvit_ */