C VERSION 2 DOES NOT USE EISPACK C C ------------------------------------------------------------------ C SUBROUTINE DILASO(OP, IOVECT, N, XL, XR, NFIG, NPERM, * NMVAL, VAL, NMVEC, MAXVEC, VEC, NBLOCK, MAXOP, MAXJ, * WORK, IND, IERR) C INTEGER N, NFIG, NPERM, NMVAL, NMVEC, MAXVEC, NBLOCK, * MAXOP, MAXJ, IND(1), IERR DOUBLE PRECISION XL, XR, VEC(NMVEC,1), VAL(NMVAL,1), WORK(1) EXTERNAL OP, IOVECT C C AUTHOR/IMPLEMENTER D.S.SCOTT-B.N.PARLETT/D.S.SCOTT C C COMPUTER SCIENCES DEPARTMENT C UNIVERSITY OF TEXAS AT AUSTIN C AUSTIN, TX 78712 C C VERSION 2 ORIGINATED APRIL 1982 C C CURRENT VERSION JUNE 1983 C C DILASO FINDS ALL THE EIGENVALUES AND EIGENVECTORS OF A LARGE C SPARSE SYMMETRIC MATRIX OUTSIDE A USER DEFINED EXCLUDED INTERVAL. C THE SUBROUTINE DILASO IS PRIMARILY A DRIVER FOR SUBROUTINE DIWLA C WHICH IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE C ORTHOGONALIZATION AND SUBROUTINE DIPPLA WHICH POST PROCESSES C THE OUTPUT OF DIWLA. HOWEVER DILASO DOES CHECK FOR INCONSISTENCIES C IN THE CALLING PARAMETERS AND DOES PREPROCESS AND USER SUPPLIED C EIGENPAIRS. C C C ON INPUT C C C OP A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE C OP(N,M,P,Q). P AND Q ARE N X M MATRICES AND Q IS C RETURNED AS THE MATRIX TIMES P. C C IOVECT A USER SUPPLIED SUBROUTINE WITH CALLING SEQUENCE C IOVECT(N,M,Q,J,K). Q IS AN N X M MATRIX. IF K = 0 C THE COLUMNS OF Q ARE STORED AS THE (J-M+1)TH THROUGH C THE JTH LANCZOS VECTORS. IF K = 1 THEN Q IS RETURNED C AS THE (J-M+1)TH THROUGH THE JTH LANCZOS VECTORS. SEE C DOCUMENTATION FOR FURTHER DETAILS AND EXAMPLES. C C N THE ORDER OF THE MATRIX. C C XL THE LEFT ENDPOINT OF THE EXCLUDED INTERVAL. C C XR THE RIGHT ENDPOINT OF THE EXCLUDED INTERVAL. C C NFIG THE NUMBER OF DECIMAL DIGITS OF ACCURACY DESIRED IN THE C EIGENVALUES. NFIG MUST BE GREATER THAN OR EQUAL TO 1. C C NPERM AN INTEGER VARIABLE WHICH SPECIFIES THE NUMBER OF USER C SUPPLIED EIGENPAIRS. IN MOST CASES NPERM WILL BE ZERO. SEE C DOCUMENTAION FOR FURTHER DETAILS OF USING NPERM GREATER C THAN ZERO. NPERM MUST NOT BE LESS THAN ZERO. C C NMVAL THE ROW DIMENSION OF THE ARRAY VAL. NMVAL MUST BE GREATER C THAN OR EQUAL TO MAXVEC. C C VAL A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW C DIMENSION NMVAL AND COLUMN DIMENSION AT LEAST 4. IF NPERM C IS GREATER THAN ZERO THEN CERTAIN INFORMATION MUST BE STORED C IN VAL. SEE DOCUMENTATION FOR DETAILS. C C NMVEC THE ROW DIMENSION OF THE ARRAY VEC. NMVEC MUST BE GREATER C THAN OR EQUAL TO N. C C MAXVEC THE MAXIMUM NUMBER OF EIGENPAIRS TO BE DETERMINED. C MAXVEC MUST NOT EXCEED THE COLUMN DIMENSION OF THE ARRAY VEC. C C VEC A TWO DIMENSIONAL DOUBLE PRECISION ARRAY OF ROW C DIMENSION NMVEC AND COLUMN DIMENSION MAXVEC. IF NPERM > 0 C THEN THE FIRST NPERM COLUMNS OF VEC MUST CONTAIN C THE USER SUPPLIED EIGENVECTORS. C C NBLOCK THE BLOCK SIZE. SEE DOCUMENTATION FOR CHOOSING C AN APPROPRIATE VALUE FOR NBLOCK. NBLOCK MUST BE GREATER C THAN ZERO .AND. LESS THAN MAXJ/6. C C MAXOP AN UPPER BOUND ON THE NUMBER OF CALLS TO THE SUBROUTINE C OP. DILASO TERMINATES WHEN MAXOP IS EXCEEDED. SEE C DOCUMENTATION FOR GUIDELINES IN CHOOSING A VALUE FOR MAXOP. C C MAXJ AN INDICATION OF THE AVAILABLE STORAGE (SEE WORK AND C DOCUMENTATION ON IOVECT). FOR THE FASTEST CONVERGENCE MAXJ C SHOULD BE AS LARGE AS POSSIBLE, ALTHOUGH IT IS USELESS TO HAVE C MAXJ LARGER THAN MAXOP*NBLOCK. C C WORK A DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST AS C LARGE AS C C 2*N*NBLOCK + MAXJ*(NBLOCK+1+MAXVEC+2) + 2*NBLOCK*NBLOCK + 3*MAXVEC C C + THE MAXIMUM OF C N*NBLOCK C AND C MAXJ*(2*NBLOCK+2) + 2*MAXVEC + 8 + (2*NBLOCK+2)*(NBLOCK+1) C C THE FIRST N*NBLOCK ELEMENTS OF WORK MUST CONTAIN THE DESIRED C STARTING VECTORS. SEE DOCUMENTATION FOR GUIDELINES IN C CHOOSING STARTING VECTORS. C C IND AN INTEGER ARRAY OF DIMENSION AT LEAST 4*MAXVEC. C C IERR AN INTEGER VARIABLE. C C C ON OUTPUT C C C NPERM THE NUMBER OF EIGENPAIRS NOW KNOWN. C C VEC THE FIRST NPERM COLUMNS OF VEC CONTAIN THE EIGENVECTORS. C C VAL THE FIRST COLUMN OF VAL CONTAINS THE CORRESPONDING C EIGENVALUES. THE SECOND COLUMN CONTAINS THE RESIDUAL NORMS OF C THE EIGENPAIRS WHICH ARE BOUNDS ON THE ACCURACY OF THE EIGEN- C VALUES. THE THIRD COLUMN CONTAINS MORE DOUBLE PRECISIONISTIC ESTIMATES C OF THE ACCURACY OF THE EIGENVALUES. THE FOURTH COLUMN CONTAIN C ESTIMATES OF THE ACCURACY OF THE EIGENVECTORS. SEE C DOCUMENTATION FOR FURTHER INFORMATION ON THESE QUANTITIES. C C WORK IF EXECUTION IS TERMINATED BEFORE ALL THE DESIRED C EIGENVALUES ARE FOUND (IERR = -2, -5, OR -6 MOD(10) ) THEN C THE FIRST N*NBLOCK ELEMENTS OF WORK CONTAIN THE BEST VECTORS C FOR RESTARTING THE ALGORITHM AND DILASO CAN BE IMMEDIATELY C RECALLED TO CONTINUE WORKING ON THE PROBLEM. C C IND IND(1) CONTAINS THE ACTUAL NUMBER OF CALLS TO OP. ON SOME C OCCASIONS THE NUMBER OF CALLS TO OP MAY BE SLIGHTLY LARGER C THAN MAXOP. C C IERR AN ERROR COMPLETION CODE. THE NORMAL COMPLETION CODE IS C ZERO. SEE THE DOCUMENTATION FOR INTERPRETATIONS OF NON-ZERO C COMPLETION CODES. C C C INTERNAL VARIABLES. C C INTEGER I, I1, I2, I3, I4, I5, I6, I7, I8, I9, I10, I11, * I12, I13, M, NBAND, NOP, NP LOGICAL RARITZ DOUBLE PRECISION DELTAL, DELTAR, EPS, TEMP, DDOT, DNRM2, DABS EXTERNAL DIPPLA, DIWLA, DMVPC, DORTQR, DAXPY, * DCOPY, DDOT, DNRM2, DSCAL, DSWAP, DLAEIG, DLABCM, DLABFC, * DLAGER, DLARAN, DVSORT C C NOP RETURNED FROM DIWLA AS THE NUMBER OF CALLS TO THE C SUBROUTINE OP. C C RARITZ RETURNED FROM DIWLA AS .TRUE. IF A FINAL RAYLEIGH-RITZ C PROCEDURE IS NEEDED. C C DELTAL RETURNED FROM DIWLA AS THE EXCLUDED EIGENVALUE CLOSEST C TO XL. C C DELTAR RETURNED FROM DIWLA AS THE EXCLUDED EIGENVALUE CLOSEST C TO XR. C C DIPPLA A SUBROUTINE FOR POST-PROCESSING THE EIGENVECTORS COMPUTED C BY DIWLA. C C DIWLA A SUBROUTINE FOR IMPLEMENTING THE LANCZOS ALGORITHM C WITH SELECTIVE ORTHOGONALIZATION. C C DMVPC A SUBROUTINE FOR COMPUTING THE RESIDUAL NORM AND C ORTHOGONALITY COEFFICIENT OF A GIVEN RITZ VECTOR. C C DORTQR A SUBROUTINE FOR ORTHONORMALIZING A BLOCK OF VECTORS C USING HOUSEHOLDER REFLECTIONS. C C DAXPY,DCOPY,DDOT,DNRM2,DSCAL,DSWAP A SUBSET OF THE BASIC LINEAR C ALGEBRA SUBPROGRAMS USED FOR VECTOR MANIPULATION. C C DLAEIG, DLABCM, DLABFC, DLAGER USED FOR BAND EIGENVALUE PROBLEMS C C DLARAN SUBROUTINE WHICH RETURNS RANDOM VECTORS C C ------------------------------------------------------------------ C C THIS SECTION CHECKS FOR INCONSISTENCY IN THE INPUT PARAMETERS. C IND(1) = 0 IERR = 0 C IF (N.LT.6*NBLOCK) IERR = 1 IF (NFIG.LE.0) IERR = IERR + 2 IF (NMVEC.LT.N) IERR = IERR + 4 IF (NPERM.LT.0) IERR = IERR + 8 IF (MAXJ.LT.6*NBLOCK) IERR = IERR + 16 IF (MAXVEC.LT.NPERM) IERR = IERR + 32 IF (MAXVEC.GT.NMVAL) IERR = IERR + 64 IF (MAXVEC.GT.MAXOP) IERR = IERR + 128 IF (XL.GT.XR) IERR = IERR + 256 IF (NBLOCK.LT.1) IERR = IERR + 512 IF (IERR.NE.0) RETURN C C ------------------------------------------------------------------ C C THIS SECTION SORTS AND ORTHONORMALIZES THE USER SUPPLIED VECTORS. C IF A USER SUPPLIED VECTOR IS ZERO OR IF DSIGNIFICANT CANCELLATION C OCCURS IN THE ORTHOGONALIZATION PROCESS THEN IERR IS SET TO -1 C AND DILASO TERMINATES. IF A USER SUPPLIED EIGENVALUE LIES INSIDE C THE EXCLUDED INTERVAL THEN IERR IS SET TO -4 AND DILASO TERMINATES. C IF (NPERM.EQ.0) GO TO 100 C C THIS SORTS THE USER SUPPLIED VALUES AND VECTORS. C CALL DVSORT(NPERM, VAL, VAL(1,2), 0, TEMP, NMVEC, N, VEC) C C THIS STORES THE NORMS OF THE VECTORS FOR LATER COMPARISON. C IT ALSO INSURES THAT THE REIDUAL NORMS ARE POSITIVE. C DO 40 I=1,NPERM VAL(I,2) = DABS(VAL(I,2)) VAL(I,3) = DNRM2(N, VEC(1,I), 1) 40 CONTINUE C C THIS CHECKS FOR VALUES OUT OF RANGE C DO 50 I=1,NPERM IF (VAL(I,1).LE.XL .OR. VAL(I,1).GE.XR) GO TO 50 IERR = -4 RETURN C 50 CONTINUE C C THIS ORTHONORMALIZES THE VECTORS C M = N*NBLOCK + 1 CALL DORTQR(NMVEC, N, NPERM, VEC, WORK(M)) M = N*NBLOCK - NPERM DO 60 I = 1, NPERM M = M + NPERM + 1 IF(DABS(WORK(M)) .GT. 0.9*VAL(I,3)) GO TO 60 IERR = -1 RETURN C 60 CONTINUE C C THIS COPIES THE RESIDUAL NORMS INTO THE CORRECT LOCATIONS IN C THE ARRAY WORK FOR LATER REFERENCE IN DIWLA. C M = 2*N*NBLOCK + 1 CALL DCOPY(NPERM, VAL(1,2), 1, WORK(M), 1) C C THIS SETS EPS TO AN APPROXIMATION OF THE MACHINE PRECISION C C ***THIS SHOULD BE REPLACED BY AN ASDSIGNMENT STATEMENT C ***IN A PRODUCTION CODE C 100 EPS = 1.0D0 DO 110 I =1,1000 EPS = 0.5D0*EPS TEMP = 1.0D0 + EPS IF(TEMP.EQ.1.0D0) GO TO 120 110 CONTINUE C C ------------------------------------------------------------------ C C THIS SECTION CALLS DIWLA WHICH IMPLEMENTS THE LANCZOS ALGORITHM C WITH SELECTIVE ORTHOGONALIZATION. C 120 NBAND = NBLOCK + 1 I1 = 1 + N*NBLOCK I2 = I1 + N*NBLOCK I3 = I2 + MAXVEC I4 = I3 + MAXVEC I5 = I4 + MAXVEC I6 = I5 + MAXJ*NBAND I7 = I6 + NBLOCK*NBLOCK I8 = I7 + NBLOCK*NBLOCK I9 = I8 + MAXJ*(MAXVEC+2) I10 = I9 + NBLOCK I11 = I10 + 2*MAXVEC + 8 I12 = I11 + MAXJ*(2*NBLOCK+1) I13 = I12 + MAXJ CALL DIWLA(OP, IOVECT, N, NBAND, XL, XR, NFIG, NPERM, VAL, * NMVEC, MAXVEC, VEC, NBLOCK, MAXOP, MAXJ, NOP, * WORK(1), WORK(I1), WORK(I2), WORK(I3), WORK(I4), * WORK(I5), WORK(I6), WORK(I7), WORK(I8), WORK(I9), * WORK(I10), WORK(I11), WORK(I12), WORK(I13), * IND, RARITZ, DELTAL, DELTAR, EPS, IERR) C C ------------------------------------------------------------------ C C THIS SECTION CALLS DIPPLA (THE POST PROCESSOR). NP IS USED C TO DIMENSION AN ARRAY IN DIPPLA. IT IS NEEDED TO LEGALIZE C CHANGING THE VALUE OF NPERM, WHICH MAY OCCUR IN DIPPLA. C IF (NPERM.EQ.0) GO TO 130 I1 = N*NBLOCK + 1 I2 = I1 + NPERM*NPERM I3 = I2 + NPERM*NPERM I4 = I3 + MAX0(N*NBLOCK, 2*NPERM*NPERM) I5 = I4 + N*NBLOCK I6 = I5 + 2*NPERM + 4 NP = NPERM CALL DIPPLA(OP, IOVECT, N, XL, XR, NP, NPERM, NOP, NMVAL, * VAL, NMVEC, VEC, NBLOCK, WORK(I1), WORK(I2), * WORK(I3), WORK(I4), WORK(I5), WORK(I6), DELTAL, * DELTAR, RARITZ, EPS, IERR) C 130 IND(1) = NOP RETURN END