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
-------