/* bandv.f -- translated by f2c (version 19961017). You must link the resulting object file with the libraries: -lf2c -lm (in that order) */ #include #include "f2c.h" /* Subroutine */ int bandv_(integer *nm, integer *n, integer *mbw, doublereal *a, doublereal *e21, integer *m, doublereal *w, doublereal *z__, integer *ierr, integer *nv, doublereal *rv, doublereal *rv6) { /* System generated locals */ integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5; doublereal d__1; /* Builtin functions */ double sqrt(doublereal), d_sign(doublereal *, doublereal *); /* Local variables */ integer maxj, maxk; doublereal norm; integer i__, j, k, r__; doublereal u, v, order; integer group=0, m1; doublereal x0=0.0, x1; integer mb, m21, ii, ij, jj, kj; doublereal uk=0.0, xu; extern doublereal pythag_(doublereal *, doublereal *), epslon_(doublereal *); integer ij1, kj1, its; doublereal eps2=0.0, eps3, eps4=0.0; /* THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC */ /* BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE */ /* ITERATION. THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS */ /* OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND */ /* COEFFICIENT MATRIX. */ /* ON INPUT */ /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ /* DIMENSION STATEMENT. */ /* N IS THE ORDER OF THE MATRIX. */ /* MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE */ /* BAND MATRIX. IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) */ /* BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT */ /* DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO */ /* SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE */ /* MATRIX. IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS */ /* OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT */ /* SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT */ /* DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS */ /* CASE, MBW=2*MB-1. */ /* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */ /* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */ /* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */ /* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */ /* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */ /* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. */ /* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */ /* EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS */ /* N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH */ /* ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF */ /* COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2 */ /* POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY, */ /* AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB */ /* POSITIONS OF THE LAST COLUMN. */ /* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */ /* E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS */ /* 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR */ /* 2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER. */ /* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */ /* EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT */ /* MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT. */ /* M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF */ /* SYSTEMS OF LINEAR EQUATIONS. */ /* W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER. */ /* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */ /* EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY */ /* MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M. */ /* Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF */ /* THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS. */ /* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */ /* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */ /* ON OUTPUT */ /* A AND W ARE UNALTERED. */ /* Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. */ /* ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. IF THE */ /* SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, */ /* Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M). */ /* IERR IS SET TO */ /* ZERO FOR NORMAL RETURN, */ /* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */ /* EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH */ /* SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR. */ /* RV AND RV6 ARE TEMPORARY STORAGE ARRAYS. NOTE THAT RV IS */ /* OF DIMENSION AT LEAST N*(2*MB-1). IF THE SUBROUTINE */ /* IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE */ /* DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON */ /* RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV. */ /* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY */ /* THIS VERSION DATED AUGUST 1983. */ /* ------------------------------------------------------------------ */ /* Parameter adjustments */ --rv6; a_dim1 = *nm; a_offset = a_dim1 + 1; a -= a_offset; z_dim1 = *nm; z_offset = z_dim1 + 1; z__ -= z_offset; --w; --rv; /* Function Body */ *ierr = 0; if (*m == 0) { goto L1001; } mb = *mbw; if (*e21 < 0.) { mb = (*mbw + 1) / 2; } m1 = mb - 1; m21 = m1 + mb; order = 1. - abs(*e21); /* .......... FIND VECTORS BY INVERSE ITERATION .......... */ i__1 = *m; for (r__ = 1; r__ <= i__1; ++r__) { its = 1; x1 = w[r__]; if (r__ != 1) { goto L100; } /* .......... COMPUTE NORM OF MATRIX .......... */ norm = 0.; i__2 = mb; for (j = 1; j <= i__2; ++j) { jj = mb + 1 - j; kj = jj + m1; ij = 1; v = 0.; i__3 = *n; for (i__ = jj; i__ <= i__3; ++i__) { v += (d__1 = a[i__ + j * a_dim1], abs(d__1)); if (*e21 >= 0.) { goto L40; } v += (d__1 = a[ij + kj * a_dim1], abs(d__1)); ++ij; L40: ; } norm = max(norm,v); /* L60: */ } if (*e21 < 0.) { norm *= .5; } /* .......... EPS2 IS THE CRITERION FOR GROUPING, */ /* EPS3 REPLACES ZERO PIVOTS AND EQUAL */ /* ROOTS ARE MODIFIED BY EPS3, */ /* EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ......... . */ if (norm == 0.) { norm = 1.; } eps2 = norm * .001 * abs(order); eps3 = epslon_(&norm); uk = (doublereal) (*n); uk = sqrt(uk); eps4 = uk * eps3; L80: group = 0; goto L120; /* .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */ L100: if ((d__1 = x1 - x0, abs(d__1)) >= eps2) { goto L80; } ++group; if (order * (x1 - x0) <= 0.) { x1 = x0 + order * eps3; } /* .......... EXPAND MATRIX, SUBTRACT EIGENVALUE, */ /* AND INITIALIZE VECTOR .......... */ L120: i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { /* Computing MIN */ i__3 = 0, i__4 = i__ - m1; ij = i__ + min(i__3,i__4) * *n; kj = ij + mb * *n; ij1 = kj + m1 * *n; if (m1 == 0) { goto L180; } i__3 = m1; for (j = 1; j <= i__3; ++j) { if (ij > m1) { goto L125; } if (ij > 0) { goto L130; } rv[ij1] = 0.; ij1 += *n; goto L130; L125: rv[ij] = a[i__ + j * a_dim1]; L130: ij += *n; ii = i__ + j; if (ii > *n) { goto L150; } jj = mb - j; if (*e21 >= 0.) { goto L140; } ii = i__; jj = mb + j; L140: rv[kj] = a[ii + jj * a_dim1]; kj += *n; L150: ; } L180: rv[ij] = a[i__ + mb * a_dim1] - x1; rv6[i__] = eps4; if (order == 0.) { rv6[i__] = z__[i__ + r__ * z_dim1]; } /* L200: */ } if (m1 == 0) { goto L600; } /* .......... ELIMINATION WITH INTERCHANGES .......... */ i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { ii = i__ + 1; /* Computing MIN */ i__3 = i__ + m1 - 1; maxk = min(i__3,*n); /* Computing MIN */ i__3 = *n - i__, i__4 = m21 - 2; maxj = min(i__3,i__4) * *n; i__3 = maxk; for (k = i__; k <= i__3; ++k) { kj1 = k; j = kj1 + *n; jj = j + maxj; i__4 = jj; i__5 = *n; for (kj = j; i__5 < 0 ? kj >= i__4 : kj <= i__4; kj += i__5) { rv[kj1] = rv[kj]; kj1 = kj; /* L340: */ } rv[kj1] = 0.; /* L360: */ } if (i__ == *n) { goto L580; } u = 0.; /* Computing MIN */ i__3 = i__ + m1; maxk = min(i__3,*n); /* Computing MIN */ i__3 = *n - ii, i__5 = m21 - 2; maxj = min(i__3,i__5) * *n; i__3 = maxk; for (j = i__; j <= i__3; ++j) { if ((d__1 = rv[j], abs(d__1)) < abs(u)) { goto L450; } u = rv[j]; k = j; L450: ; } j = i__ + *n; jj = j + maxj; if (k == i__) { goto L520; } kj = k; i__3 = jj; i__5 = *n; for (ij = i__; i__5 < 0 ? ij >= i__3 : ij <= i__3; ij += i__5) { v = rv[ij]; rv[ij] = rv[kj]; rv[kj] = v; kj += *n; /* L500: */ } if (order != 0.) { goto L520; } v = rv6[i__]; rv6[i__] = rv6[k]; rv6[k] = v; L520: if (u == 0.) { goto L580; } i__5 = maxk; for (k = ii; k <= i__5; ++k) { v = rv[k] / u; kj = k; i__3 = jj; i__4 = *n; for (ij = j; i__4 < 0 ? ij >= i__3 : ij <= i__3; ij += i__4) { kj += *n; rv[kj] -= v * rv[ij]; /* L540: */ } if (order == 0.) { rv6[k] -= v * rv6[i__]; } /* L560: */ } L580: ; } /* .......... BACK SUBSTITUTION */ /* FOR I=N STEP -1 UNTIL 1 DO -- .......... */ L600: i__2 = *n; for (ii = 1; ii <= i__2; ++ii) { i__ = *n + 1 - ii; maxj = min(ii,m21); if (maxj == 1) { goto L620; } ij1 = i__; j = ij1 + *n; jj = j + (maxj - 2) * *n; i__5 = jj; i__4 = *n; for (ij = j; i__4 < 0 ? ij >= i__5 : ij <= i__5; ij += i__4) { ++ij1; rv6[i__] -= rv[ij] * rv6[ij1]; /* L610: */ } L620: v = rv[i__]; if (abs(v) >= eps3) { goto L625; } /* .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM ..... ..... */ if (order == 0.) { *ierr = -r__; } v = d_sign(&eps3, &v); L625: rv6[i__] /= v; /* L630: */ } xu = 1.; if (order == 0.) { goto L870; } /* .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */ /* MEMBERS OF GROUP .......... */ if (group == 0) { goto L700; } i__2 = group; for (jj = 1; jj <= i__2; ++jj) { j = r__ - group - 1 + jj; xu = 0.; i__4 = *n; for (i__ = 1; i__ <= i__4; ++i__) { /* L640: */ xu += rv6[i__] * z__[i__ + j * z_dim1]; } i__4 = *n; for (i__ = 1; i__ <= i__4; ++i__) { /* L660: */ rv6[i__] -= xu * z__[i__ + j * z_dim1]; } /* L680: */ } L700: norm = 0.; i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { /* L720: */ norm += (d__1 = rv6[i__], abs(d__1)); } if (norm >= .1) { goto L840; } /* .......... IN-LINE PROCEDURE FOR CHOOSING */ /* A NEW STARTING VECTOR .......... */ if (its >= *n) { goto L830; } ++its; xu = eps4 / (uk + 1.); rv6[1] = eps4; i__2 = *n; for (i__ = 2; i__ <= i__2; ++i__) { /* L760: */ rv6[i__] = xu; } rv6[its] -= eps4 * uk; goto L600; /* .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */ L830: *ierr = -r__; xu = 0.; goto L870; /* .......... NORMALIZE SO THAT SUM OF SQUARES IS */ /* 1 AND EXPAND TO FULL ORDER .......... */ L840: u = 0.; i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { /* L860: */ u = pythag_(&u, &rv6[i__]); } xu = 1. / u; L870: i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { /* L900: */ z__[i__ + r__ * z_dim1] = rv6[i__] * xu; } x0 = x1; /* L920: */ } L1001: return 0; } /* bandv_ */