00011 {
00012     /* System generated locals */
00013     integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5;
00014     doublereal d__1;
00015 
00016     /* Builtin functions */
00017     double sqrt(doublereal), d_sign(doublereal *, doublereal *);
00018 
00019     /* Local variables */
00020     static integer maxj, maxk;
00021     static doublereal norm;
00022     static integer i__, j, k, r__;
00023     static doublereal u, v, order;
00024     static integer group, m1;
00025     static doublereal x0, x1;
00026     static integer mb, m21, ii, ij, jj, kj;
00027     static doublereal uk, xu;
00028     extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal 
00029             *);
00030     static integer ij1, kj1, its;
00031     static doublereal eps2, eps3, eps4;
00032 
00033 
00034 
00035 /*     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC */
00036 /*     BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE 
00037 */
00038 /*     ITERATION.  THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS */
00039 /*     OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND */
00040 /*     COEFFICIENT MATRIX. */
00041 
00042 /*     ON INPUT */
00043 
00044 /*        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
00045 /*          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
00046 /*          DIMENSION STATEMENT. */
00047 
00048 /*        N IS THE ORDER OF THE MATRIX. */
00049 
00050 /*        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE */
00051 /*          BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) */
00052 /*          BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT 
00053 */
00054 /*          DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO */
00055 /*          SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE */
00056 /*          MATRIX.  IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS */
00057 /*          OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT */
00058 /*          SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT */
00059 /*          DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS */
00060 /*          CASE, MBW=2*MB-1. */
00061 
00062 /*        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
00063 /*          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL */
00064 /*          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
00065 /*          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
00066 /*          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
00067 /*          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. */
00068 /*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
00069 /*          EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS */
00070 /*          N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH */
00071 /*          ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF 
00072 */
00073 /*          COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2 */
00074 /*          POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY, */
00075 /*          AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB */
00076 /*          POSITIONS OF THE LAST COLUMN. */
00077 /*          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
00078 
00079 /*        E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS */
00080 /*            0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR */
00081 /*            2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER. */
00082 /*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
00083 /*          EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT */
00084 /*          MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT. */
00085 
00086 /*        M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF */
00087 /*          SYSTEMS OF LINEAR EQUATIONS. */
00088 
00089 /*        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. 
00090 */
00091 /*          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
00092 /*          EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY */
00093 /*          MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M. */
00094 
00095 /*        Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF */
00096 /*          THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS. 
00097 */
00098 
00099 /*        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
00100 /*          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
00101 
00102 /*     ON OUTPUT */
00103 
00104 /*        A AND W ARE UNALTERED. */
00105 
00106 /*        Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. */
00107 /*          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO.  IF THE */
00108 /*          SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, */
00109 /*          Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M). */
00110 
00111 /*        IERR IS SET TO */
00112 /*          ZERO       FOR NORMAL RETURN, */
00113 /*          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */
00114 /*                     EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH */
00115 /*                     SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR. */
00116 
00117 /*        RV AND RV6 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RV IS */
00118 /*          OF DIMENSION AT LEAST N*(2*MB-1).  IF THE SUBROUTINE */
00119 /*          IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE */
00120 /*          DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON */
00121 /*          RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV. */
00122 
00123 /*     CALLS PYTHAG FOR  DSQRT(A*A + B*B) . */
00124 
00125 /*     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
00126 /*     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY 
00127 */
00128 
00129 /*     THIS VERSION DATED AUGUST 1983. */
00130 
00131 /*     ------------------------------------------------------------------ 
00132 */
00133 
00134     /* Parameter adjustments */
00135     --rv6;
00136     a_dim1 = *nm;
00137     a_offset = a_dim1 + 1;
00138     a -= a_offset;
00139     z_dim1 = *nm;
00140     z_offset = z_dim1 + 1;
00141     z__ -= z_offset;
00142     --w;
00143     --rv;
00144 
00145     /* Function Body */
00146     *ierr = 0;
00147     if (*m == 0) {
00148         goto L1001;
00149     }
00150     mb = *mbw;
00151     if (*e21 < 0.) {
00152         mb = (*mbw + 1) / 2;
00153     }
00154     m1 = mb - 1;
00155     m21 = m1 + mb;
00156     order = 1. - abs(*e21);
00157 /*     .......... FIND VECTORS BY INVERSE ITERATION .......... */
00158     i__1 = *m;
00159     for (r__ = 1; r__ <= i__1; ++r__) {
00160         its = 1;
00161         x1 = w[r__];
00162         if (r__ != 1) {
00163             goto L100;
00164         }
00165 /*     .......... COMPUTE NORM OF MATRIX .......... */
00166         norm = 0.;
00167 
00168         i__2 = mb;
00169         for (j = 1; j <= i__2; ++j) {
00170             jj = mb + 1 - j;
00171             kj = jj + m1;
00172             ij = 1;
00173             v = 0.;
00174 
00175             i__3 = *n;
00176             for (i__ = jj; i__ <= i__3; ++i__) {
00177                 v += (d__1 = a[i__ + j * a_dim1], abs(d__1));
00178                 if (*e21 >= 0.) {
00179                     goto L40;
00180                 }
00181                 v += (d__1 = a[ij + kj * a_dim1], abs(d__1));
00182                 ++ij;
00183 L40:
00184                 ;
00185             }
00186 
00187             norm = max(norm,v);
00188 /* L60: */
00189         }
00190 
00191         if (*e21 < 0.) {
00192             norm *= .5;
00193         }
00194 /*     .......... EPS2 IS THE CRITERION FOR GROUPING, */
00195 /*                EPS3 REPLACES ZERO PIVOTS AND EQUAL */
00196 /*                ROOTS ARE MODIFIED BY EPS3, */
00197 /*                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
00198 . */
00199         if (norm == 0.) {
00200             norm = 1.;
00201         }
00202         eps2 = norm * .001 * abs(order);
00203         eps3 = epslon_(&norm);
00204         uk = (doublereal) (*n);
00205         uk = sqrt(uk);
00206         eps4 = uk * eps3;
00207 L80:
00208         group = 0;
00209         goto L120;
00210 /*     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */
00211 L100:
00212         if ((d__1 = x1 - x0, abs(d__1)) >= eps2) {
00213             goto L80;
00214         }
00215         ++group;
00216         if (order * (x1 - x0) <= 0.) {
00217             x1 = x0 + order * eps3;
00218         }
00219 /*     .......... EXPAND MATRIX, SUBTRACT EIGENVALUE, */
00220 /*                AND INITIALIZE VECTOR .......... */
00221 L120:
00222         i__2 = *n;
00223         for (i__ = 1; i__ <= i__2; ++i__) {
00224 /* Computing MIN */
00225             i__3 = 0, i__4 = i__ - m1;
00226             ij = i__ + min(i__3,i__4) * *n;
00227             kj = ij + mb * *n;
00228             ij1 = kj + m1 * *n;
00229             if (m1 == 0) {
00230                 goto L180;
00231             }
00232 
00233             i__3 = m1;
00234             for (j = 1; j <= i__3; ++j) {
00235                 if (ij > m1) {
00236                     goto L125;
00237                 }
00238                 if (ij > 0) {
00239                     goto L130;
00240                 }
00241                 rv[ij1] = 0.;
00242                 ij1 += *n;
00243                 goto L130;
00244 L125:
00245                 rv[ij] = a[i__ + j * a_dim1];
00246 L130:
00247                 ij += *n;
00248                 ii = i__ + j;
00249                 if (ii > *n) {
00250                     goto L150;
00251                 }
00252                 jj = mb - j;
00253                 if (*e21 >= 0.) {
00254                     goto L140;
00255                 }
00256                 ii = i__;
00257                 jj = mb + j;
00258 L140:
00259                 rv[kj] = a[ii + jj * a_dim1];
00260                 kj += *n;
00261 L150:
00262                 ;
00263             }
00264 
00265 L180:
00266             rv[ij] = a[i__ + mb * a_dim1] - x1;
00267             rv6[i__] = eps4;
00268             if (order == 0.) {
00269                 rv6[i__] = z__[i__ + r__ * z_dim1];
00270             }
00271 /* L200: */
00272         }
00273 
00274         if (m1 == 0) {
00275             goto L600;
00276         }
00277 /*     .......... ELIMINATION WITH INTERCHANGES .......... */
00278         i__2 = *n;
00279         for (i__ = 1; i__ <= i__2; ++i__) {
00280             ii = i__ + 1;
00281 /* Computing MIN */
00282             i__3 = i__ + m1 - 1;
00283             maxk = min(i__3,*n);
00284 /* Computing MIN */
00285             i__3 = *n - i__, i__4 = m21 - 2;
00286             maxj = min(i__3,i__4) * *n;
00287 
00288             i__3 = maxk;
00289             for (k = i__; k <= i__3; ++k) {
00290                 kj1 = k;
00291                 j = kj1 + *n;
00292                 jj = j + maxj;
00293 
00294                 i__4 = jj;
00295                 i__5 = *n;
00296                 for (kj = j; i__5 < 0 ? kj >= i__4 : kj <= i__4; kj += i__5) {
00297                     rv[kj1] = rv[kj];
00298                     kj1 = kj;
00299 /* L340: */
00300                 }
00301 
00302                 rv[kj1] = 0.;
00303 /* L360: */
00304             }
00305 
00306             if (i__ == *n) {
00307                 goto L580;
00308             }
00309             u = 0.;
00310 /* Computing MIN */
00311             i__3 = i__ + m1;
00312             maxk = min(i__3,*n);
00313 /* Computing MIN */
00314             i__3 = *n - ii, i__5 = m21 - 2;
00315             maxj = min(i__3,i__5) * *n;
00316 
00317             i__3 = maxk;
00318             for (j = i__; j <= i__3; ++j) {
00319                 if ((d__1 = rv[j], abs(d__1)) < abs(u)) {
00320                     goto L450;
00321                 }
00322                 u = rv[j];
00323                 k = j;
00324 L450:
00325                 ;
00326             }
00327 
00328             j = i__ + *n;
00329             jj = j + maxj;
00330             if (k == i__) {
00331                 goto L520;
00332             }
00333             kj = k;
00334 
00335             i__3 = jj;
00336             i__5 = *n;
00337             for (ij = i__; i__5 < 0 ? ij >= i__3 : ij <= i__3; ij += i__5) {
00338                 v = rv[ij];
00339                 rv[ij] = rv[kj];
00340                 rv[kj] = v;
00341                 kj += *n;
00342 /* L500: */
00343             }
00344 
00345             if (order != 0.) {
00346                 goto L520;
00347             }
00348             v = rv6[i__];
00349             rv6[i__] = rv6[k];
00350             rv6[k] = v;
00351 L520:
00352             if (u == 0.) {
00353                 goto L580;
00354             }
00355 
00356             i__5 = maxk;
00357             for (k = ii; k <= i__5; ++k) {
00358                 v = rv[k] / u;
00359                 kj = k;
00360 
00361                 i__3 = jj;
00362                 i__4 = *n;
00363                 for (ij = j; i__4 < 0 ? ij >= i__3 : ij <= i__3; ij += i__4) {
00364                     kj += *n;
00365                     rv[kj] -= v * rv[ij];
00366 /* L540: */
00367                 }
00368 
00369                 if (order == 0.) {
00370                     rv6[k] -= v * rv6[i__];
00371                 }
00372 /* L560: */
00373             }
00374 
00375 L580:
00376             ;
00377         }
00378 /*     .......... BACK SUBSTITUTION */
00379 /*                FOR I=N STEP -1 UNTIL 1 DO -- .......... */
00380 L600:
00381         i__2 = *n;
00382         for (ii = 1; ii <= i__2; ++ii) {
00383             i__ = *n + 1 - ii;
00384             maxj = min(ii,m21);
00385             if (maxj == 1) {
00386                 goto L620;
00387             }
00388             ij1 = i__;
00389             j = ij1 + *n;
00390             jj = j + (maxj - 2) * *n;
00391 
00392             i__5 = jj;
00393             i__4 = *n;
00394             for (ij = j; i__4 < 0 ? ij >= i__5 : ij <= i__5; ij += i__4) {
00395                 ++ij1;
00396                 rv6[i__] -= rv[ij] * rv6[ij1];
00397 /* L610: */
00398             }
00399 
00400 L620:
00401             v = rv[i__];
00402             if (abs(v) >= eps3) {
00403                 goto L625;
00404             }
00405 /*     .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM .....
00406 ..... */
00407             if (order == 0.) {
00408                 *ierr = -r__;
00409             }
00410             v = d_sign(&eps3, &v);
00411 L625:
00412             rv6[i__] /= v;
00413 /* L630: */
00414         }
00415 
00416         xu = 1.;
00417         if (order == 0.) {
00418             goto L870;
00419         }
00420 /*     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */
00421 /*                MEMBERS OF GROUP .......... */
00422         if (group == 0) {
00423             goto L700;
00424         }
00425 
00426         i__2 = group;
00427         for (jj = 1; jj <= i__2; ++jj) {
00428             j = r__ - group - 1 + jj;
00429             xu = 0.;
00430 
00431             i__4 = *n;
00432             for (i__ = 1; i__ <= i__4; ++i__) {
00433 /* L640: */
00434                 xu += rv6[i__] * z__[i__ + j * z_dim1];
00435             }
00436 
00437             i__4 = *n;
00438             for (i__ = 1; i__ <= i__4; ++i__) {
00439 /* L660: */
00440                 rv6[i__] -= xu * z__[i__ + j * z_dim1];
00441             }
00442 
00443 /* L680: */
00444         }
00445 
00446 L700:
00447         norm = 0.;
00448 
00449         i__2 = *n;
00450         for (i__ = 1; i__ <= i__2; ++i__) {
00451 /* L720: */
00452             norm += (d__1 = rv6[i__], abs(d__1));
00453         }
00454 
00455         if (norm >= .1) {
00456             goto L840;
00457         }
00458 /*     .......... IN-LINE PROCEDURE FOR CHOOSING */
00459 /*                A NEW STARTING VECTOR .......... */
00460         if (its >= *n) {
00461             goto L830;
00462         }
00463         ++its;
00464         xu = eps4 / (uk + 1.);
00465         rv6[1] = eps4;
00466 
00467         i__2 = *n;
00468         for (i__ = 2; i__ <= i__2; ++i__) {
00469 /* L760: */
00470             rv6[i__] = xu;
00471         }
00472 
00473         rv6[its] -= eps4 * uk;
00474         goto L600;
00475 /*     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */
00476 L830:
00477         *ierr = -r__;
00478         xu = 0.;
00479         goto L870;
00480 /*     .......... NORMALIZE SO THAT SUM OF SQUARES IS */
00481 /*                1 AND EXPAND TO FULL ORDER .......... */
00482 L840:
00483         u = 0.;
00484 
00485         i__2 = *n;
00486         for (i__ = 1; i__ <= i__2; ++i__) {
00487 /* L860: */
00488             u = pythag_(&u, &rv6[i__]);
00489         }
00490 
00491         xu = 1. / u;
00492 
00493 L870:
00494         i__2 = *n;
00495         for (i__ = 1; i__ <= i__2; ++i__) {
00496 /* L900: */
00497             z__[i__ + r__ * z_dim1] = rv6[i__] * xu;
00498         }
00499 
00500         x0 = x1;
00501 /* L920: */
00502     }
00503 
00504 L1001:
00505     return 0;
00506 } /* bandv_ */