00012 {
00013     /* System generated locals */
00014     integer a_dim1, a_offset, z_dim1, z_offset, rm1_dim1, rm1_offset, i__1, 
00015             i__2, i__3;
00016     doublereal d__1, d__2;
00017 
00018     /* Builtin functions */
00019     double sqrt(doublereal);
00020 
00021     /* Local variables */
00022     extern /* Subroutine */ int cdiv_(doublereal *, doublereal *, doublereal *
00023             , doublereal *, doublereal *, doublereal *);
00024     static doublereal norm;
00025     static integer i__, j, k, l, s;
00026     static doublereal t, w, x, y;
00027     static integer n1;
00028     static doublereal normv;
00029     static integer ii;
00030     static doublereal ilambd;
00031     static integer ip, mp, ns, 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 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 REAL 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 /*        A CONTAINS THE HESSENBERG MATRIX. */
00059 
00060 /*        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
00061 /*          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE */
00062 /*          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  HQR, */
00063 /*          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. */
00064 
00065 /*        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE */
00066 /*          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS */
00067 /*          SPECIFIED BY SETTING SELECT(J) TO .TRUE.. */
00068 
00069 /*        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
00070 /*          COLUMNS REQUIRED TO STORE THE EIGENVECTORS TO BE FOUND. */
00071 /*          NOTE THAT TWO COLUMNS ARE REQUIRED TO STORE THE */
00072 /*          EIGENVECTOR CORRESPONDING TO A COMPLEX EIGENVALUE. */
00073 
00074 /*     ON OUTPUT */
00075 
00076 /*        A AND WI ARE UNALTERED. */
00077 
00078 /*        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED 
00079 */
00080 /*          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. */
00081 
00082 /*        SELECT MAY HAVE BEEN ALTERED.  IF THE ELEMENTS CORRESPONDING */
00083 /*          TO A PAIR OF CONJUGATE COMPLEX EIGENVALUES WERE EACH */
00084 /*          INITIALLY SET TO .TRUE., THE PROGRAM RESETS THE SECOND OF */
00085 /*          THE TWO ELEMENTS TO .FALSE.. */
00086 
00087 /*        M IS THE NUMBER OF COLUMNS ACTUALLY USED TO STORE */
00088 /*          THE EIGENVECTORS. */
00089 
00090 /*        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS. */
00091 /*          IF THE NEXT SELECTED EIGENVALUE IS REAL, THE NEXT COLUMN */
00092 /*          OF Z CONTAINS ITS EIGENVECTOR.  IF THE EIGENVALUE IS */
00093 /*          COMPLEX, THE NEXT TWO COLUMNS OF Z CONTAIN THE REAL AND */
00094 /*          IMAGINARY PARTS OF ITS EIGENVECTOR.  THE EIGENVECTORS ARE */
00095 /*          NORMALIZED SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. */
00096 /*          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. */
00097 
00098 /*        IERR IS SET TO */
00099 /*          ZERO       FOR NORMAL RETURN, */
00100 /*          -(2*N+1)   IF MORE THAN MM COLUMNS OF Z ARE NECESSARY */
00101 /*                     TO STORE THE EIGENVECTORS CORRESPONDING TO */
00102 /*                     THE SPECIFIED EIGENVALUES. */
00103 /*          -K         IF THE ITERATION CORRESPONDING TO THE K-TH */
00104 /*                     VALUE FAILS, */
00105 /*          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR. */
00106 
00107 /*        RM1, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RM1 
00108 */
00109 /*          IS SQUARE OF DIMENSION N BY N AND, AUGMENTED BY TWO COLUMNS */
00110 /*          OF Z, IS THE TRANSPOSE OF THE CORRESPONDING ALGOL B ARRAY. */
00111 
00112 /*     THE ALGOL PROCEDURE GUESSVEC APPEARS IN INVIT IN LINE. */
00113 
00114 /*     CALLS CDIV FOR COMPLEX DIVISION. */
00115 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00116 
00117 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00118 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00119 */
00120 
00121 /*     THIS VERSION DATED AUGUST 1983. */
00122 
00123 /*     ------------------------------------------------------------------ 
00124 */
00125 
00126     /* Parameter adjustments */
00127     --rv2;
00128     --rv1;
00129     rm1_dim1 = *n;
00130     rm1_offset = rm1_dim1 + 1;
00131     rm1 -= rm1_offset;
00132     --select;
00133     --wi;
00134     --wr;
00135     a_dim1 = *nm;
00136     a_offset = a_dim1 + 1;
00137     a -= a_offset;
00138     z_dim1 = *nm;
00139     z_offset = z_dim1 + 1;
00140     z__ -= z_offset;
00141 
00142     /* Function Body */
00143     *ierr = 0;
00144     uk = 0;
00145     s = 1;
00146 /*     .......... IP = 0, REAL EIGENVALUE */
00147 /*                     1, FIRST OF CONJUGATE COMPLEX PAIR */
00148 /*                    -1, SECOND OF CONJUGATE COMPLEX PAIR .......... */
00149     ip = 0;
00150     n1 = *n - 1;
00151 
00152     i__1 = *n;
00153     for (k = 1; k <= i__1; ++k) {
00154         if (wi[k] == 0. || ip < 0) {
00155             goto L100;
00156         }
00157         ip = 1;
00158         if (select[k] && select[k + 1]) {
00159             select[k + 1] = FALSE_;
00160         }
00161 L100:
00162         if (! select[k]) {
00163             goto L960;
00164         }
00165         if (wi[k] != 0.) {
00166             ++s;
00167         }
00168         if (s > *mm) {
00169             goto L1000;
00170         }
00171         if (uk >= k) {
00172             goto L200;
00173         }
00174 /*     .......... CHECK FOR POSSIBLE SPLITTING .......... */
00175         i__2 = *n;
00176         for (uk = k; uk <= i__2; ++uk) {
00177             if (uk == *n) {
00178                 goto L140;
00179             }
00180             if (a[uk + 1 + uk * a_dim1] == 0.) {
00181                 goto L140;
00182             }
00183 /* L120: */
00184         }
00185 /*     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK */
00186 /*                (HESSENBERG) MATRIX .......... */
00187 L140:
00188         norm = 0.;
00189         mp = 1;
00190 
00191         i__2 = uk;
00192         for (i__ = 1; i__ <= i__2; ++i__) {
00193             x = 0.;
00194 
00195             i__3 = uk;
00196             for (j = mp; j <= i__3; ++j) {
00197 /* L160: */
00198                 x += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00199             }
00200 
00201             if (x > norm) {
00202                 norm = x;
00203             }
00204             mp = i__;
00205 /* L180: */
00206         }
00207 /*     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION */
00208 /*                AND CLOSE ROOTS ARE MODIFIED BY EPS3 .......... */
00209         if (norm == 0.) {
00210             norm = 1.;
00211         }
00212         eps3 = epslon_(&norm);
00213 /*     .......... GROWTO IS THE CRITERION FOR THE GROWTH .......... */
00214         ukroot = (doublereal) uk;
00215         ukroot = sqrt(ukroot);
00216         growto = .1 / ukroot;
00217 L200:
00218         rlambd = wr[k];
00219         ilambd = wi[k];
00220         if (k == 1) {
00221             goto L280;
00222         }
00223         km1 = k - 1;
00224         goto L240;
00225 /*     .......... PERTURB EIGENVALUE IF IT IS CLOSE */
00226 /*                TO ANY PREVIOUS EIGENVALUE .......... */
00227 L220:
00228         rlambd += eps3;
00229 /*     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- .......... */
00230 L240:
00231         i__2 = km1;
00232         for (ii = 1; ii <= i__2; ++ii) {
00233             i__ = k - ii;
00234             if (select[i__] && (d__1 = wr[i__] - rlambd, abs(d__1)) < eps3 && 
00235                     (d__2 = wi[i__] - ilambd, abs(d__2)) < eps3) {
00236                 goto L220;
00237             }
00238 /* L260: */
00239         }
00240 
00241         wr[k] = rlambd;
00242 /*     .......... PERTURB CONJUGATE EIGENVALUE TO MATCH .......... */
00243         ip1 = k + ip;
00244         wr[ip1] = rlambd;
00245 /*     .......... FORM UPPER HESSENBERG A-RLAMBD*I (TRANSPOSED) */
00246 /*                AND INITIAL REAL VECTOR .......... */
00247 L280:
00248         mp = 1;
00249 
00250         i__2 = uk;
00251         for (i__ = 1; i__ <= i__2; ++i__) {
00252 
00253             i__3 = uk;
00254             for (j = mp; j <= i__3; ++j) {
00255 /* L300: */
00256                 rm1[j + i__ * rm1_dim1] = a[i__ + j * a_dim1];
00257             }
00258 
00259             rm1[i__ + i__ * rm1_dim1] -= rlambd;
00260             mp = i__;
00261             rv1[i__] = eps3;
00262 /* L320: */
00263         }
00264 
00265         its = 0;
00266         if (ilambd != 0.) {
00267             goto L520;
00268         }
00269 /*     .......... REAL EIGENVALUE. */
00270 /*                TRIANGULAR DECOMPOSITION WITH INTERCHANGES, */
00271 /*                REPLACING ZERO PIVOTS BY EPS3 .......... */
00272         if (uk == 1) {
00273             goto L420;
00274         }
00275 
00276         i__2 = uk;
00277         for (i__ = 2; i__ <= i__2; ++i__) {
00278             mp = i__ - 1;
00279             if ((d__1 = rm1[mp + i__ * rm1_dim1], abs(d__1)) <= (d__2 = rm1[
00280                     mp + mp * rm1_dim1], abs(d__2))) {
00281                 goto L360;
00282             }
00283 
00284             i__3 = uk;
00285             for (j = mp; j <= i__3; ++j) {
00286                 y = rm1[j + i__ * rm1_dim1];
00287                 rm1[j + i__ * rm1_dim1] = rm1[j + mp * rm1_dim1];
00288                 rm1[j + mp * rm1_dim1] = y;
00289 /* L340: */
00290             }
00291 
00292 L360:
00293             if (rm1[mp + mp * rm1_dim1] == 0.) {
00294                 rm1[mp + mp * rm1_dim1] = eps3;
00295             }
00296             x = rm1[mp + i__ * rm1_dim1] / rm1[mp + mp * rm1_dim1];
00297             if (x == 0.) {
00298                 goto L400;
00299             }
00300 
00301             i__3 = uk;
00302             for (j = i__; j <= i__3; ++j) {
00303 /* L380: */
00304                 rm1[j + i__ * rm1_dim1] -= x * rm1[j + mp * rm1_dim1];
00305             }
00306 
00307 L400:
00308             ;
00309         }
00310 
00311 L420:
00312         if (rm1[uk + uk * rm1_dim1] == 0.) {
00313             rm1[uk + uk * rm1_dim1] = eps3;
00314         }
00315 /*     .......... BACK SUBSTITUTION FOR REAL VECTOR */
00316 /*                FOR I=UK STEP -1 UNTIL 1 DO -- .......... */
00317 L440:
00318         i__2 = uk;
00319         for (ii = 1; ii <= i__2; ++ii) {
00320             i__ = uk + 1 - ii;
00321             y = rv1[i__];
00322             if (i__ == uk) {
00323                 goto L480;
00324             }
00325             ip1 = i__ + 1;
00326 
00327             i__3 = uk;
00328             for (j = ip1; j <= i__3; ++j) {
00329 /* L460: */
00330                 y -= rm1[j + i__ * rm1_dim1] * rv1[j];
00331             }
00332 
00333 L480:
00334             rv1[i__] = y / rm1[i__ + i__ * rm1_dim1];
00335 /* L500: */
00336         }
00337 
00338         goto L740;
00339 /*     .......... COMPLEX EIGENVALUE. */
00340 /*                TRIANGULAR DECOMPOSITION WITH INTERCHANGES, */
00341 /*                REPLACING ZERO PIVOTS BY EPS3.  STORE IMAGINARY */
00342 /*                PARTS IN UPPER TRIANGLE STARTING AT (1,3) ..........
00343  */
00344 L520:
00345         ns = *n - s;
00346         z__[(s - 1) * z_dim1 + 1] = -ilambd;
00347         z__[s * z_dim1 + 1] = 0.;
00348         if (*n == 2) {
00349             goto L550;
00350         }
00351         rm1[rm1_dim1 * 3 + 1] = -ilambd;
00352         z__[(s - 1) * z_dim1 + 1] = 0.;
00353         if (*n == 3) {
00354             goto L550;
00355         }
00356 
00357         i__2 = *n;
00358         for (i__ = 4; i__ <= i__2; ++i__) {
00359 /* L540: */
00360             rm1[i__ * rm1_dim1 + 1] = 0.;
00361         }
00362 
00363 L550:
00364         i__2 = uk;
00365         for (i__ = 2; i__ <= i__2; ++i__) {
00366             mp = i__ - 1;
00367             w = rm1[mp + i__ * rm1_dim1];
00368             if (i__ < *n) {
00369                 t = rm1[mp + (i__ + 1) * rm1_dim1];
00370             }
00371             if (i__ == *n) {
00372                 t = z__[mp + (s - 1) * z_dim1];
00373             }
00374             x = rm1[mp + mp * rm1_dim1] * rm1[mp + mp * rm1_dim1] + t * t;
00375             if (w * w <= x) {
00376                 goto L580;
00377             }
00378             x = rm1[mp + mp * rm1_dim1] / w;
00379             y = t / w;
00380             rm1[mp + mp * rm1_dim1] = w;
00381             if (i__ < *n) {
00382                 rm1[mp + (i__ + 1) * rm1_dim1] = 0.;
00383             }
00384             if (i__ == *n) {
00385                 z__[mp + (s - 1) * z_dim1] = 0.;
00386             }
00387 
00388             i__3 = uk;
00389             for (j = i__; j <= i__3; ++j) {
00390                 w = rm1[j + i__ * rm1_dim1];
00391                 rm1[j + i__ * rm1_dim1] = rm1[j + mp * rm1_dim1] - x * w;
00392                 rm1[j + mp * rm1_dim1] = w;
00393                 if (j < n1) {
00394                     goto L555;
00395                 }
00396                 l = j - ns;
00397                 z__[i__ + l * z_dim1] = z__[mp + l * z_dim1] - y * w;
00398                 z__[mp + l * z_dim1] = 0.;
00399                 goto L560;
00400 L555:
00401                 rm1[i__ + (j + 2) * rm1_dim1] = rm1[mp + (j + 2) * rm1_dim1] 
00402                         - y * w;
00403                 rm1[mp + (j + 2) * rm1_dim1] = 0.;
00404 L560:
00405                 ;
00406             }
00407 
00408             rm1[i__ + i__ * rm1_dim1] -= y * ilambd;
00409             if (i__ < n1) {
00410                 goto L570;
00411             }
00412             l = i__ - ns;
00413             z__[mp + l * z_dim1] = -ilambd;
00414             z__[i__ + l * z_dim1] += x * ilambd;
00415             goto L640;
00416 L570:
00417             rm1[mp + (i__ + 2) * rm1_dim1] = -ilambd;
00418             rm1[i__ + (i__ + 2) * rm1_dim1] += x * ilambd;
00419             goto L640;
00420 L580:
00421             if (x != 0.) {
00422                 goto L600;
00423             }
00424             rm1[mp + mp * rm1_dim1] = eps3;
00425             if (i__ < *n) {
00426                 rm1[mp + (i__ + 1) * rm1_dim1] = 0.;
00427             }
00428             if (i__ == *n) {
00429                 z__[mp + (s - 1) * z_dim1] = 0.;
00430             }
00431             t = 0.;
00432             x = eps3 * eps3;
00433 L600:
00434             w /= x;
00435             x = rm1[mp + mp * rm1_dim1] * w;
00436             y = -t * w;
00437 
00438             i__3 = uk;
00439             for (j = i__; j <= i__3; ++j) {
00440                 if (j < n1) {
00441                     goto L610;
00442                 }
00443                 l = j - ns;
00444                 t = z__[mp + l * z_dim1];
00445                 z__[i__ + l * z_dim1] = -x * t - y * rm1[j + mp * rm1_dim1];
00446                 goto L615;
00447 L610:
00448                 t = rm1[mp + (j + 2) * rm1_dim1];
00449                 rm1[i__ + (j + 2) * rm1_dim1] = -x * t - y * rm1[j + mp * 
00450                         rm1_dim1];
00451 L615:
00452                 rm1[j + i__ * rm1_dim1] = rm1[j + i__ * rm1_dim1] - x * rm1[j 
00453                         + mp * rm1_dim1] + y * t;
00454 /* L620: */
00455             }
00456 
00457             if (i__ < n1) {
00458                 goto L630;
00459             }
00460             l = i__ - ns;
00461             z__[i__ + l * z_dim1] -= ilambd;
00462             goto L640;
00463 L630:
00464             rm1[i__ + (i__ + 2) * rm1_dim1] -= ilambd;
00465 L640:
00466             ;
00467         }
00468 
00469         if (uk < n1) {
00470             goto L650;
00471         }
00472         l = uk - ns;
00473         t = z__[uk + l * z_dim1];
00474         goto L655;
00475 L650:
00476         t = rm1[uk + (uk + 2) * rm1_dim1];
00477 L655:
00478         if (rm1[uk + uk * rm1_dim1] == 0. && t == 0.) {
00479             rm1[uk + uk * rm1_dim1] = eps3;
00480         }
00481 /*     .......... BACK SUBSTITUTION FOR COMPLEX VECTOR */
00482 /*                FOR I=UK STEP -1 UNTIL 1 DO -- .......... */
00483 L660:
00484         i__2 = uk;
00485         for (ii = 1; ii <= i__2; ++ii) {
00486             i__ = uk + 1 - ii;
00487             x = rv1[i__];
00488             y = 0.;
00489             if (i__ == uk) {
00490                 goto L700;
00491             }
00492             ip1 = i__ + 1;
00493 
00494             i__3 = uk;
00495             for (j = ip1; j <= i__3; ++j) {
00496                 if (j < n1) {
00497                     goto L670;
00498                 }
00499                 l = j - ns;
00500                 t = z__[i__ + l * z_dim1];
00501                 goto L675;
00502 L670:
00503                 t = rm1[i__ + (j + 2) * rm1_dim1];
00504 L675:
00505                 x = x - rm1[j + i__ * rm1_dim1] * rv1[j] + t * rv2[j];
00506                 y = y - rm1[j + i__ * rm1_dim1] * rv2[j] - t * rv1[j];
00507 /* L680: */
00508             }
00509 
00510 L700:
00511             if (i__ < n1) {
00512                 goto L710;
00513             }
00514             l = i__ - ns;
00515             t = z__[i__ + l * z_dim1];
00516             goto L715;
00517 L710:
00518             t = rm1[i__ + (i__ + 2) * rm1_dim1];
00519 L715:
00520             cdiv_(&x, &y, &rm1[i__ + i__ * rm1_dim1], &t, &rv1[i__], &rv2[i__]
00521                     );
00522 /* L720: */
00523         }
00524 /*     .......... ACCEPTANCE TEST FOR REAL OR COMPLEX */
00525 /*                EIGENVECTOR AND NORMALIZATION .......... */
00526 L740:
00527         ++its;
00528         norm = 0.;
00529         normv = 0.;
00530 
00531         i__2 = uk;
00532         for (i__ = 1; i__ <= i__2; ++i__) {
00533             if (ilambd == 0.) {
00534                 x = (d__1 = rv1[i__], abs(d__1));
00535             }
00536             if (ilambd != 0.) {
00537                 x = pythag_(&rv1[i__], &rv2[i__]);
00538             }
00539             if (normv >= x) {
00540                 goto L760;
00541             }
00542             normv = x;
00543             j = i__;
00544 L760:
00545             norm += x;
00546 /* L780: */
00547         }
00548 
00549         if (norm < growto) {
00550             goto L840;
00551         }
00552 /*     .......... ACCEPT VECTOR .......... */
00553         x = rv1[j];
00554         if (ilambd == 0.) {
00555             x = 1. / x;
00556         }
00557         if (ilambd != 0.) {
00558             y = rv2[j];
00559         }
00560 
00561         i__2 = uk;
00562         for (i__ = 1; i__ <= i__2; ++i__) {
00563             if (ilambd != 0.) {
00564                 goto L800;
00565             }
00566             z__[i__ + s * z_dim1] = rv1[i__] * x;
00567             goto L820;
00568 L800:
00569             cdiv_(&rv1[i__], &rv2[i__], &x, &y, &z__[i__ + (s - 1) * z_dim1], 
00570                     &z__[i__ + s * z_dim1]);
00571 L820:
00572             ;
00573         }
00574 
00575         if (uk == *n) {
00576             goto L940;
00577         }
00578         j = uk + 1;
00579         goto L900;
00580 /*     .......... IN-LINE PROCEDURE FOR CHOOSING */
00581 /*                A NEW STARTING VECTOR .......... */
00582 L840:
00583         if (its >= uk) {
00584             goto L880;
00585         }
00586         x = ukroot;
00587         y = eps3 / (x + 1.);
00588         rv1[1] = eps3;
00589 
00590         i__2 = uk;
00591         for (i__ = 2; i__ <= i__2; ++i__) {
00592 /* L860: */
00593             rv1[i__] = y;
00594         }
00595 
00596         j = uk - its + 1;
00597         rv1[j] -= eps3 * x;
00598         if (ilambd == 0.) {
00599             goto L440;
00600         }
00601         goto L660;
00602 /*     .......... SET ERROR -- UNACCEPTED EIGENVECTOR .......... */
00603 L880:
00604         j = 1;
00605         *ierr = -k;
00606 /*     .......... SET REMAINING VECTOR COMPONENTS TO ZERO .......... 
00607 */
00608 L900:
00609         i__2 = *n;
00610         for (i__ = j; i__ <= i__2; ++i__) {
00611             z__[i__ + s * z_dim1] = 0.;
00612             if (ilambd != 0.) {
00613                 z__[i__ + (s - 1) * z_dim1] = 0.;
00614             }
00615 /* L920: */
00616         }
00617 
00618 L940:
00619         ++s;
00620 L960:
00621         if (ip == -1) {
00622             ip = 0;
00623         }
00624         if (ip == 1) {
00625             ip = -1;
00626         }
00627 /* L980: */
00628     }
00629 
00630     goto L1001;
00631 /*     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR */
00632 /*                SPACE REQUIRED .......... */
00633 L1000:
00634     if (*ierr != 0) {
00635         *ierr -= *n;
00636     }
00637     if (*ierr == 0) {
00638         *ierr = -((*n << 1) + 1);
00639     }
00640 L1001:
00641     *m = s - 1 - abs(ip);
00642     return 0;
00643 } /* invit_ */