* Date: Tue, 21 Jul 92 15:49:01 -0700 * From: Cheryl Carey C***************************************************************** C***************************************************************** C AN ITERATIVE BLOCK LANCZOS METHOD FOR THE C SOLUTION OF LARGE SPARSE SYMMETRIC C EIGENPROBLEMS C by C Richard Ray Underwood C C C if further details concerning the content of the code are C required, you are advised to refer to R. R. Underwood's C PhD thesis - "An Iterative Block Lanczos Method for the C Solution of Large Sparse Symmetric Eigenproblems", 1975. C***************************************************************** C***************************************************************** IMPLICIT REAL*8 (A-H,O-Z) C================================================================= C sample main program. The values for N, Q, PINIT, R, MMAX and C EPS are all problem dependent and should be altered C accordingly, as too should the dimensions of the arrays D C and X. C================================================================= DIMENSION D(25),X(20000) EXTERNAL AX INTEGER Q,PINIT,R N = 300 Q=12 PINIT =4 R=4 MMAX=5000 EPS=1D-03 M=0 C C------ C MINVAL is used to compute the 4 least eigenvalues C of the matrix A = diag(-1,-1/2,...,-1/300) to an C approximate precision of 10**(-3). 12 vectors are C allowed for the block lanczos method and an initial C block size of 4 is chosen. C------ CALL MINVAL(N,Q,PINIT,R,MMAX,EPS,AX,M,D,X,IECODE) C PRINT 1001,M,IECODE,(D(I),I=1,M) 1001 FORMAT(/' =>M,IECODE=',I4,1X,I4/(' =>E ',5D23.15)) STOP END C C C C SUBROUTINE MINVAL(N,Q,PINIT,R,MMAX,EPS,OP,M,D,X,IECODE) C================================================================= C this subroutine is the main subroutine implementing the C iterative block lanczos method for computing the C eigenvalues and eigenvectors of symmetric matrices. C================================================================= IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,PINIT,R,MMAX REAL*8 EPS EXTERNAL OP DIMENSION D(Q),X(N,Q) INTEGER IECODE DIMENSION E(25),C(25,25) DIMENSION U(1296),V(1296) INTEGER P,S,PS C C------ C description of parameters: C N: integer variable. The order of the symmetric C matrix A whose eigenvalues and eigenvectors are C being computed. The value of N should be less than C or equal to 1296 and greater than or equal to 2. C C Q: integer variable. The number of vectors of length C N contained in the array X. The value of Q should C be less than or equal to 25, at least one greater C than the value of R and less than or equal to N. C C PINIT: integer variable. The initial block size to be used C in the block lanczos method. If PINIT is negative, C then -PINIT is used for the block size and columns C M+L,...,M+(-PINIT) of the array X are assumed to be C initialized to the initial matrix used to start the C block lanczos method. If the subroutine terminates C with a value of M less than R, then PINIT is assigned C a value -P where P is the final block size chosen. C In this circumstance, columns M+1,...M+P will contain C the most recent set of eigenvector approximations C which can be used to restart the subroutine if desired. C C R: integer variable. The number of eigenvalues and C eigenvectors being computed. That is, MINVAL attempts C to compute accurate approximations to the R least C eigenvalues and eigenvectors of the matrix A. The value C of R should be greater than zero and less than Q. C C MMAX: integer variable. The maximum number of matrix-vector C products A*X (where X is a vector) that are allowed C during one call of this subroutine to complete its task C of computing R eigenvalues and eigenvectors. Unless the C problem indicates otherwise, MMAX should be given a very C large value. C C EPS: REAL*8 variable. Initially, EPS should contain a value C indicating the relative precision to which MINVAL will C attempt to compute the eigenvalues and eigenvectors of A. C For eigenvalues less in modulus than 1, EPS will be an C absolute tolerance. Because of the way this method works, C it may happen that the later eigenvalues cannot be C computed to the same relative precision as those less in C value. C C OP: subroutine name. The actual argument corresponding to OP C should be the name of a subroutine used to define the C matrix A. This subroutine should have three arguments C N, U, and V, say, where N is an integer variable giving C the order of A, and U and V are two one-dimensional C arrays of length N. If W denotes the vector of order N C stored in U, then the statement C CALL OP(N,U,V) C should result in the vector A*W being computed and stored C in V. The contents of U can be modified by this call. C C M: integer variable. M gives the number of eigenvalues and C eigenvectors already computed. Thus, initially, M should C be zero. If M is greater than zero, then columns one C through M of the array X are assumed to contain the C computed approximations to the M least eigenvalues and C eigenvectors of A. On exit, M contains a value equal to C the total number of eigenvalues and eigenvectors C computed including any already computed when MINVAL was C entered. Thus, on exit, the first M elements of D and the C first M columns of X will contain approximations to the C M least eigenvalues of A. C C D: REAL*8 array. D contains the computed eigenvalues. D C should be a one-dimensional array with at least Q C elements. C C X: REAL*8 array. X contains the computed eigenvectors. X C should be an array containing at least N*Q elements. X C is used not only to store the eigenvectors computed by C MINVAL, but also as working storage for the block lanczos C method. On exit, the first N*M elements of X contain the C eigenvector approximations - the first vector in the C first N elements, the second in the second N elements, C etc... C C C IECODE:integer variable. The value of IECODE indicates whether C MINVAL terminated successfully, and if not, the reason C why. C C IECODE=0 : successful termination. C IECODE=1 : the value of N is less than 2. C IECODE=2 : the value of N exceeds 1296. C IECODE=3 : the value of R is less than 1. C IECODE=4 : the value of Q is less than or equal to R. C IECODE=5 : the value of Q is greater than 25. C IECODE=6 : the value of Q exceeds N. C IECODE=7 : the value of MMAX was exceeded before R C eigenvalues and eigenvectors were C computed. C C Note that the subroutine has been designed to allow initial C approximations to the eigenvectors corresponding to the least C eigenvalues to be utilised if they are known (by storing them C in X and assigning PINIT minus the value of their number). C Furthermore, it has also been designed to allow restarting if C it stops with IECODE=7. Thus, the user of this program can C restart it after examining any partial results without loss of C previous work. C------ C C------ C check that the initial values of the subroutine parameters C are in range. C------ IF (N.LT.2) GO TO 901 IF (N.GT.1296) GO TO 902 IF (R.LT.1) GO TO 903 IF (Q.LE.R) GO TO 904 IF (Q.GT.25) GO TO 905 IF (Q.GT.N) GO TO 906 C C------ C choose initial values for the block size P, the number of C steps that the block lanczos method is carried out, and C choose an initial N-by-P orthonormal matrix X1 used to start C the block lanczos method. C------ P = PINIT IF (P.LT.0) P = -P S = (Q-M)/P IF (S.GT.2) GO TO 100 S = 2 P = Q/2 100 IF (PINIT.LT.0) GO TO 150 DO 120 K = 1,P CALL RANDOM(N,Q,K,X) 120 CONTINUE 150 IF (M.GT.0) GO TO 200 CALL ORTHG(N,Q,M,P,C,X) C C------ C rotate the initial N-by-P matrix X1 so that C X1'*A*X1=diag(D1,D2,...,DP) C where DI is stored in D(I), I=1,...,P. C------ CALL SECTN(N,Q,M,P,OP,X,C,D,U,V) ERRC = 0.D0 200 ITER = 0 IMM=0 C C------ C the main body of the subroutine starts here. IMM C counts the number of matrix-vector products computed, C which is the number of times the subroutine named by C OP is called. ERRC measures the accumulated error in C the eigenvalues and eigenvectors. C------ C 300 IF (M.GE.R) GO TO 900 IF (IMM.GT.MMAX) GO TO 907 ITER = ITER + 1 PS = P*S C C------ C BKLANC carries out the block lanczos method and stores C the resulting block tridiagonal matrix MS in C and the C N-by-PS orthonormal matrix XS in X. The initial N-by-P C orthonormal matrix is assumed to be stored in columns C M+1 through M+PS of X. The residuals for these vectors C and the eigenvalue approximations in D are computed and C stored in E. C------ CALL BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V) C C------ C EIGEN solves the eigenproblem for MS, storing the eigenvalues C in elements M+1 through M+PS of D and the eigenvectors in the C first P*S rows and columns of C (overwriting MS, possibly). C------ CALL EIGEN(Q,M,P,PS,C,D) C C------ C CNVTST determines how many of the eigenvalues and eigenvectors C have converged using the error estimates stored in E. The number C that have converged is stored in NCONV. If NCONV=0, then none C have converged. C------ CALL CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV) write(6,*) nconv C C------ C PCH chooses new values for P and S, the block size and the C number of steps for the block lanczos subprogram, respectively. C------ CALL PCH(N,Q,M,R,NCONV,P,S) C C------ C ROTATE computes the eigenvectors of the restricted matrix C using XS stored in X and the eigenvectors of MS stored in C. C These vectors serve both as eigenvector approximations and C to form the matrix used to start the block lanczos method in C the next iteration. C------ CALL ROTATE(N,Q,M,PS,NCONV+P,C,X) C M = M + NCONV IMM=IMM+P*S PRINT 1001,ITER,IMM,P,PS 1001 FORMAT(' =>ITER,IMM,P,PS =',4I5) GO TO 300 C C------ C this is the end of the main body of the subroutine. Now set C the value of IECODE and EXIT. C------ 900 IECODE = 0 RETURN 901 IECODE = 1 RETURN 902 IECODE = 2 RETURN 903 IECODE = 3 RETURN 904 IECODE = 4 RETURN 905 IECODE = 5 RETURN 906 IECODE = 6 RETURN 907 IECODE = 7 PINIT = -P C RETURN END C C C C SUBROUTINE BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V) C==================================================================== C this subroutine implements the block lanczos method C with reorthogonalization. BKLANC computes a block C tridiagonal matrix MS which it stores in rows and C columns M+1 through M+P*S of the array C, and an C orthonormal matrix XS which it stores in columns M+1 C through M+P*S of the N-by-Q array X. MS is a symmetric C matrix with P-by-P symmetric matrices M(1),...,M(S) on C its diagonal and P-by-P upper triangular matrices C R(2),...,R(S) along its lower diagonal. Since MS is C symmetric and banded, only its lower triangle (P+1 C diagonals) is stored in C. XS is composed of S N-by-P C orthonormal matrices X(1),...,X(S) where X(1) is given C and should be stored in columns M+1 through M+P of X. C Furthermore, X(1) is assumed to satisfy X(1)*A*X(1) = C diag(D(M+1),D(M+2),...,D(M+P)), and if M>0, then X(1) is C assumed to be orthogonal to the vectors stored in columns C 1 through M of X. OP is the name of an external subroutine C used to define the matrix A. During the first step, the C subroutine ERR is called and the quantities EJ are computed C where EJ=||A*X1J-D(M+J)*X1J||, X1J is the J-th column of x(1), C and ||*|| denotes the Euclidean norm. EJ is stored in E(M+J), C J=1,2,...,P. U and V are auxilliary vectors used by OP. C==================================================================== IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P,S DIMENSION D(Q),C(Q,Q),X(N,Q) DIMENSION E(Q),U(N),V(N) C MP1 = M+1 MPPS = M+P*S DO 90 L = 1,S LL = M+(L-1)*P+1 LU = M + L*P DO 70 K = LL,LU DO 10 I = 1,N 10 U(I) = X(I,K) CALL OP(N,U,V) IF (L.GT.1) GO TO 19 DO 12 I=K,LU 12 C(I,K) = 0.D0 C(K,K) = D(K) DO 14 I = 1,N 14 V(I) = V(I) - D(K)*X(I,K) GO TO 61 19 DO 30 I = K,LU T = 0 DO 20 J = 1,N 20 T = T+V(J)*X(J,I) 30 C(I,K) = T IT = K-P DO 60 I = 1,N T = 0 DO 40 J = IT,K 40 T = T + X(I,J)*C(K,J) IF (K.EQ.LU) GO TO 60 KP1 = K+1 DO 50 J = KP1,LU 50 T = T+X(I,J)*C(J,K) 60 V(I) = V(I)-T 61 IF (L.EQ.S) GO TO 70 DO 63 I = 1,N 63 X(I,K+P) = V(I) 70 CONTINUE IF (L.EQ.1) CALL ERR(N,Q,M,P,X,E) IF (L.EQ.S) GO TO 90 CALL ORTHG(N,Q,LU,P,C,X) IL = LU+1 IT = LU DO 80 J = 1,P IT = IT+1 DO 80 I = IL,IT 80 C(I,IT-P) = C(I,IT) 90 CONTINUE C RETURN END C C C C SUBROUTINE PCH(N,Q,M,R,NCONV,P,S) C==================================================================== C based on the values of N, Q, M, R and NCONV, PCH chooses new C values for P and S, the block size and number of steps for the C block lanczos method. The strategy used here is to choose P to C be the smaller of the two following values: C 1) the previous block size C and, 2) the number of values left to be computed. S is chosen C as large as possible subject to the constraints imposed by the C limits of storage. In any event, S is greater than or equal to C 2. N is the order of the problem and Q is the number of vectors C available for storing eigenvectors and applying the block C lanczos method. M is the number of eigenvalues and eigenvectors C that have already been computed and R is the required number. C Finally, NCONV is the number of eigenvalues and eigenvectors C that have converged in the current iteration. C==================================================================== IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,R,NCONV,P,S INTEGER PT,ST C MT = M+NCONV PT = R-MT IF (PT.GT.P) PT = P IF (PT.GT.0) GO TO 101 P = 0 RETURN 101 CONTINUE ST = (Q-MT)/PT IF (ST.GT.2) GO TO 110 ST = 2 PT = (Q-MT)/2 110 P = PT S = ST C RETURN END C C C C SUBROUTINE ERR(N,Q,M,P,X,E) C================================================ C ERR computes the Euclidean lengths of the C vectors stored in the columns M+P+1 through C M+P+P of the N-by-Q array X and stores them C in elements M+1 through M+P of E. C================================================ IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P DIMENSION X(N,Q),E(Q) C MP1 = M+P+1 MPP = M+P+P DO 200 K = MP1,MPP T = 0.D0 DO 100 I = 1,N 100 T = T+X(I,K)**2 200 E(K-P) = DSQRT(T) C RETURN END C C C C SUBROUTINE CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV) C=========================================================== C CNVTST determines which of the P eigenvalues stored C in elements M+1 through M+P of D have converged. ERRC C is a measure of the accumulated error in the M C previously computed eigenvalues and eigenvectors. ERRC C is updated if more approximations have converged. The C norms of the residual vectors are stored in elements C M+1 through M+P of E. EPS is the precision to which we C are computing the approximations. Finally, NCONV is the C number that have converged. If NCONV=0, then none have C converged. C============================================================ IMPLICIT REAL*8 (A-H,O-Z) INTEGER Q,M,P REAL*8 EPS DIMENSION D(Q),E(Q) DATA CHEPS/2.22D-16/ C K = 0 DO 100 I = 1,P T = DABS(D(M+I)) IF (T.LT.1.D0) T = 1.D0 IF (E(M+I).GT.T*(EPS+10.D0*N*CHEPS)+ERRC) GO TO 200 100 K = I 200 NCONV = K IF (K.EQ.0) RETURN T = 0.D0 DO 300 I = 1,K T = T+E(M+I)**2 300 CONTINUE ERRC = DSQRT(ERRC**2+T) C RETURN END C C C C SUBROUTINE EIGEN(Q,M,P,PS,C,D) C======================================================= C EIGEN solves the eigenproblem for the symmetric C matrix MS of order PS stored in rows and columns C M+1 through M+PS of C. The eigenvalues of MS are C stored in elements M+1 through M+PS of D and the C eigenvactors are stored in rows and columns 1 C through PS of C possibly overwriting MS. EIGEN C simply re-stores MS in a manner acceptable to the C subroutines TRED2 and TQL2. These two routines are C available through eispack. C======================================================= IMPLICIT REAL*8 (A-H,O-Z) INTEGER M,P,Q,PS,PP1 DIMENSION C(Q,Q),D(Q) DIMENSION DD(25),V(25) C DO 150 I=1,PS LIM = I-P IF (I.LE.P) LIM = 1 IF (LIM.LE.1) GO TO 130 LM1 = LIM-1 DO 120 J = 1,LM1 120 C(I,J) = 0.D0 130 DO 140 J = LIM,I 140 C(I,J) = C(I+M,J+M) 150 CONTINUE C C------ C call eispack routines here C------ CALL TRED2(Q,PS,C,DD,V,C) CALL TQL2(Q,PS,DD,V,C,IERR) C PRINT 601,PS,(DD(I),I=1,PS) 601 FORMAT(' => ORDER =',I4,/,(' EIGENVALUES =',10d10.4)) PP1 = P+1 DO 155 J = 1,PP1 155 PRINT 602,J,(C(I,J),I=1,PS) 602 FORMAT(' => J =',I4,/,(' EIGENVECTORS =',10d10.4)) DO 160 I = 1,PS 160 D(M+I) = DD(I) C RETURN END C C C C SUBROUTINE SECTN(N,Q,M,P,OP,X,C,D,U,V) C============================================================== C SECTN transforms the N-by-P orthonormal matrix X1, C say, stored in columns M+1 through M+P of the N-by-Q C array X so that X1'*A*X1 = diag(D1,D2,...,DP), where C ' denotes transpose and A is a symmetric matrix of C order N defined by the subroutine OP. The values D1,... C ,DP are stored in elements M+1 through M+P of D. SECTN C forms the matrix X1'*A*X1 = CP, storing CP in the array C C. The values D1,D2,...,DP and the eigenvectors QP of CP C are computed by EIGEN and stored in D and C respectively. C ROTATE then carries out the transformation X1<=X1*QP. C============================================================== IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P EXTERNAL OP DIMENSION X(N,Q),C(Q,Q),D(Q),U(N),V(N) C ICOL1 = M DO 300 J = 1,P ICOL1 = ICOL1 + 1 DO 100 I = 1,N 100 U(I) = X(I,ICOL1) CALL OP(N,U,V) ICOL2 = M DO 300 I = 1,J ICOL2 = ICOL2+1 T = 0.D0 DO 200 K = 1,N 200 T = T + V(K)*X(K,ICOL2) 300 C(ICOL1,ICOL2) = T CALL EIGEN(Q,M,P,P,C,D) CALL ROTATE(N,Q,M,P,P,C,X) C RETURN END C C C C SUBROUTINE ROTATE(N,Q,M,PS,L,C,X) C====================================================== C ROTATE computes the first L columns of the matrix C XS*QS where XS is an N-by-PS orthonormal matrix C stored in columns M+1 through M+PS of the N-by-Q C array X and QS is a PS-by-PS orthonormal matrix C stored in rows and columns 1 through Ps of the C array C. The result is stored in columns M+1 C through M+L of X overwriting part of XS. C====================================================== IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,PS,L DIMENSION C(Q,Q),X(N,Q) DIMENSION V(25) C DO 300 I = 1,N DO 200 K = 1,L T = 0.D0 DO 100 J = 1,PS 100 T = T+X(I,M+J)*C(J,K) 200 V(K) = T DO 300 K = 1,L 300 X(I,M+K) = V(K) C RETURN END C C C C SUBROUTINE ORTHG(N,Q,F,P,B,X) C========================================================= C ORTHG reorthogonalizes the N-by-P matrix Z stored C in columns F+1 through F+P of the N-by-Q array X C with respect to the vectors stored in the first F C columns of X and then decomposes the resulting C matrix into the product of an N-by-P orthonormal C matrix XORTH, say, and a P-by-P upper triangular C matrix R. XORTH is stored over Z and the upper C triangle of R is stored in rows and columns F+1 C through F+P of the Q-by-Q array B. A stable variant C of the Gram-Schmidt orthogonalization method is C utilised. C========================================================= IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,F,P DIMENSION B(Q,Q),X(N,Q) INTEGER FP1,FPP LOGICAL ORIG C IF (P.EQ.0) RETURN FP1 = F+1 FPP = F+P DO 50 K = FP1,FPP ORIG = .TRUE. KM1 = K-1 10 T = 0.D0 IF (KM1.LT.1) GO TO 25 DO 20 I = 1,KM1 S = 0.D0 DO 15 J = 1,N 15 S = S+X(J,I)*X(J,K) IF (ORIG.AND.I.GT.F) B(I,K) = S T = T + S*S DO 20 J = 1,N 20 X(J,K) = X(J,K) - S*X(J,I) 25 S = 0.D0 DO 30 J = 1,N 30 S = S + X(J,K)*X(J,K) T = T+S IF (S.GT.T/100.D0) GO TO 40 ORIG = .FALSE. GO TO 10 40 S = DSQRT(S) B(K,K) = S IF (S.NE.0) S = 1.D0/S DO 50 J = 1,N 50 X(J,K) = S*X(J,K) C RETURN END C C C C SUBROUTINE RANDOM(N,Q,L,X) C====================================================== C RANDOM computes and stores a sequence of N C pseudo-random numbers in the L-th column of the C N-by-Q array X. RANDOM generates two sequences of C pseudo-random numbers, filling an array with one C sequence and using the second to access the array C in a random fashion. C====================================================== IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,L,FT,X1,F1,F2 INTEGER A,C,X0 DIMENSION X(N,Q) DIMENSION T(100) DATA F1/71416/,F2/27183/ DATA A/6821/,C/5327/,X0/5328/ C DO 100 I = 1,100 X1 = A *X0+C IF (X1.GE.10000) X1 = X1 - 10000 T(I) = X1/9999.D0 - 0.5D0 100 X0 = X1 DO 200 I = 1,N FT = F1+F2 IF (FT.GE.1000000) FT = FT-1000000 F1 = F2 F2 = FT K = FT/1.D6*100 + 1 X(I,L) = T(K) X1 = A*X0 + C IF (X1.GE.10000) X1 = X1 - 10000 T(K) = X1/9999D0 - 0.5D0 200 X0 = X1 C RETURN END C C C C SUBROUTINE AX(N,U,V) C====================================================== C user dependent subroutine. It is defined to be C the actual argument corresponding to OP and is C used to define the matrix A. This subroutine C should have three arguments N, U and V, say, C where N is an integer variable giving the order C of A, and U and V are two one-dimensional arrays C of length N. C====================================================== IMPLICIT REAL*8(A-H,O-Z) INTEGER N DIMENSION U(N),V(N) C C------ C in this example, the matrix A is defined to be C A = diag(-1,-1/2,...,-1/300). C------ DO 100 I=1,N V(I)=(-1.D0/I)*U(I) 100 CONTINUE C RETURN END C C C C