PROGRAM LANCZOS_PRINCIPLE5
C=======================================================================
C FORMATION OF ORTHONORMAL VECTOR AND TRI-DAIGONAL MATRIX BY
C LANCZOS PRINCIPAL ALONG WITH
C SHIFT-INVERT LANCZOS METHOD
C EIGENVALUES ARE SOLVED BY BISECTION METHOD
C***** COMPUTATION OF [A]{U}1 BY ORDINALLY SIMULTANEOUS SOLUTION *****
C I.E. [A]{X}={Z}
C 2011/FEB/4 EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXN=50, EPS=1.D-15, MXIGEN=50 )
C------- LANCZOS MEMORY
DIMENSION A(MXN,MXN), U(MXN,MXIGEN), ALPHA(MXIGEN), BETA(MXIGEN),
* R(MXIGEN), Q(MXIGEN,MXIGEN), ABAR(MXN,MXN),
* ABARSAVE(MXN,MXN), Z(MXN)
C------- BISECTION METHOD MEMORY
DIMENSION F(MXIGEN), EIGENVEC(MXN,MXIGEN), EIGEN(MXIGEN)
DIMENSION TINVERSE(MXIGEN,MXIGEN), Q2(MXIGEN),
* T(MXIGEN,MXIGEN)
C=======================================================================
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=======================================================================
OPEN (1,FILE='LANCZOS-PRINCIPLE5.OUT', STATUS='UNKNOWN')
WRITE (1,*) '** SOLUTION OF SHIFTED-INVERSE LANCZOS METHOD **'
WRITE (1,*) '** INVERSE BY [A]{X}={Z}, SHIFT-PARAM = DELTA **'
WRITE (1,*) '******* m <= n ******* see NLANCZOS'
WRITE (1,*) 'EIGENVALUES BY BISECTION2'
C=======================================================================
C--------- GIVEN [A]
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
C=======================================================================
C-------- SHIFT PARAMETER, DELTA
DELTA = 0.25D0
WRITE (1,*) 'SHIFT PARAMETER =', DELTA
C=======================================================================
C-------- [A-BAR]
DO I = 1 , N
DO J = 1 , N
ABAR(I,J) = A(I,J)
END DO
ABAR(I,I) = ABAR(I,I) - DELTA
END DO
DO I = 1 , N
DO J = 1 , N
ABARSAVE(I,J) = ABAR(I,J)
END DO
END DO
C=======================================================================
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
C=======================================================================
C--------- NUMBER OF IGENVALUES
NLANCZOS = 10
WRITE (1,*) 'NUMBER OF IGENVALUES =',NLANCZOS
C=======================================================================
DO IGEN = 1 , NLANCZOS
C
DO I = 1 , N
DO J = 1 , N
ABAR(I,J) = ABARSAVE(I,J)
END DO
END DO
C
DO I = 1 , N
Z(I) = U(I,IGEN)
END DO
CALL SYSTEM ( MXN, N , ABAR , Z )
C
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
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
WRITE (1,*) 'ALPHA BETA'
DO I = 1 , NLANCZOS - 1
WRITE (1,*) ALPHA(I), BETA(I)
END DO
WRITE (1,*) ALPHA(NLANCZOS)
C=======================================================================
C--- EVALUATION OF IGENVALUES AND VECTORS BY BISECTION METHOD ---
CALL BSECTION2 ( EPS,MXIGEN,NLANCZOS,ALPHA,BETA,F,EIGEN )
WRITE (1,*)
WRITE (1,*) 'BISECTION METHOD WITH [T]'
WRITE (1,*) 'MODE SUBSPACE-EIGENVALUES'
DO I = 1 , NLANCZOS
WRITE (1,*) I, EIGEN(I)
END DO
C=======================================================================
CALL VECCOMP (EPS,MXIGEN,NLANCZOS,Q,EIGEN,TINVERSE,
* ALPHA,BETA, Q2, T )
WRITE (1,*) 'EIGEN VECTORS IN SUBSPACE'
DO I = 1 , NLANCZOS
WRITE (1,*) (Q(I,J), J= 1 , NLANCZOS)
END DO
WRITE (1,*) 'MODE REAL-SPACE-EIGENVALUES'
DO I = 1 , NLANCZOS
EIGEN(I) = DELTA + 1.D0 / EIGEN(I)
WRITE (1,*) I, EIGEN(I)
END DO
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 SYSTEM ( MXN, N , A , C )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A (MXN,MXN) , C (MXN)
N1 = N - 1
DO K = 1, N1
L = K + 1
DO I = L , N
P = A (I,K) / A (K,K)
IF ( P .NE. 0. ) THEN
DO J = L , N
A (I,J) = A (I,J) - P * A ( K , J )
END DO
C ( I ) = C ( I) - P * C ( K )
END IF
END DO
END DO
C---- BACK SUBSTITUTION
C (N) = C (N) / A (N,N)
DO K = 1 , N1
I = N - K
L = I + 1
P = C ( I )
DO J = L , N
P = P - C (J) * A (I,J)
END DO
C ( I ) = P / A (I,I)
END DO
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 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 BSECTION2 ( EPS,MXENGN,NEIGEN,ALPHA,BETA,F,EIGEN )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION F(MXENGN), EIGEN(MXENGN),
* U(MXENGN,MXENGN),ALPHA(MXENGN), BETA(MXENGN)
C-------- COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST
BMAX= DABS(ALPHA(1))+BETA(1)
DO I=2,NEIGEN-1
BMAX = DMAX1( BMAX, DABS(ALPHA(I))+BETA(I)+BETA(I-1) )
END DO
BMAX= DMAX1(BMAX, DABS(ALPHA(NEIGEN))+BETA(NEIGEN-1) )
C=========================================================================
DO MODE = 1, NEIGEN
A1 = 0.D0
A2 = BMAX
DENOMI = DABS ( DMAX1 ( A1, A2 ) )
IF ( DENOMI .EQ. 0.D0 ) DENOMI = BMAX
C--- BEGIN OF DO WHILE -----
DO WHILE ( DABS(A2-A1)/DENOMI .GT. EPS )
A3 = 0.5D0*(A1+A2)
C------- FORMING STURM SEQUENCES
C------- P-0, P-1, P-2, ......., P-NEIGEN. NOTE: P(0)=1.
F(1) = ALPHA(1)-A3
F(2) = (ALPHA(2)-A3)*F(1) - BETA(1)**2
DO J = 3 , NEIGEN
F(J) = (ALPHA(J)-A3)*F(J-1) - BETA(J-1)**2*F(J-2)
END DO
C-------- COUNT NUMBER OF SIGN-CHANGES IN F(I)
NA3 = 0
P = 1.D0
DO J = 1 , NEIGEN
IF ( P*F(J) .LT. 0.D0 ) THEN
NA3 = NA3 + 1
P = -P
END IF
END DO
C-------- NARROW DOWN LIMTS WHERE EIGENVALUE OF MODE EXISTS
IF ( NA3 .LT. MODE ) THEN
A1=A3
ELSE
A2=A3
END IF
DENOMI = DABS ( DMAX1 ( A1, A2 ) )
IF ( DENOMI .EQ. 0.D0 ) DENOMI = BMAX
END DO
C--- END OF DO WHILE -----
EIGEN(MODE)=A3
END DO
RETURN
END