PROGRAM LANCZOS_PRINCIPLE7_WITH_NEWTON
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 ****** NEWTON-RAPHSON METHOD IMPLEMENTED IN BISECTION ALGORITHM ******
C************* CRUSHED [T] AND DECOMPOSITION METHOD USED ***************
C 2011/FEB/9 EIJI FUKUMORI
C [A] = MATRIX [A]
C [Q] = EIGENVECTORS IN SUBSPACE
C [V] = EIGENVECTORS IN REALSPACE
C [U] = LANCZOS ORTHOGONAL VECTORS
C [T] = LANCZOS [T]
C {P} = STURM COLUMNS
C {R} = {r} IN LANCZOS METHOD
C {EIGEN} = EIGENVALUES IN SUB-SPACE ALSO IN REAL SPACE
C EQUATION TO BE SOLVED: [T]{new_Q2}={old_Q2}
C------------- DECOMPOSITION CRUSHED [T] INTO [L][D][L]t ---------------
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXN=200, EPS=1.D-14, MXIGEN=200 )
C-------------------------- LANCZOS MEMORY -----------------------------
DIMENSION A(MXN,MXN), U(MXN,MXIGEN), ALPHA(MXIGEN), BETA(MXIGEN),
* R(MXIGEN), Z(MXN), V(MXN,MXIGEN)
C---------------------- BISECTION METHOD MEMORY ------------------------
DIMENSION P(MXIGEN), EIGEN(MXIGEN)
CHARACTER*8 MODENAME(MXIGEN)
C------------------ INVERSE POWER METHOD MEMORY ------------------------
DIMENSION T(MXIGEN,2), Q(MXIGEN,MXIGEN),Q2(MXIGEN)
C=======================================================================
OPEN (1,FILE='LANCZOS-PRINCIPLE7.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 IN INPUT'
WRITE (1,*) 'EIGENVALUES BY BISECTION2/NEWTON-RAPHSON'
C=======================================================================
CALL INPUT ( MXN, N, DELTA, NLANCZOS, A )
C=======================================================================
C---------------- EVALUATION OF INITIAL VECTOR {U}1 --------------------
CALL INITLVEC ( MXN,MXIGEN, N, A, U )
WRITE (1,*) 'INITIAL VECTOR {U}1'
DO I = 1 , N
WRITE (1,*) U(I,1)
END DO
WRITE (1,*)
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--------- PRINTS ALPHA AND BETA
WRITE (1,*) 'LANCZOS {ALPHA} AND {BETA}'
DO I = 1 , NLANCZOS - 1
WRITE (1,*) ALPHA(I), BETA(I)
END DO
WRITE (1,*) ALPHA(NLANCZOS)
WRITE (1,*)
C=======================================================================
C--------------------- EVALUATION OF IGENVALUES ------------------------
CALL BSECTION3 ( EPS,MXIGEN,NLANCZOS,ALPHA,BETA,P,EIGEN )
WRITE (1,*) 'EIGENVALUES BY BISECTION METHOD USING ALPHA AND BETA'
WRITE (1,*) 'MODE SUBSPACE-EIGENVALUES {EIGEN}'
DO I = 1 , NLANCZOS
WRITE (1,*) I, EIGEN(I)
END DO
C=======================================================================
C-- COMPUTATION OF EIGENVECTORS IN SUB-SPACE BY INVERSE POWER METHOD --
CALL VECCOMP (EPS,MXIGEN,NLANCZOS,Q,EIGEN,ALPHA,BETA, Q2, T )
WRITE (1,*)
WRITE (1,*) 'EIGEN VECTORS IN SUBSPACE [Q]'
DO I = 1 , NLANCZOS
WRITE (1,*) (Q(I,J), J= 1 , NLANCZOS)
END DO
C=======================================================================
WRITE (1,*)
WRITE (1,*) 'CONFIRMATION OF SUB-SPACE EI-VECTOR ORTHOGONALITY'
DOTMAX = 0.D0
DO MODEI = 1 , NLANCZOS-1
DO MODEJ = MODEI+1 , NLANCZOS
SUM = 0.D0
DO K = 1 , NLANCZOS
SUM = SUM + Q(K,MODEI)*Q(K,MODEJ)
END DO
DOTMAX = DMAX1 ( DOTMAX, ABS(SUM) )
END DO
END DO
WRITE (1,*) 'MAXIMUM ERROR IN DOT PRODUCT =', DOTMAX
C=======================================================================
C------------- COMPUTATION OF EIGENVALUESS IN REAL-SPACE ---------------
WRITE (1,*)
C--------- EIGENVALUES IN REAL SPACE -------
WRITE (1,*) 'MODE REAL-SPACE-EIGENVALUES {EIGEN}'
DO I = 1 , NLANCZOS
EIGEN(I) = DELTA + 1.D0 / EIGEN(I)
WRITE (1,*) I, EIGEN(I)
END DO
WRITE (1,*)
C=======================================================================
C------------- COMPUTATION OF EIGENVECTORS IN REAL-SPACE ---------------
C----------------------------- [V]=[U][Q] ------------------------------
DO I = 1 , N
DO J = 1 , NLANCZOS
V(I,J) = 0.D0
DO K = 1 , NLANCZOS
V(I,J) = V(I,J) + U(I,K)*Q(K,J)
END DO
END DO
END DO
WRITE (1,*) 'EIGEN VECTORS IN REAL SPACE [V]'
DO I = 1 , N
WRITE (1,*) (V(I,J), J= 1 , NLANCZOS)
END DO
C=======================================================================
WRITE (1,*)
WRITE (1,*) 'CONFIRMATION OF REAL-SPACE EI-VECTOR ORTHOGONALITY'
DOTMAX = 0.D0
DO MODEI = 1 , NLANCZOS-1
DO MODEJ = MODEI+1 , NLANCZOS
SUM = 0.D0
DO K = 1 , N
SUM = SUM + Q(K,MODEI)*Q(K,MODEJ)
END DO
DOTMAX = DMAX1 ( DOTMAX, ABS(SUM) )
END DO
END DO
WRITE (1,*) 'MAXIMUM ERROR IN DOT PRODUCT =', DOTMAX
C=======================================================================
CLOSE (1)
STOP 'NORMAL TERMINATION'
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
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 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]'
IF ( N .LE. 20 ) THEN
DO I = 1 , N
WRITE (1,*) (A(I,J),J=1,N)
END DO
END IF
WRITE (1,*)
WRITE (1,*) 'SHIFT PARAMETER =', DELTA
WRITE (1,*)
WRITE (1,*) 'NUMBER OF EIGENVALUES TO BE EVALUATED=',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
RETURN
END
C
C
SUBROUTINE VECCOMP ( EPS,MXIGEN,NLANCZOS,Q,EIGEN,ALPHA,BETA,Q2,T )
C- THIS EVALUATES EIGENVECTORS IN SUB-SPACE BASED ON KNOWN EIGENVALUES -
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ALPHA(MXIGEN), BETA(MXIGEN),EIGEN(MXIGEN)
C----------------- ARRAY FOR EIGENVECTOR EVALUATION --------------------
DIMENSION T(MXIGEN,2),Q(MXIGEN,MXIGEN), Q2(MXIGEN)
LOGICAL DCSN4
C-------- [Q] = EIGENVECTORS
C-------- [T] = LANCZOS TRIDIAGONAL MATRIX (CRUSHED ARRAY)
C-------- {Q2} = DUMMY VECTOR FOR EIGENVECTOR
C-------- {EIGEN} = EIGENVALUES
C-------- AVL = AVERAGE VECTOR LENGTH
AVL = DSQRT ( 1.D0/NLANCZOS )
C----------------------- INITIAL VALUES FOR [Q] ------------------------
DO I=1,NLANCZOS
DO J=1,NLANCZOS
Q(I,J)=0.0D0
END DO
Q(I,I) = 1.D0
END DO
C---------------- EIGENVECTOR COMPUTATION STARTS HERE ------------------
DO MODE = 1 , NLANCZOS
C------- COMPUTATION OF INVERSE MATRIX OF [[T]-LAMDA(1+EPS)[I]] --------
BETA(NLANCZOS) = 0.D0
DO I = 1 , NLANCZOS
T(I,1) = ALPHA(I) - EIGEN(MODE)*(1.D0+EPS)
T(I,2) = BETA(I)
END DO
C------------ DECOMPOSITION OF MATRIX [T] INTO [L][D][L]t --------------
CALL LLDECOMPVT ( T, NLANCZOS, 2, 2, MXIGEN )
C-------------------- INITIALIZATION OF RHS {Q2} -----------------------
DO J = 1 , NLANCZOS
Q2(J) = Q(J,MODE)
END DO
C-------- COMPUTATION OF EIGENVECTORS BY INVERSE POWER METHOD ----------
DCSN4 = .TRUE.
DO WHILE ( DCSN4 )
DIFFMAX = 0.D0
C---------------------- EVALUATION OF NEW {Q2} -------------------------
CALL SYSLDLVT ( MXIGEN, 2, NLANCZOS, 2, T, Q2 )
C----------------- VECTOR LENGTH COMPUTATION OF {Q2} -------------------
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
Q2(I) = Q(I,MODE)
END DO
DCSN4 = ( DIFFMAX .GT. EPS )
C------------- END OF DO WHILE ----------
END DO
C------------- END OF DO MODE -----------
END DO
RETURN
END
C
C
SUBROUTINE LLDECOMPVT ( A, NNODE, IBAND, MXW, MXN )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MXN,MXW)
C-------I = 1
DO J = 2, IBAND
A(1,J) = A(1,J)/A(1,1)
END DO
C------ I >= 2
DO I = 2 , NNODE
SUM = 0.D0
DO M = 2, MIN0(I,IBAND)
II = I - M + 1
SUM = SUM + A(II,M)**2*A(II,1)
END DO
A(I,1) = A(I,1) - SUM
DO J = 2, MIN0(IBAND,NNODE-I+1)
SUM = 0.D0
DO M = 2, MIN0( I,IBAND, IBAND-J+1 )
II = I - M + 1
L = J + M - 1
SUM = SUM + A(II,1)*A(II,M)*A(II,L)
END DO
A(I,J) = (A(I,J)-SUM)/A(I,1)
END DO
END DO
RETURN
END
C
C
SUBROUTINE SYSLDLVT ( MXN, MXW, NNODE, NBWDTH, GSMTX, V1 )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION GSMTX(MXN,MXW), V1(MXN)
C------- SOLUTION OF [A]{X}={B} --------
C----- SEQUENCE 1: [L]{X'}={B}
DO I=2,NNODE
DO J = 2, MIN0(I,NBWDTH)
V1(I) = V1(I) - GSMTX(I-J+1,J)*V1(I-J+1)
END DO
END DO
C------ SEQUENCE 2 [D]{X"'}={X"}
DO I = 1 , NNODE
V1(I) = V1(I) / GSMTX(I,1)
END DO
C------ SEQUENCE 3 [L]TRANSPOSE{X'}={X"'} NOTE THAT L(I,I) = 1
DO I = NNODE-1 , 1 , -1
DO J = 2 , MIN0(NNODE-I+1,NBWDTH)
V1(I) = V1(I) - GSMTX(I,J)*V1(I+J-1)
END DO
END DO
RETURN
END
C
C
SUBROUTINE BSECTION3 ( EPS,MXIGEN,NLANCZOS,ALPHA,BETA,P,EIGEN )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ALPHA(MXIGEN), BETA(MXIGEN), P(MXIGEN), EIGEN(MXIGEN)
DIMENSION RESULT(100,100)
LOGICAL DCSN1, DCSN2, DCSN3, DCSN4, DCSN5
MAXITA = 1000
C------------ COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST ------------
BETA(NLANCZOS) = 0.D0
BANDMAX= DABS(ALPHA(1))+BETA(1)
DO I=2,NLANCZOS-1
BANDMAX = DMAX1( BANDMAX, DABS(ALPHA(I))+BETA(I)+BETA(I-1) )
END DO
BANDMAX = DMAX1(BANDMAX, DABS(ALPHA(NLANCZOS))+BETA(NLANCZOS-1) )
C------------ NORMALIZATION OF ALPHA AND BETA BY BANDMAX ---------------
DO I = 1 , NLANCZOS
ALPHA(I) = ALPHA(I)/BANDMAX
BETA (I) = BETA (I)/BANDMAX
END DO
C-------------- ERROR CHECK VALUE WITHIN BISECTION METHOD --------------
EPSBI =EPS**0.1
C-------------------- COMPUTATION OF EIGENVALUES -----------------------
ITAMAX = 0
DO MODE = 1, NLANCZOS
BOTTOLMT =-BANDMAX/BANDMAX
UPPERLMT = BANDMAX/BANDMAX
ITA = 0
C------------ INITIAL SETTING FOR CONVERGENCE PARAMETERS ---------------
N1 = 0
N2 = NLANCZOS
BASE = UPPERLMT
DCSN1 = .FALSE.
DCSN2 = .TRUE.
DCSN3 = .FALSE.
DCSN4 = .FALSE.
DCSN5 = .TRUE.
C------- BEGIN OF DO WHILE -----
DO WHILE ( DCSN2 .AND. DCSN5 )
ITA = ITA + 1
ITAMAX = MAX0(ITAMAX , ITA)
DCSN5 = ITA .LT. MAXITA
C------------ PREDICTION OF NEW LAMBDA BY BISECTION METHOD -------------
FLAMBDA = 0.5D0*( BOTTOLMT + UPPERLMT )
C---------- PREDICTION OF NEW LAMBDA BY NEWTON-RAPHSON METHOD ----------
IF ( DCSN4 ) THEN
IF ( DCSN1 .AND. DCSN3 ) FLAMBDA = FLNEWTON
END IF
C----------------------- FORMING STURM COLUMNS ------------------------
P(1) = ALPHA(1)-FLAMBDA
P(2) = (ALPHA(2)-FLAMBDA)*P(1) - BETA(1)**2
DO J = 3 , NLANCZOS
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 , NLANCZOS
IF ( PRODUCT*P(J) .LT. 0.D0 ) THEN
KOUNT = KOUNT + 1
PRODUCT = - PRODUCT
END IF
END DO
C--------- NARROW DOWN LIMITS WHERE MODEth EIGENVALUE EXISTS -----------
IF ( KOUNT .LT. MODE ) THEN
BOTTOLMT = FLAMBDA
N1 = KOUNT
ELSE
UPPERLMT = FLAMBDA
N2 = KOUNT
END IF
C----------------------------- MAGNITUDE -------------------------------
FMAG = DMAX1 (DABS(UPPERLMT), DABS(BOTTOLMT))
C--------------------- CHECK LOCATION OF LIMITS ------------------------
IF ( N2 .EQ. N1 ) EXIT
IF ( FLAMBDA .EQ. BASE ) EXIT
C-------------------- CONVERGENCE CHECK: BISECTION ---------------------
DCSN1 = ( N2-N1 .EQ. 1 )
IF ( DCSN1 ) THEN
DCSN2 = DABS(UPPERLMT-BOTTOLMT)/FMAG .GT. EPS
DCSN3 = .FALSE.
DCSN4 = DABS(UPPERLMT-BOTTOLMT)/FMAG .LE. EPSBI
C--------------- PREDICTION BY NEWTON-RAPHSON METHOD -------------------
DPDLAMBDA = (P(NLANCZOS)-PN)/(FLAMBDA-BASE)
IF ( DABS(DPDLAMBDA) .NE. 0.D0 ) THEN
DLAMBDA = - P(NLANCZOS) / DPDLAMBDA
FLNEWTON = FLAMBDA + DLAMBDA
DCSN3 = (FLNEWTON-BOTTOLMT)*(FLNEWTON-UPPERLMT) .LE. 0.D0
END IF
C----------------------- CONVERGENCE CHECK -----------------------------
IF ( DCSN3 ) DCSN2 = ABS(DLAMBDA)/FMAG .GT. EPS
END IF
C------------------ ADVANCEMENT OF PN AND LAMBDA -----------------------
PN = P(NLANCZOS)
BASE = FLAMBDA
C----------------- END OF DO WHILE --------------------
END DO
EIGEN(MODE) = FLAMBDA*BANDMAX
IF ( DCSN3 ) EIGEN(MODE) = FLNEWTON*BANDMAX
C----------------- END OF DO MODE --------------------
END DO
DO I = 1 , NLANCZOS
ALPHA(I) = ALPHA(I)*BANDMAX
BETA (I) = BETA (I)*BANDMAX
END DO
RETURN
END