00012 {
00013     /* System generated locals */
00014     integer z_dim1, z_offset, i__1, i__2, i__3;
00015     doublereal d__1, d__2, d__3, d__4;
00016 
00017     /* Builtin functions */
00018     double sqrt(doublereal);
00019 
00020     /* Local variables */
00021     static doublereal norm;
00022     static integer i__, j, p, q, r__, s;
00023     static doublereal u, v, order;
00024     static integer group;
00025     static doublereal x0, x1;
00026     static integer ii, jj, ip;
00027     static doublereal uk, xu;
00028     extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 
00029             *);
00030     static integer tag, its;
00031     static doublereal eps2, eps3, eps4;
00032 
00033 
00034 
00035 /*     THIS SUBROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH- */
00036 /*     NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. */
00037 /*     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */
00038 
00039 /*     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL */
00040 /*     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, */
00041 /*     USING INVERSE ITERATION. */
00042 
00043 /*     ON INPUT */
00044 
00045 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00046 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00047 /*          DIMENSION STATEMENT. */
00048 
00049 /*        N IS THE ORDER OF THE MATRIX. */
00050 
00051 /*        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
00052 
00053 /*        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
00054 /*          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY. */
00055 
00056 /*        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E, */
00057 /*          WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. */
00058 /*          E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN */
00059 /*          THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM */
00060 /*          OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST CONTAIN */
00061 /*          0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0D0 */
00062 /*          IF THE EIGENVALUES ARE IN DESCENDING ORDER.  IF  BISECT, */
00063 /*          TRIDIB, OR  IMTQLV  HAS BEEN USED TO FIND THE EIGENVALUES, */
00064 /*          THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE. */
00065 
00066 /*        M IS THE NUMBER OF SPECIFIED EIGENVALUES. */
00067 
00068 /*        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. 
00069 */
00070 
00071 /*        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
00072 /*          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
00073 /*          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
00074 /*          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC. 
00075 */
00076 
00077 /*     ON OUTPUT */
00078 
00079 /*        ALL INPUT ARRAYS ARE UNALTERED. */
00080 
00081 /*        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS. */
00082 /*          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. */
00083 
00084 /*        IERR IS SET TO */
00085 /*          ZERO       FOR NORMAL RETURN, */
00086 /*          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */
00087 /*                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS. */
00088 
00089 /*        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS. */
00090 
00091 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00092 
00093 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00094 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00095 */
00096 
00097 /*     THIS VERSION DATED AUGUST 1983. */
00098 
00099 /*     ------------------------------------------------------------------ 
00100 */
00101 
00102     /* Parameter adjustments */
00103     --rv6;
00104     --rv4;
00105     --rv3;
00106     --rv2;
00107     --rv1;
00108     --e2;
00109     --e;
00110     --d__;
00111     z_dim1 = *nm;
00112     z_offset = z_dim1 + 1;
00113     z__ -= z_offset;
00114     --ind;
00115     --w;
00116 
00117     /* Function Body */
00118     *ierr = 0;
00119     if (*m == 0) {
00120         goto L1001;
00121     }
00122     tag = 0;
00123     order = 1. - e2[1];
00124     q = 0;
00125 /*     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX .......... */
00126 L100:
00127     p = q + 1;
00128 
00129     i__1 = *n;
00130     for (q = p; q <= i__1; ++q) {
00131         if (q == *n) {
00132             goto L140;
00133         }
00134         if (e2[q + 1] == 0.) {
00135             goto L140;
00136         }
00137 /* L120: */
00138     }
00139 /*     .......... FIND VECTORS BY INVERSE ITERATION .......... */
00140 L140:
00141     ++tag;
00142     s = 0;
00143 
00144     i__1 = *m;
00145     for (r__ = 1; r__ <= i__1; ++r__) {
00146         if (ind[r__] != tag) {
00147             goto L920;
00148         }
00149         its = 1;
00150         x1 = w[r__];
00151         if (s != 0) {
00152             goto L510;
00153         }
00154 /*     .......... CHECK FOR ISOLATED ROOT .......... */
00155         xu = 1.;
00156         if (p != q) {
00157             goto L490;
00158         }
00159         rv6[p] = 1.;
00160         goto L870;
00161 L490:
00162         norm = (d__1 = d__[p], abs(d__1));
00163         ip = p + 1;
00164 
00165         i__2 = q;
00166         for (i__ = ip; i__ <= i__2; ++i__) {
00167 /* L500: */
00168 /* Computing MAX */
00169             d__3 = norm, d__4 = (d__1 = d__[i__], abs(d__1)) + (d__2 = e[i__],
00170                      abs(d__2));
00171             norm = max(d__3,d__4);
00172         }
00173 /*     .......... EPS2 IS THE CRITERION FOR GROUPING, */
00174 /*                EPS3 REPLACES ZERO PIVOTS AND EQUAL */
00175 /*                ROOTS ARE MODIFIED BY EPS3, */
00176 /*                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
00177 . */
00178         eps2 = norm * .001;
00179         eps3 = epslon_(&norm);
00180         uk = (doublereal) (q - p + 1);
00181         eps4 = uk * eps3;
00182         uk = eps4 / sqrt(uk);
00183         s = p;
00184 L505:
00185         group = 0;
00186         goto L520;
00187 /*     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */
00188 L510:
00189         if ((d__1 = x1 - x0, abs(d__1)) >= eps2) {
00190             goto L505;
00191         }
00192         ++group;
00193         if (order * (x1 - x0) <= 0.) {
00194             x1 = x0 + order * eps3;
00195         }
00196 /*     .......... ELIMINATION WITH INTERCHANGES AND */
00197 /*                INITIALIZATION OF VECTOR .......... */
00198 L520:
00199         v = 0.;
00200 
00201         i__2 = q;
00202         for (i__ = p; i__ <= i__2; ++i__) {
00203             rv6[i__] = uk;
00204             if (i__ == p) {
00205                 goto L560;
00206             }
00207             if ((d__1 = e[i__], abs(d__1)) < abs(u)) {
00208                 goto L540;
00209             }
00210 /*     .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF */
00211 /*                E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY ......
00212 .... */
00213             xu = u / e[i__];
00214             rv4[i__] = xu;
00215             rv1[i__ - 1] = e[i__];
00216             rv2[i__ - 1] = d__[i__] - x1;
00217             rv3[i__ - 1] = 0.;
00218             if (i__ != q) {
00219                 rv3[i__ - 1] = e[i__ + 1];
00220             }
00221             u = v - xu * rv2[i__ - 1];
00222             v = -xu * rv3[i__ - 1];
00223             goto L580;
00224 L540:
00225             xu = e[i__] / u;
00226             rv4[i__] = xu;
00227             rv1[i__ - 1] = u;
00228             rv2[i__ - 1] = v;
00229             rv3[i__ - 1] = 0.;
00230 L560:
00231             u = d__[i__] - x1 - xu * v;
00232             if (i__ != q) {
00233                 v = e[i__ + 1];
00234             }
00235 L580:
00236             ;
00237         }
00238 
00239         if (u == 0.) {
00240             u = eps3;
00241         }
00242         rv1[q] = u;
00243         rv2[q] = 0.;
00244         rv3[q] = 0.;
00245 /*     .......... BACK SUBSTITUTION */
00246 /*                FOR I=Q STEP -1 UNTIL P DO -- .......... */
00247 L600:
00248         i__2 = q;
00249         for (ii = p; ii <= i__2; ++ii) {
00250             i__ = p + q - ii;
00251             rv6[i__] = (rv6[i__] - u * rv2[i__] - v * rv3[i__]) / rv1[i__];
00252             v = u;
00253             u = rv6[i__];
00254 /* L620: */
00255         }
00256 /*     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */
00257 /*                MEMBERS OF GROUP .......... */
00258         if (group == 0) {
00259             goto L700;
00260         }
00261         j = r__;
00262 
00263         i__2 = group;
00264         for (jj = 1; jj <= i__2; ++jj) {
00265 L630:
00266             --j;
00267             if (ind[j] != tag) {
00268                 goto L630;
00269             }
00270             xu = 0.;
00271 
00272             i__3 = q;
00273             for (i__ = p; i__ <= i__3; ++i__) {
00274 /* L640: */
00275                 xu += rv6[i__] * z__[i__ + j * z_dim1];
00276             }
00277 
00278             i__3 = q;
00279             for (i__ = p; i__ <= i__3; ++i__) {
00280 /* L660: */
00281                 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00282             }
00283 
00284 /* L680: */
00285         }
00286 
00287 L700:
00288         norm = 0.;
00289 
00290         i__2 = q;
00291         for (i__ = p; i__ <= i__2; ++i__) {
00292 /* L720: */
00293             norm += (d__1 = rv6[i__], abs(d__1));
00294         }
00295 
00296         if (norm >= 1.) {
00297             goto L840;
00298         }
00299 /*     .......... FORWARD SUBSTITUTION .......... */
00300         if (its == 5) {
00301             goto L830;
00302         }
00303         if (norm != 0.) {
00304             goto L740;
00305         }
00306         rv6[s] = eps4;
00307         ++s;
00308         if (s > q) {
00309             s = p;
00310         }
00311         goto L780;
00312 L740:
00313         xu = eps4 / norm;
00314 
00315         i__2 = q;
00316         for (i__ = p; i__ <= i__2; ++i__) {
00317 /* L760: */
00318             rv6[i__] *= xu;
00319         }
00320 /*     .......... ELIMINATION OPERATIONS ON NEXT VECTOR */
00321 /*                ITERATE .......... */
00322 L780:
00323         i__2 = q;
00324         for (i__ = ip; i__ <= i__2; ++i__) {
00325             u = rv6[i__];
00326 /*     .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE */
00327 /*                WAS PERFORMED EARLIER IN THE */
00328 /*                TRIANGULARIZATION PROCESS .......... */
00329             if (rv1[i__ - 1] != e[i__]) {
00330                 goto L800;
00331             }
00332             u = rv6[i__ - 1];
00333             rv6[i__ - 1] = rv6[i__];
00334 L800:
00335             rv6[i__] = u - rv4[i__] * rv6[i__ - 1];
00336 /* L820: */
00337         }
00338 
00339         ++its;
00340         goto L600;
00341 /*     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */
00342 L830:
00343         *ierr = -r__;
00344         xu = 0.;
00345         goto L870;
00346 /*     .......... NORMALIZE SO THAT SUM OF SQUARES IS */
00347 /*                1 AND EXPAND TO FULL ORDER .......... */
00348 L840:
00349         u = 0.;
00350 
00351         i__2 = q;
00352         for (i__ = p; i__ <= i__2; ++i__) {
00353 /* L860: */
00354             u = pythag_(&u, &rv6[i__]);
00355         }
00356 
00357         xu = 1. / u;
00358 
00359 L870:
00360         i__2 = *n;
00361         for (i__ = 1; i__ <= i__2; ++i__) {
00362 /* L880: */
00363             z__[i__ + r__ * z_dim1] = 0.;
00364         }
00365 
00366         i__2 = q;
00367         for (i__ = p; i__ <= i__2; ++i__) {
00368 /* L900: */
00369             z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00370         }
00371 
00372         x0 = x1;
00373 L920:
00374         ;
00375     }
00376 
00377     if (q < *n) {
00378         goto L100;
00379     }
00380 L1001:
00381     return 0;
00382 } /* tinvit_ */