*DECK TINVIT SUBROUTINE TINVIT (NM, N, D, E, E2, M, W, IND, Z, IERR, RV1, RV2, + RV3, RV4, RV6) C***BEGIN PROLOGUE TINVIT C***PURPOSE Compute the eigenvectors of symmetric tridiagonal matrix C corresponding to specified eigenvalues, using inverse C iteration. C***LIBRARY SLATEC (EISPACK) C***CATEGORY D4C3 C***TYPE SINGLE PRECISION (TINVIT-S) C***KEYWORDS EIGENVECTORS, EISPACK C***AUTHOR Smith, B. T., et al. C***DESCRIPTION C C This subroutine is a translation of the inverse iteration tech- C nique in the ALGOL procedure TRISTURM by Peters and Wilkinson. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). C C This subroutine finds those eigenvectors of a TRIDIAGONAL C SYMMETRIC matrix corresponding to specified eigenvalues, C using inverse iteration. C C On Input C C NM must be set to the row dimension of the two-dimensional C array parameter, Z, as declared in the calling program C dimension statement. NM is an INTEGER variable. C C N is the order of the matrix. N is an INTEGER variable. C N must be less than or equal to NM. C C D contains the diagonal elements of the symmetric tridiagonal C matrix. D is a one-dimensional REAL array, dimensioned D(N). C C E contains the subdiagonal elements of the symmetric C tridiagonal matrix in its last N-1 positions. E(1) is C arbitrary. E is a one-dimensional REAL array, dimensioned C E(N). C C E2 contains the squares of the corresponding elements of E, C with zeros corresponding to negligible elements of E. C E(I) is considered negligible if it is not larger than C the product of the relative machine precision and the sum C of the magnitudes of D(I) and D(I-1). E2(1) must contain C 0.0e0 if the eigenvalues are in ascending order, or 2.0e0 C if the eigenvalues are in descending order. If BISECT, C TRIDIB, or IMTQLV has been used to find the eigenvalues, C their output E2 array is exactly what is expected here. C E2 is a one-dimensional REAL array, dimensioned E2(N). C C M is the number of specified eigenvalues for which eigenvectors C are to be determined. M is an INTEGER variable. C C W contains the M eigenvalues in ascending or descending order. C W is a one-dimensional REAL array, dimensioned W(M). C C IND contains in its first M positions the submatrix indices C associated with the corresponding eigenvalues in W -- C 1 for eigenvalues belonging to the first submatrix from C the top, 2 for those belonging to the second submatrix, etc. C If BISECT or TRIDIB has been used to determine the C eigenvalues, their output IND array is suitable for input C to TINVIT. IND is a one-dimensional INTEGER array, C dimensioned IND(M). C C On Output C C ** All input arrays are unaltered.** C C Z contains the associated set of orthonormal eigenvectors. C Any vector which fails to converge is set to zero. C Z is a two-dimensional REAL array, dimensioned Z(NM,M). C C IERR is an INTEGER flag set to C Zero for normal return, C -J if the eigenvector corresponding to the J-th C eigenvalue fails to converge in 5 iterations. C C RV1, RV2 and RV3 are one-dimensional REAL arrays used for C temporary storage. They are used to store the main diagonal C and the two adjacent diagonals of the triangular matrix C produced in the inverse iteration process. RV1, RV2 and C RV3 are dimensioned RV1(N), RV2(N) and RV3(N). C C RV4 and RV6 are one-dimensional REAL arrays used for temporary C storage. RV4 holds the multipliers of the Gaussian C elimination process. RV6 holds the approximate eigenvectors C in this process. RV4 and RV6 are dimensioned RV4(N) and C RV6(N). C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- C system Routines - EISPACK Guide, Springer-Verlag, C 1976. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 760101 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE TINVIT C INTEGER I,J,M,N,P,Q,R,S,II,IP,JJ,NM,ITS,TAG,IERR,GROUP INTEGER IND(*) REAL D(*),E(*),E2(*),W(*),Z(NM,*) REAL RV1(*),RV2(*),RV3(*),RV4(*),RV6(*) REAL U,V,UK,XU,X0,X1,EPS2,EPS3,EPS4,NORM,ORDER C C***FIRST EXECUTABLE STATEMENT TINVIT IERR = 0 IF (M .EQ. 0) GO TO 1001 TAG = 0 ORDER = 1.0E0 - E2(1) Q = 0 C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX .......... 100 P = Q + 1 C DO 120 Q = P, N IF (Q .EQ. N) GO TO 140 IF (E2(Q+1) .EQ. 0.0E0) GO TO 140 120 CONTINUE C .......... FIND VECTORS BY INVERSE ITERATION .......... 140 TAG = TAG + 1 S = 0 C DO 920 R = 1, M IF (IND(R) .NE. TAG) GO TO 920 ITS = 1 X1 = W(R) IF (S .NE. 0) GO TO 510 C .......... CHECK FOR ISOLATED ROOT .......... XU = 1.0E0 IF (P .NE. Q) GO TO 490 RV6(P) = 1.0E0 GO TO 870 490 NORM = ABS(D(P)) IP = P + 1 C DO 500 I = IP, Q 500 NORM = MAX(NORM, ABS(D(I)) + ABS(E(I))) C .......... EPS2 IS THE CRITERION FOR GROUPING, C EPS3 REPLACES ZERO PIVOTS AND EQUAL C ROOTS ARE MODIFIED BY EPS3, C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .......... EPS2 = 1.0E-3 * NORM EPS3 = NORM 502 EPS3 = 0.5E0*EPS3 IF (NORM + EPS3 .GT. NORM) GO TO 502 UK = SQRT(REAL(Q-P+5)) EPS3 = UK * EPS3 EPS4 = UK * EPS3 UK = EPS4 / UK S = P 505 GROUP = 0 GO TO 520 C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... 510 IF (ABS(X1-X0) .GE. EPS2) GO TO 505 GROUP = GROUP + 1 IF (ORDER * (X1 - X0) .LE. 0.0E0) X1 = X0 + ORDER * EPS3 C .......... ELIMINATION WITH INTERCHANGES AND C INITIALIZATION OF VECTOR .......... 520 V = 0.0E0 C DO 580 I = P, Q RV6(I) = UK IF (I .EQ. P) GO TO 560 IF (ABS(E(I)) .LT. ABS(U)) GO TO 540 C .......... WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF C E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......... XU = U / E(I) RV4(I) = XU RV1(I-1) = E(I) RV2(I-1) = D(I) - X1 RV3(I-1) = 0.0E0 IF (I .NE. Q) RV3(I-1) = E(I+1) U = V - XU * RV2(I-1) V = -XU * RV3(I-1) GO TO 580 540 XU = E(I) / U RV4(I) = XU RV1(I-1) = U RV2(I-1) = V RV3(I-1) = 0.0E0 560 U = D(I) - X1 - XU * V IF (I .NE. Q) V = E(I+1) 580 CONTINUE C IF (U .EQ. 0.0E0) U = EPS3 RV1(Q) = U RV2(Q) = 0.0E0 RV3(Q) = 0.0E0 C .......... BACK SUBSTITUTION C FOR I=Q STEP -1 UNTIL P DO -- .......... 600 DO 620 II = P, Q I = P + Q - II RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I) V = U U = RV6(I) 620 CONTINUE C .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS C MEMBERS OF GROUP .......... IF (GROUP .EQ. 0) GO TO 700 J = R C DO 680 JJ = 1, GROUP 630 J = J - 1 IF (IND(J) .NE. TAG) GO TO 630 XU = 0.0E0 C DO 640 I = P, Q 640 XU = XU + RV6(I) * Z(I,J) C DO 660 I = P, Q 660 RV6(I) = RV6(I) - XU * Z(I,J) C 680 CONTINUE C 700 NORM = 0.0E0 C DO 720 I = P, Q 720 NORM = NORM + ABS(RV6(I)) C IF (NORM .GE. 1.0E0) GO TO 840 C .......... FORWARD SUBSTITUTION .......... IF (ITS .EQ. 5) GO TO 830 IF (NORM .NE. 0.0E0) GO TO 740 RV6(S) = EPS4 S = S + 1 IF (S .GT. Q) S = P GO TO 780 740 XU = EPS4 / NORM C DO 760 I = P, Q 760 RV6(I) = RV6(I) * XU C .......... ELIMINATION OPERATIONS ON NEXT VECTOR C ITERATE .......... 780 DO 820 I = IP, Q U = RV6(I) C .......... IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE C WAS PERFORMED EARLIER IN THE C TRIANGULARIZATION PROCESS .......... IF (RV1(I-1) .NE. E(I)) GO TO 800 U = RV6(I-1) RV6(I-1) = RV6(I) 800 RV6(I) = U - RV4(I) * RV6(I-1) 820 CONTINUE C ITS = ITS + 1 GO TO 600 C .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... 830 IERR = -R XU = 0.0E0 GO TO 870 C .......... NORMALIZE SO THAT SUM OF SQUARES IS C 1 AND EXPAND TO FULL ORDER .......... 840 U = 0.0E0 C DO 860 I = P, Q 860 U = U + RV6(I)**2 C XU = 1.0E0 / SQRT(U) C 870 DO 880 I = 1, N 880 Z(I,R) = 0.0E0 C DO 900 I = P, Q 900 Z(I,R) = RV6(I) * XU C X0 = X1 920 CONTINUE C IF (Q .LT. N) GO TO 100 1001 RETURN END