PROGRAM LANCZOS_PRINCIPLE6
C=======================================================================
C FORMATION OF ORTHONORMAL VECTOR AND TRI-DAIGONAL MATRIX BY
C LANCZOS PRINCIPAL ALONG WITH
C SHIFT-INVERT LANCZOS METHOD
C AND
C MATRIX DECOMPOSITION OF MATRIX [A]
C EIGENVALUES ARE SOLVED BY BISECTION METHOD
C ***** COMPUTATION OF [A]{U}1 BY [L][D][L]t SOLUTION *****
C 2011/FEB/9 EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXN=200, EPS=1.D-13, MXIGEN=200 )
C------- LANCZOS MEMORY
DIMENSION A(MXN,MXN), U(MXN,MXIGEN), ALPHA(MXIGEN), BETA(MXIGEN),
* R(MXIGEN), Q(MXIGEN,MXIGEN), Z(MXN)
C------- BISECTION METHOD MEMORY
DIMENSION P(MXIGEN), EIGENVEC(MXN,MXIGEN), EIGEN(MXIGEN)
DIMENSION TINVERSE(MXIGEN,MXIGEN), Q2(MXIGEN),
* T(MXIGEN,MXIGEN)
C=======================================================================
OPEN (1,FILE='LANCZOS-PRINCIPLE6.OUT', STATUS='UNKNOWN')
WRITE (1,*) '** SOLUTION OF SHIFTED-INVERSE LANCZOS METHOD **'
WRITE (1,*) '** INVERSE BY [L][D][L]t, SHIFT-PARAM = DELTA **'
WRITE (1,*) '******* m <= n ******* see NLANCZOS'
WRITE (1,*) 'EIGENVALUES BY BISECTION2'
C=======================================================================
CALL INPUT ( MXN, N, DELTA, NLANCZOS, A )
C=======================================================================
C--------- EVALUATION OF INITIAL VECTOR {U}1
CALL INITLVEC ( MXN,MXIGEN, N, A, U )
C=======================================================================
C--------- IMPLEMENTATION OF SHIFT PARAMETER(DELTA) INTO [A]
DO I = 1 , N
DO J = 1 , N
END DO
A(I,I) = A(I,I) - DELTA
END DO
C=======================================================================
C--------- DECOMPOSITION OF [A] INTO [L] AND [D]
CALL DECOMP ( MXN, N, A )
C=======================================================================
C--------- EVALUATION OF [U], ALPHA, AND BETA.
CALL LANCZOS ( MXIGEN,MXN,N,NLANCZOS,A,Z,R,U,ALPHA,BETA )
C=======================================================================
C--------- EVALUATION OF IGENVALUES AND VECTORS ---
CALL BSECTION2 ( EPS,MXIGEN,NLANCZOS,ALPHA,BETA,P,EIGEN )
C=======================================================================
WRITE (1,*) 'EIGENVALUES BY BISECTION METHOD USING ALPHA AND BETA'
WRITE (1,*) 'MODE SUBSPACE-EIGENVALUES'
DO I = 1 , NLANCZOS
WRITE (1,*) I, EIGEN(I)
END DO
CALL VECCOMP (EPS,MXIGEN,NLANCZOS,Q,EIGEN,TINVERSE,
* ALPHA,BETA, Q2, T )
WRITE (1,*)
WRITE (1,*) 'EIGEN VECTORS IN SUBSPACE'
DO I = 1 , NLANCZOS
WRITE (1,*) (Q(I,J), J= 1 , NLANCZOS)
END DO
WRITE (1,*)
C=======================================================================
C--------- EIGENVALUES IN REAL SPACE -------
WRITE (1,*) 'MODE REAL-SPACE-EIGENVALUES'
DO I = 1 , NLANCZOS
EIGEN(I) = DELTA + 1.D0 / EIGEN(I)
WRITE (1,*) I, EIGEN(I)
END DO
WRITE (1,*)
C=======================================================================
C---------- [EIGENVEC]=[U][Q] ------------
DO I = 1 , N
DO J = 1 , NLANCZOS
EIGENVEC(I,J) = 0.D0
DO K = 1 , NLANCZOS
EIGENVEC(I,J) = EIGENVEC(I,J) + U(I,K)*Q(K,J)
END DO
END DO
END DO
WRITE (1,*) 'EIGEN VECTORS IN REAL SPACE'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , NLANCZOS)
END DO
C=======================================================================
CLOSE (1)
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE DECOMP ( MXN, NNODE, A )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MXN,MXN)
C------- COMPUTATION OF UPPER L(I,J) --------
C-------I = 1
DO J = 2, NNODE
A(1,J) = A(1,J)/A(1,1)
END DO
C------- I >= 2
DO I = 2 , NNODE
SUM = 0.D0
DO M = 1, I-1
SUM = SUM + A(M,I)**2*A(M,M)
END DO
A(I,I) = A(I,I) - SUM
DO J = I+1, NNODE
SUM = 0.D0
DO M = 1, I-1
SUM = SUM + A(M,I)*A(M,J)*A(M,M)
END DO
A(I,J) = (A(I,J)-SUM)/A(I,I)
END DO
END DO
C ----------- MAKE LOWER [L] -----------
DO I = 1 , NNODE
DO J = I+1, NNODE
A(J,I) = A(I,J)
END DO
END DO
RETURN
END
C
C
SUBROUTINE SYSTEMDP ( MXN, NNODE, A, B )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MXN,MXN), B(MXN)
C======= STEP 1: OBTAIN {X'} BY [L]{X'}={B}
C------- NOTE THAT B(1) IS KNOWN
DO I = 2 , NNODE
SUM = 0.D0
DO J = 1, I-1
SUM = SUM + A(I,J)*B(J)
END DO
B(I) = B(I) - SUM
END DO
C======= STEP 2: OBTAIN {X''} BY [D]{X''}={X}
DO I = 1 , NNODE
B(I) = B(I) / A(I,I)
END DO
C======= STEP 3: OBTAIN {X} BY [L]T{X}={X''}
C------- NOTE THAT B(NNODE) IS KNOWN
DO I = NNODE-1, 1, -1
SUM = 0.D0
DO J = I+1 , NNODE
SUM = SUM + A(I,J)*B(J)
END DO
B(I) = B(I) - SUM
END DO
RETURN
END
C
C
SUBROUTINE LANCZOS ( MXIGEN,MXN,N,NLANCZOS,A,Z,R,U,ALPHA,BETA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXN), U(MXN,MXIGEN), ALPHA(MXIGEN), BETA(MXIGEN),
* R(MXIGEN), Z(MXN)
C
DO IGEN = 1 , NLANCZOS
DO I = 1 , N
Z(I) = U(I,IGEN)
END DO
C--------- COMPUTATION OF {Z} OF [A]{Z}={U}1 BY [L][D][L]T DECOMP. ----
CALL SYSTEMDP ( MXN, N, A, Z )
C--------- COMPUTATION OF ALPHA VALUE
ALPHA(IGEN) = 0.D0
DO I = 1 , N
ALPHA(IGEN) = ALPHA(IGEN) + Z(I)*U(I,IGEN)
END DO
IF ( IGEN .EQ. N ) EXIT
C--------- COMPUTATION OF {R}
DO I = 1 , N
R(I) = Z(I) - ALPHA(IGEN)*U(I,IGEN)
END DO
IF ( IGEN .GT. 1 ) THEN
DO I = 1 , N
R(I) = R(I) - BETA(IGEN-1)*U(I,IGEN-1)
END DO
END IF
C-------- RE-ORTHOGONALIZATION
DO I = 1 , IGEN
DOTPRDCT = 0.0D0
DO J = 1 , N
DOTPRDCT = DOTPRDCT + U(J,I)*R(J)
END DO
DO J = 1 , N
R(J) = R(J) - DOTPRDCT*U(J,I)
END DO
END DO
C-------- NEW NORMALIZED OTHOGONAL VECTOR
BETA(IGEN) = 0.D0
DO I = 1 , N
BETA(IGEN) = BETA(IGEN) + R(I)*R(I)
END DO
BETA(IGEN) = DSQRT(BETA(IGEN))
DO I = 1 , N
U(I,IGEN+1) = R(I)/BETA(IGEN)
END DO
C
END DO
C--------- PRINTS ALPHA AND BETA
WRITE (1,*) 'ALPHA BETA'
DO I = 1 , NLANCZOS - 1
WRITE (1,*) ALPHA(I), BETA(I)
END DO
WRITE (1,*) ALPHA(NLANCZOS)
WRITE (1,*)
RETURN
END
C
C
SUBROUTINE INPUT ( MXN, N, DELTA, NLANCZOS, A )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXN)
C------- MATRIX [A]
C------ INITIALIZATION
C--------- SIZE OF MATRIX [A]
N = 10
DO I = 1 , N
DO J = 1 , N
A(I,J) = N - I + 1
IF ( J .GT. I ) A(I,J) = N - J + 1
END DO
END DO
C-------- SHIFT PARAMETER, DELTA
DELTA = 0.2D0
C--------- NUMBER OF IGENVALUES
NLANCZOS = 10
C---------------- ECHO PRINT ------------------
WRITE (1,*) 'SIZE OF MATRIX [A] =', N
WRITE (1,*) 'MATRIX [A]'
DO I = 1 , N
WRITE (1,*) (A(I,J),J=1,N)
END DO
WRITE (1,*)
WRITE (1,*) 'SHIFT PARAMETER =', DELTA
WRITE (1,*)
WRITE (1,*) 'NUMBER OF EIGENVALUES =',NLANCZOS
WRITE (1,*)
RETURN
END
C
C
SUBROUTINE INITLVEC ( MXN,MXIGEN, N, A, U )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MXN,MXN), U(MXN,MXIGEN)
C-------- GENERATION OF INITIAL VECOTR {U}1 THE LENGTH = 1
VECL = 0.D0
DO I = 1 , N
VECL = VECL + A(I,I)*A(I,I)
END DO
VECL = DSQRT(VECL)
DO I = 1 , N
U(I,1) = A(I,I)/VECL
END DO
WRITE (1,*) 'INITIAL VECTOR {U}1'
DO I = 1 , N
WRITE (1,*) U(I,1)
END DO
WRITE (1,*)
RETURN
END
C
C
SUBROUTINE VECCOMP (EPS,MXIGEN,NLANCZOS,Q,EIGEN,TINVERSE,
* ALPHA,BETA, Q2,T )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ALPHA(MXIGEN), BETA(MXIGEN),T(MXIGEN,MXIGEN),
* EIGEN(MXIGEN),Q(MXIGEN,MXIGEN), TINVERSE(MXIGEN,MXIGEN),
* Q2(MXIGEN)
LOGICAL DCSN4
C------- THIS EVALUATES EIGENVECTORS IN SUB-SPACE BASED ON
C------- EIGENVALUES FOUND IN BISECTION METHOD.
C------- AVL = AVERAGE VECTOR LENGTH
AVL = DSQRT ( 1.D0/NLANCZOS )
C------------- COMPUTATION OF EIGENVECTORS --------------
C----------- INITIAL SETTING FOR [Q] AND {Q2} -----------
DO I=1,NLANCZOS
DO J=1,NLANCZOS
Q(I,J)=0.0D0
END DO
Q(I,I) = 1.D0
END DO
C++++++++++++++++++++++++++++
DO MODE = 1 , NLANCZOS
C------- COMPUTATION OF INVERSE MATRIX OF [[T]-LAMDA(1+EPS)[I]] --------
DO I = 1 , NLANCZOS
DO J = 1 , NLANCZOS
T(I,J) = 0.D0
END DO
T(I,I) = ALPHA(I) - EIGEN(MODE)*(1.0D0+EPS)
END DO
DO I = 1 , NLANCZOS-1
T(I,I+1) = BETA(I)
T(I+1,I) = BETA(I)
END DO
C----------- DECOMPOSITION OF [T]
CALL DECOMP ( MXIGEN, NLANCZOS, T )
C----------- COMPUTATION OF [T] INVERSE
C::::::::::::::::::::::::::::::::::::::::::::::::
DO I = 1 , NLANCZOS
DO J = 1 , NLANCZOS
Q2(J) = 0.D0
END DO
Q2(I) = 1.D0
CALL SYSTEMDP ( MXIGEN, NLANCZOS, T, Q2 )
DO J = 1 , NLANCZOS
TINVERSE(J,I) = Q2(J)
END DO
END DO
C::::::::::::::::::::::::::::::::::::::::::::::::
C--------- COMPUTATION OF EIGENVECTORS BY POWER METHOD ----------
DCSN4 = .TRUE.
C*****************************************
DO WHILE ( DCSN4 )
DIFFMAX = 0.D0
C------------- COMPUTATION OF RIGHT HAND SIDE {Q2}---------
DO I = 1 , NLANCZOS
Q2(I)=0.D0
DO K = 1,NLANCZOS
Q2(I) = Q2(I) + TINVERSE(I,K)*Q(K,MODE)
END DO
END DO
C------------- COMPUTATION OF VECTOR LENGTH --------------
VECLNGH = 0.D0
DO I = 1 , NLANCZOS
VECLNGH = VECLNGH + Q2(I)*Q2(I)
END DO
VECLNGH = DSQRT(VECLNGH)
C-------------- IMPROVED EIGENVECTOR AND ERROR COMPUTATION -------------
DO I=1,NLANCZOS
DIFFATI = ( DABS(Q(I,MODE))-DABS(Q2(I)/VECLNGH) ) / AVL
DIFFMAX = DMAX1 ( DIFFMAX, DIFFATI )
Q(I,MODE) = Q2(I)/VECLNGH
END DO
DCSN4 = ( DIFFMAX .GT. EPS )
C------------- BELOW END OF DO WHILE ----------
END DO
C*****************************************
END DO
C++++++++++++++++++++++++++++
RETURN
END
C
C
SUBROUTINE BSECTION2 ( EPS,MXENGN,NEIGEN,ALPHA,BETA,P,EIGEN )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(MXENGN), EIGEN(MXENGN),ALPHA(MXENGN),BETA(MXENGN)
LOGICAL DCSN2, DCSN5
MAXITA = 1000
C-------- COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST
BANDMAX= DABS(ALPHA(1))+BETA(1)
DO I=2,NEIGEN-1
BANDMAX = DMAX1( BANDMAX, DABS(ALPHA(I))+BETA(I)+BETA(I-1) )
END DO
BANDMAX= DMAX1(BANDMAX, DABS(ALPHA(NEIGEN))+BETA(NEIGEN-1) )
DO I = 1 , NEIGEN
ALPHA (I) = ALPHA(I)/BANDMAX
BETA (I) = BETA (I)/BANDMAX
END DO
C-------- COMPUTATION OF IGENVALUES ----------
DO MODE = 1, NEIGEN
ITA = 0
BOTTOLMT =-1.D0
UPPERLMT = 1.D0
FMAG = 1.D0
C------------ INITIAL SETTING FOR CONVERGENCE PARAMETERS ---------------
DCSN2 = .TRUE.
DCSN5 = .TRUE.
C*** BEGIN OF DO WHILE ***
DO WHILE ( DCSN2 .AND. DCSN5 )
DCSN2 = DABS(UPPERLMT-BOTTOLMT)/FMAG .GT. EPS
ITA = ITA + 1
DCSN5 = ITA .LT. MAXITA
FLAMBDA = 0.5D0*(BOTTOLMT+UPPERLMT)
C------- FORMING STURM SEQUENCES
C------- P(0), P(1), P(2), ......., P(NEIGEN). NOTE: P(0)=1.
P(1) = ALPHA(1)-FLAMBDA
P(2) = (ALPHA(2)-FLAMBDA)*P(1) - BETA(1)**2
DO J = 3 , NEIGEN
P(J) = (ALPHA(J)-FLAMBDA)*P(J-1) - BETA(J-1)**2*P(J-2)
END DO
C------- COUNT NUMBER OF SIGN-CHANGES IN P(I)
KOUNT = 0
PRODUCT = 1.D0
DO J = 1 , NEIGEN
IF ( PRODUCT*P(J) .LT. 0.D0 ) THEN
KOUNT = KOUNT + 1
PRODUCT = - PRODUCT
END IF
END DO
C------- NARROW DOWN LIMTS WHERE EIGENVALUE OF MODE EXISTS
IF ( KOUNT .LT. MODE ) THEN
BOTTOLMT = FLAMBDA
ELSE
UPPERLMT = FLAMBDA
END IF
FMAG = DMAX1 (DABS(UPPERLMT), DABS(BOTTOLMT))
IF ( FMAG .EQ. 0.D0 ) FMAG = 1.D0
END DO
C*** END OF DO WHILE ***
EIGEN(MODE) = FLAMBDA*BANDMAX
END DO
DO I = 1 , NEIGEN
ALPHA (I) = ALPHA(I)*BANDMAX
BETA (I) = BETA (I)*BANDMAX
END DO
RETURN
END