PROGRAM LANCZOS_PRINCIPLE2
C=======================================================================
C FORMATION OF ORTHONORMAL VECTOR AND TRI-DAIGONAL MATRIX BY
C LANCZOS PRINCIPAL
C EIGENVALUES ARE SOLVED BY THE FOLLOWING METHODS:
C (1) BISECTION METHOD
C (2) JACOBS METHOD
C (3) POWER METHOD
C NOTE: *******BISECTION METHOD EMPLOYS [T] MATRIX*********
C 2011/JAN/25 EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXN=50, EPS=1.D-13 )
C------- LANCZOS MEMORY
DIMENSION A(MXN,MXN), U(MXN,MXN), ALPHA(MXN), BETA(MXN),
* AU1(MXN), R(MXN), AU(MXN,MXN), UTAU(MXN,MXN), T(MXN,MXN)
C------- POWER METHOD MEMORY
DIMENSION X(MXN),Y(MXN)
C------- BISECTION METHOD MEMORY
DIMENSION F(MXN)
C------- COMMON
DIMENSION ASAVE(MXN,MXN),TSAVE(MXN,MXN),EIGENVEC(MXN,MXN),
* EIGEN(MXN), Q2(MXN), TINVERSE(MXN,MXN)
C=======================================================================
OPEN (1,FILE='LANCZOS-PRINCIPLE2.OUT', STATUS='UNKNOWN')
WRITE (1,*) '** SOLUTION OF LANCZOS METHOD **'
WRITE (1,*) '******* m = n *******'
WRITE (1,*) 'EIGENVALUES BY BISECTION1, JACOBS, POWER METHODS'
C=======================================================================
C------ INITIALIZATION
C--------- SIZE OF MATRIX [A]
N = 10
C--------- GIVEN [A]
DO I = 1 , N
DO J = 1 , N
K = N-J+1
IF ( I .GT. J ) K = N-I+1
A(I,J) = K
END DO
END DO
WRITE (1,*)'MATRIX [A]'
DO I = 1 , N
WRITE (1,*) (A(I,J),J=1,N)
END DO
C-------- INITIAL VECTOR {U}1 THE LENGTH = 1
ABSU1 = 0.D0
DO I = 1 , N
ABSU1 = + ABSU1 + A(I,I)**2
END DO
ABSU1 = DSQRT (ABSU1)
DO I = 1 , N
U(I,1) = A(I,I)/ABSU1
END DO
WRITE (1,*) 'INITIAL VECTOR {U}1'
WRITE (1,*) (U(I,1),I=1,N)
C******************** STARTS LANCZOS_PRINCIPLE HERE *****
DO IGEN = 1 , N
C
DO I = 1 , N
AU1(I) = 0.D0
DO J = 1 , N
AU1(I) = AU1(I) + A(I,J)*U(J,IGEN)
END DO
END DO
C
ALPHA(IGEN) = 0.D0
DO I = 1 , N
ALPHA(IGEN) = ALPHA(IGEN) + AU1(I)*U(I,IGEN)
END DO
IF ( IGEN .EQ. N ) EXIT
C
DO I = 1 , N
R(I) = AU1(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******************** END OF LANCZOS_PRINCIPLE *********************
WRITE (1,*) 'RESULT OF [U] BY LANCZOS PRINCIPAL'
DO I = 1 , N
WRITE (1,*) ( U(I,J), J = 1 , N )
END DO
C------------------ COMPUTATION OF [A]*[U] = [AU]
DO I = 1 , N
DO J = 1 , N
AU(I,J) = 0.D0
DO K = 1 , N
AU(I,J) = AU(I,J) + A(I,K)*U(K,J)
END DO
END DO
END DO
C------- COMPUTATION OF [U]TANSPOSE[AU] = [U]TANSPOSE[A][U]
DO I = 1 , N
DO J = 1 , N
UTAU(I,J) = 0.D0
DO K = 1 , N
UTAU(I,J) = UTAU(I,J) + U(K,I)*AU(K,J)
END DO
END DO
END DO
C------------- PRINTING [U]TANSPOSE[A][U]
WRITE (1,*) 'RESULT OF [U]TANSPOSE[A][U]'
DO I = 1 , N
WRITE (1,*) ( UTAU(I,J), J = 1 , N )
END DO
C******** FORMATION OF [T] BY MAKING USE OF -------------------------
C---------------------- ALPHAS AND BETAS; I.E.LANCZOS METHOD ********
DO I = 1 , N
DO J = 1 , N
T(I,J) = 0.D0
END DO
T(I,I) = ALPHA(I)
END DO
DO I = 1 , N-1
T(I,I+1) = BETA(I)
T(I+1,I) = BETA(I)
END DO
C------- SAVE
DO I = 1 , N
DO J = 1 , N
TSAVE(I,J) = T(I,J)
ASAVE(I,J) = A(I,J)
END DO
END DO
C-------- PRINTING [T]
WRITE (1,*) 'DIRECT FORMATION OF [T] MATRIX BY LANCZOS METHOD'
DO I = 1 , N
WRITE (1,*) ( T(I,J), J = 1 , N )
END DO
C=======================================================================
C COMPUTATION OF IGENVALUES AND IGEN-VECTORS BY VARIOUS METHODS
C=======================================================================
C------- BISECTION METHOD
CALL BSECTION1 ( EPS, MXN, N, T, F, EIGEN )
WRITE (1,*)
WRITE (1,*) 'BISECTION METHOD WITH [T]'
WRITE (1,*) 'MODE EIGENVALUES'
DO I = 1 , N
WRITE (1,*) I, EIGEN(I)
END DO
C=======================================================================
CALL VECCOMP (EPS,MXN,N,EIGENVEC,EIGEN,TINVERSE,ALPHA,BETA,Q2,T )
WRITE (1,*) 'EIGEN VECTORS'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , N)
END DO
C=======================================================================
C------- JACOBS METHOD
DO I = 1 , N
DO J = 1 , N
T(I,J) = TSAVE(I,J)
END DO
END DO
MAXTRTN = 20000
C------------------- COMPUTATION WITH [T]
CALL JACOBS ( MXN, EPS, MAXTRTN, N, T, EIGENVEC, EIGEN )
WRITE (1,*)
WRITE (1,*) 'JACOBS METHOD WITH [T]'
WRITE (1,*) 'MODE EIGENVALUES'
DO I = 1 , N
WRITE (1,*) I, EIGEN(I)
END DO
WRITE (1,*) 'EIGEN VECTORS'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , N)
END DO
DO I = 1 , N
DO J = 1 , N
T(I,J) = TSAVE(I,J)
END DO
END DO
C------------------- COMPUTATION WITH [A]
CALL JACOBS ( MXN, EPS, MAXTRTN, N, A, EIGENVEC, EIGEN )
WRITE (1,*)
WRITE (1,*) 'JACOBS METHOD WITH [A]'
WRITE (1,*) 'MODE EIGENVALUES'
DO I = 1 , N
WRITE (1,*) I, EIGEN(I)
END DO
WRITE (1,*) 'EIGEN VECTORS'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , N)
END DO
DO I = 1 , N
DO J = 1 , N
A(I,J) = ASAVE(I,J)
END DO
END DO
C=======================================================================
C------- POWER METHOD
C------------------- COMPUTATION WITH [A]
CALL POWERMTD ( EPS, MXN, N, A, X, Y, EIGENVEC, EIGEN )
WRITE (1,*)
WRITE (1,*) 'POWER METHOD WITH [A]'
WRITE (1,*) 'MODE EIGENVALUE'
DO I = 1 , N
WRITE (1,*) I , EIGEN(I)
END DO
WRITE (1,*) 'EIGEN VECTORS'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , N)
END DO
C------------------- COMPUTATION WITH [T]
CALL POWERMTD ( EPS, MXN, N, T, X, Y, EIGENVEC, EIGEN )
WRITE (1,*)
WRITE (1,*) 'POWER METHOD WITH [T]'
WRITE (1,*) 'MODE EIGENVALUE'
DO I = 1 , N
WRITE (1,*) I , EIGEN(I)
END DO
WRITE (1,*) 'EIGEN VECTORS'
DO I = 1 , N
WRITE (1,*) (EIGENVEC(I,J), J= 1 , N)
END DO
C=======================================================================
CLOSE (1)
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE BSECTION1 ( EPS,MXENGN, NEIGEN, T, F, EIGEN )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION T(MXENGN,MXENGN), F(MXENGN), EIGEN(MXENGN)
C-------- COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST
BMAX= DMAX1(DABS(T(1,1))+DABS(T(1,2)),
* DABS(T(NEIGEN,NEIGEN))+DABS(T(NEIGEN,NEIGEN-1)) )
DO I=2,NEIGEN-1
BMAX = DMAX1(BMAX,DABS(T(I,I))+
* DABS(T(I,I+1))+DABS(T(I,I-1)))
END DO
C=========================================================================
DO MODE = 1, NEIGEN
A1 = 0.D0
A2 = BMAX
C--- BEGIN OF DO WHILE -----
DO WHILE ( DABS(A2-A1)/BMAX .GT. EPS )
A3 = (A1+A2)/2.0D0
C------- FORMING STURM SEQUENCES
C------- P-0, P-1, P-2, ......., P-NEIGEN
F(1) = T(1,1)-A3
F(2) = (T(2,2)-A3)*F(1) - T(1,2)**2
DO J = 3 , NEIGEN
F(J) = (T(J,J)-A3)*F(J-1) - T(J-1,J)**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 EIGENVALUES OF MODE
IF ( NA3 .LT. MODE ) THEN
A1=A3
ELSE
A2=A3
END IF
END DO
C--- END OF DO WHILE -----
EIGEN(MODE)=A3
END DO
RETURN
END
C
C
SUBROUTINE POWERMTD ( EPS, MXN, N, A, X, Y, EIGENVEC, EIGEN )
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXITERA=100 )
DIMENSION A(MXN,MXN), X(MXN),Y(MXN),EIGENVEC(MXN,MXN),EIGEN(MXN)
DO MODE = 1 , N
C
IF ( MODE .GT. 1 ) THEN
DO I = 1 , N
DO J = 1 , N
A(I,J) = A(I,J) -
* EIGEN(MODE-1)*EIGENVEC(I,MODE-1)*EIGENVEC(J,MODE-1)
END DO
END DO
END IF
C
EIGEN(MODE) = 1.D0
DO I = 1 , N
X(I) = 0.D0
END DO
X(1) = 1.D0
DIFF = 1.D0
NITERA = 0
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
DO WHILE ( DIFF .GT. EPS .AND. NITERA .LE. MXITERA )
NITERA = NITERA + 1
C-------- COMPUTATION OF [A]{X}={Y}
DO I = 1 , N
SUM = 0.D0
DO J = 1 , N
SUM = SUM + A(I,J)*X(J)
END DO
Y(I) = SUM
END DO
C-------- COMPUTATION OF EIGENVALUE
SUM = 0.D0
DO I = 1 , N
SUM = SUM + Y(I)*Y(I)
END DO
EIGEN(MODE) = DSQRT(SUM)
C-------- COMPUTATION OF ERROR
DIFF = 0.D0
DO I = 1 , N
XNEW = Y(I) / EIGEN(MODE)
DIFF = DIFF + DABS( X(I) - XNEW )
X(I) = XNEW
END DO
DIFF = DIFF / N
END DO
C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
DO I = 1 , N
EIGENVEC(I,MODE) = X(I)
END DO
C
END DO
RETURN
END
C
C
SUBROUTINE JACOBS ( MXN, EPS, MAXTRTN, NNODE, A, U, FLAMBDA )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION A(MXN,MXN),U(MXN,MXN), FLAMBDA(MXN)
C==================== INITIALIZATION OF EIGEN VECTORS ==================
DO I=1,NNODE
DO J=1,NNODE
U(I,J) = 0.0D0
END DO
U(I,I) = 1.0D0
END DO
C======================INITIALIZATION FOR DO WHILE =====================
ITERATION = 0
AMXOFF = EPS
DMX = EPS
C================== ITERATION PROCEDURE STARTS HERE ====================
DO WHILE ( (AMXOFF/DMX.GT.EPS) .AND. (ITERATION.LT.MAXTRTN) )
DO I = 1 , NNODE
IF ( DABS(A(I,I)) .GT. DMX ) DMX = DABS(A(I,I))
END DO
C--------- FIND MAXIMUM IN OFF DIAGONAL ELEMENTS OF A(I,J) -------------
AMXOFF = 0.D0
DO I = 1 , NNODE-1
DO J = I+1 , NNODE
IF ( DABS(A(I,J)) .GT. AMXOFF ) THEN
AMXOFF = DABS(A(I,J))
K = I
L = J
END IF
END DO
END DO
IF ( AMXOFF .EQ. 0.D0 ) EXIT
C------------------ COMPUTATION OF ROTATIONAL ANGLE --------------------
C------------------ R IMPLIES RADIUS OF MOHR CIRCLE --------------------
C--------- CENTER IMPLIES CENTER COORDINATE OF MOHR CIRCLE -------------
ITERATION = ITERATION + 1
DIFFHALF = (A(K,K) - A(L,L))/2.D0
IF ( DIFFHALF .EQ. 0.D0 ) THEN
R = DABS(A(K,L))
COSINE = 1.D0/DSQRT(2.D0)
SINE = 1.D0/DSQRT(2.D0)
IF ( A(K,L) .LT. 0.D0 ) SINE = -SINE
IF ( R .EQ. 0.D0 ) THEN
COSINE = 1.D0
SINE = 0.D0
END IF
ELSE
R = DSQRT ( DIFFHALF**2 + A(K,L)**2 )
COSINE = DSQRT ( DABS(DIFFHALF)/(2.D0*R) + 0.5D0 )
SINE = A(K,L)/(2.D0*R*COSINE)
IF ( DIFFHALF .LT. 0.D0 ) COSINE = -COSINE
END IF
C------------------------- ORTHOGNALIZATION ----------------------------
IF ( R .GT. 0.D0 ) THEN
DO I = 1 , NNODE
AW = U(I,K)
AZ = U(I,L)
U(I,K) = AW*COSINE + AZ*SINE
U(I,L) = -AW*SINE + AZ*COSINE
IF ( ( I .NE. K ) .AND. ( I .NE. L ) ) THEN
AW = A(I,K)
AZ = A(I,L)
A(I,K) = AW*COSINE + AZ*SINE
A(I,L) = -AW*SINE + AZ*COSINE
A(K,I) = A(I,K)
A(L,I) = A(I,L)
END IF
END DO
A(K,L) = 0.D0
A(L,K) = 0.D0
CENTER = 0.5D0*(A(K,K)+A(L,L))
IF ( DIFFHALF .LT. 0.D0 ) R = -R
A(K,K) = CENTER + R
A(L,L) = CENTER - R
C----------------- LINES BELOW ALSO WORKS FINE.------------------------
C AKK = A(K,K)
C ALL = A(L,L)
C CROSS = 2.D0*A(K,L)*SINE*COSINE
C A(K,K) = AKK*COSINE**2 + ALL*SINE**2 + CROSS
C A(L,L) = AKK*SINE**2 + ALL*COSINE**2 - CROSS
C----------------------------------------------------------------------
END IF
END DO
C======================== EIGENVALUES = A(I,I) ========================
DO I = 1 , NNODE
FLAMBDA(I) = 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