C SAMPLE DRIVER PROGRAM FOR THE LASO EIGENVALUE EXTRACTION PACKAGE. C THE FOLLOWING PROGRAM CALLS THE LASO PACKAGE 3 TIMES TO SOLVE C A SIMPLE TEST PROBLEM. SINCE THIS PROGRAM WAS EXECUTED ON C A DEC 20, IT IS WRITTEN IN DOUBLE PRECISION AND CALLS THE C DOUBLE PRECISION VERSION OF LASO. THE MATRIX USED IS C C A = DIAG(2, 2, 4, 6, 8, 20, 21, 22, ..., 114) C C OF DIMENSION 100. THE FIVE SMALLEST EIGENVALUES ARE COMPUTED IN C ALL THREE RUNS. THE FIRST CALL IS TO DNLASO TO COMPUTE THE C FIVE SMALLEST EIGENVALUES DIRECTLY. THE SECOND CALL IS TO DILASO C TO COMPUTE ALL THE EIGENVALUES LESS THAN 10 AND GREATER THAN C 1000, WHERE 1000 WAS CHOSEN SO THAT THERE WERE NO EIGENVALUES C BIGGER THAN IT. THIS APPROACH IS HANDY WHEN THE NUMBER OF C EIGENVALUES LESS THAN 10 IS UNKNOWN ALTHOUGH A SMALL OVERHEAD C IS PAID FOR PERIODICALLY LOOKING FOR EIGENVALUES BIGGER THAN 1000. C THE THIRD CALL IS TO DILASO WITH A MODIFIED MATRIX MULTIPLY C SUBROUTINE (OP). IT IS ASSUMED THAT A-5I HAS BEEN FACTORED C SO THAT THE OP SUBROUTINE USES (A-5I) INVERSE. THIS TRANSFORMS C THE EIGENVALUES SO THAT THE DESIRED EIGENVALUES (THOSE BETWEEN C ZERO AND TEN) ARE NOW THOSE LESS THAN -0.2 AND GREATER THAN C +0.2. THIS TRANSFORMATION ACCELERATES THE CONVERGENCE C SIGNIFICANTLY AT THE COST OF A MATRIX FACTORIZATION. WHETHER THIS C IS WORTHWHILE WILL DEPEND ON THE RELATIVE COSTS. C C INPUT PARAMETERS TO LASO. C C OP: MATRIX MULTIPLY SUBROUTINE. IN THE FIRST TWO CALLS THIS C JUST MULTIPLIES BY A. C C OP1: MATRIX MULTIPLY SUBROUTINE FOR THE THIRD CALL. IT USES C (A-5I) INVERSE. C C IOVECT: SUBROUTINE FOR STORING LANCZOS VECTORS. THE ONE USED C HERE STORES THE VECTORS IN CORE. C C N: THE DIMENSION OF THE MATRIX (100). C C NVAL: -5 (FIND THE FIVE SMALLEST) USED IN DNLASO ONLY. C C XL: THE LEFT ENDPOINT FOR DILASO (10 IN THE FIRST CALL AND C -0.2 IN THE SECOND). C C XR: THE RIGHT ENDPOINT FOR DILASO (1000 IN THE FIRST CALL AND C +0.2 IN THE SECOND). C C NFIG: THE NUMBER OF DECIMAL DIGITS OF ACCURACY REQUESTED IN C THE EIGENVALUES (5 IN ALL CASES). C C NPERM: THE NUMBER OF KNOWN EIGENPAIRS ON ENTRY (ZERO IN ALL CASES). C C NMVAL: THE ROW DIMENSION OF VAL (10 IN ALL CASES). C C VAL: ARRAY THAT RETURNS EIGENVALUES AND ERROR ESTIMATES. C C NMVEC: ROW DIMENSION OF VEC (110 IN ALL CASES). C C MAXVEC: COLUMN DIMENSION OF VEC (10, USED ONLY BY DILASO). C C VEC: ARRAY THAT RETURNS THE EIGENVECTORS. C C NBLOCK: THE BLOCK SIZE (SET TO THREE IN ALL CASES) C C MAXOP: THE MAXIMUM NUMBER OF MATRIX MULTIPLIES ALLOWED BEFORE C TERMINATION. (SET TO 40 IN ALL CASES). C C MAXJ: THE MAXIMUM SIZE OF A KRYLOV SUBSPACE USED (SEE WORK) C SET TO 50 IN ALL CASES. C C WORK: WORK SPACE. THE FIRST N*NBLOCK LOCATIONS ARE SET TO C THE STARTING VECTORS. IN ALL CASES THE FIRST VECTOR C IS SET TO ALL 1'S AND THE OTHER TWO ARE SET TO ZERO C TO BE REPLACED BY RANDOM VECTORS. THE MINIMUM SIZE C OF WORK IS C C NBLOCK*(3*N + 2*NBLOCK) + MAXJ*(3*NBLOCK + ABS(NVAL) + 6) C + 3*ABS(NVAL) C C FOR DNLASO (WHICH IS 1933) AND C C NBLOCK*(3*N + 2*NBLOCK) + MAXJ*(3*NBLOCK + MAXVEC + 6) C + 3*MAXVEC C C FOR DILASO (WHICH IS 2280). C C IND: INTEGER WORK ARRAY. C C IERR: ERROR RETURN CODE. C LOGICAL BOOL(6) COMMON/BB/BOOL DOUBLE PRECISION A(100),VAL(10,4),VEC(110,10),WORK(2500), 1 QS(110,50), XL, XR DIMENSION IND(20) COMMON /MULT/A COMMON /STORE/QS EXTERNAL OP,OP1,IOVECT C C SET UP A ARRAY A(1)=2.0D0 A(2)=2.0D0 A(3)=4.0D0 A(4)=6.0D0 A(5)=8.0D0 DO 10 I=6,100 A(I)=DBLE(FLOAT(14+I)) 10 CONTINUE C C SET THE FIXED PARAMETERS N=100 NFIG=5 NMVAL=10 NMVEC=110 NBLOCK=3 MAXOP=40 MAXJ=50 C C SET PARAMETERS FOR CALL TO DNLASO NVAL=-5 NPERM = 0 C C SET STARTING VECTORS DO 20 I=1,100 WORK(I)=1.0D0 WORK(N+I)=0.0D0 WORK(2*N+I)=0.0D0 20 CONTINUE C C CALL DNLASO CALL DNLASO(OP,IOVECT,N,NVAL,NFIG,NPERM,NMVAL,VAL,NMVEC, 1 VEC,NBLOCK,MAXOP,MAXJ,WORK,IND,IERR) C C WRITE OUTPUT WRITE(6,1000)IERR,IND(1) 1000 FORMAT(//'1 IERR =',I5,' NUMBER OF MATRIX ACCESSES =',I6) WRITE(6,2000)((VAL(I,J),J=1,4),I=1,NPERM) 2000 FORMAT(//7X,' EIGENVALUE',10X,'RESIDUAL NORM',2X,'VALUE ERROR', 1 4X,'VECTOR ERROR'//(D25.15,3D15.5)) WRITE(6,3000)(J,(VEC(I,J),I=1,6),J=1,NPERM) 3000 FORMAT(//' FIRST SIX COMPONENTS OF THE EIGENVECTORS'// 1 (I5,6D12.3)) C C SET PARAMETERS FOR FIRST CALL TO DILASO XL=10.0D0 XR=1000.0D0 NPERM=0 MAXVEC=10 DO 30 I=1,100 WORK(I)=1.0D0 WORK(N+I)=0.0D0 WORK(2*N+I)=0.0D0 30 CONTINUE CALL DILASO(OP,IOVECT,N,XL,XR,NFIG,NPERM,NMVAL,VAL,NMVEC, 1 MAXVEC,VEC,NBLOCK,MAXOP,MAXJ,WORK,IND,IERR) C C WRITE OUTPUT WRITE(6,1000)IERR,IND(1) WRITE(6,2000)((VAL(I,J),J=1,4),I=1,NPERM) WRITE(6,3000)(J,(VEC(I,J),I=1,6),J=1,NPERM) C C SET PARAMETERS FOR SECOND CALL TO DILASO NPERM=0 XL=-0.2D0 XR=+0.2D0 DO 60 I=1,N WORK(I)=1.0D0 WORK(N+I)=0.0D0 WORK(2*N+I)=0.0D0 60 CONTINUE CALL DILASO(OP1,IOVECT,N,XL,XR,NFIG,NPERM,NMVAL,VAL,NMVEC, 1 MAXVEC,VEC,NBLOCK,MAXOP,MAXJ,WORK,IND,IERR) C WRITE(6,1000)IERR,IND(1) WRITE(6,4000) 4000 FORMAT(//'********* NOTE *********'/ 1 ' ALL THE OUTPUT HERE REFERS TO THE TRANSFORMED PROBLEM'/ 2 ' IT IS NECESSARY TO BACK TRANSFORM THE EIGENVALUES TO'/ 3 ' OBTAIN THE EIGENVALUES OF THE ORIGINAL PROBLEM') WRITE(6,2000)((VAL(I,J),J=1,4),I=1,NPERM) WRITE(6,3000)(J,(VEC(I,J),I=1,6),J=1,NPERM) STOP END C C SUBROUTINE OP(N,M,P,Q) C C THIS IS THE MATRIX MULTIPLY FOR THE MATRIX A. DOUBLE PRECISION P(N,M),Q(N,M),A(100) COMMON /MULT/A DO 20 I=1,M DO 10 J=1,N Q(J,I)=A(J)*P(J,I) 10 CONTINUE 20 CONTINUE RETURN END C C SUBROUTINE OP1(N,M,P,Q) C C THIS IS THE MATRIX MULTIPLY FOR THE TRANSFORMED A. DOUBLE PRECISION P(N,M),Q(N,M),A(100) COMMON /MULT/A DO 20 I=1,M DO 10 J=1,N Q(J,I)=P(J,I)/(A(J)-5.0D0) 10 CONTINUE 20 CONTINUE RETURN END C C SUBROUTINE IOVECT(N,M,Q,J,K) C C THIS SUBROUTINE STORES AND RECALLS VECTORS. IT STORES THE VECTORS C IN CORE ALTHOUGH MOST REAL PROBLEMS WOULD USE A DISK. DOUBLE PRECISION Q(N,M),QS(110,50) COMMON /STORE/QS IF(K.EQ.1)GO TO 30 DO 20 L=1,M L1=J-M+L DO 10 I=1,N QS(I,L1)=Q(I,L) 10 CONTINUE 20 CONTINUE RETURN C 30 DO 50 L=1,M L1=J-M+L DO 40 I=1,N Q(I,L)=QS(I,L1) 40 CONTINUE 50 CONTINUE RETURN END -------