PROGRAM LANCZOS_PRINCIPLE_BASIC3
C=======================================================================
C FORMATION OF ORTHONORMAL VECTOR AND TRI-DAIGONAL MATRIX BY
C LANCZOS PRINCIPAL
C EIGENVALUES ARE SOLVED BY BISECTION METHOD
C NOTE: BISECTION METHOD EMPLOYS ALPHAS AND BETAS
C NO [T] MATRIX ASSEMBLED
C***********************************************************************
C NEWTON-RAPHSON METHOD IS UTILIZED IN BISECTION METHOD TO SPEED UP
C EIGENVALUE COMPUTATION
C***********************************************************************
C 2011/JAN/25 EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( MXN=100, EPS=1.D-13 )
C------- LANCZOS MEMORY
DIMENSION A(MXN,MXN), U(MXN,MXN), ALPHA(MXN), BETA(MXN),
* AU1(MXN), R(MXN), F(MXN), EIGEN(MXN), P(MXN)
C=======================================================================
OPEN (1,FILE='LANCZOS-BASIC3.OUT', STATUS='UNKNOWN')
WRITE (1,*) '** SOLUTION OF LANCZOS METHOD/BISECTION/NEWTON **'
WRITE (1,*) '******* m = n *******'
WRITE (1,*) 'EIGENVALUES BY BISECTION3'
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=======================================================================
C COMPUTATION OF IGENVALUES BY BISECTION METHOD METHOD
C=======================================================================
C------- BISECTION METHOD
CALL BSECTION3 (EPS,MXN,N, EIGEN,P,ALPHA,BETA )
WRITE (1,*)
WRITE (1,*) 'BISECTION METHOD WITH ALPHA AND BETA VALUES'
WRITE (1,*) 'MODE EIGENVALUES'
DO I = 1 , N
WRITE (1,*) I, EIGEN(I)
END DO
C=======================================================================
CLOSE (1)
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE BSECTION3 (EPS,MXIGEN,NLANCZOS,EIGEN,P,ALPHA,BETA )
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ALPHA(MXIGEN), BETA(MXIGEN), P(MXIGEN), EIGEN(MXIGEN)
LOGICAL DCSN1, DCSN2, DCSN3
C------------ COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST ------------
WRITE (1,*)
WRITE (1,*) 'ERROR JUDGMENT VALUE(EPS)=', EPS
WRITE (1,*)
WRITE (1,*) 'MODE METHOD EIGENVALUE'
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) )
ENTRYMX = 0.D0
DO I = 1 , NLANCZOS
ENTRYMX = DMAX1 ( ENTRYMX, DABS(ALPHA(I)), DABS(BETA(I)) )
END DO
ENTRYMX = ENTRYMX
DO I = 1 , NLANCZOS
ALPHA(I) = ALPHA(I)/ENTRYMX
BETA (I) = BETA (I)/ENTRYMX
END DO
C-------------------- COMPUTATION OF EIGENVALUES -----------------------
DO MODE = 1, NLANCZOS
BOTTOLMT =-BANDMAX/ENTRYMX
UPPERLMT = BANDMAX/ENTRYMX
C------------ INITIAL SETTING FOR CONVERGENCE PARAMETERS ---------------
N1 = 0
N2 = NLANCZOS
BASE = UPPERLMT
DCSN1 = .FALSE.
DCSN2 = .TRUE.
DCSN3 = .FALSE.
C------- BEGIN OF DO WHILE -----
DO WHILE ( DCSN2 )
C------------ PREDICTION OF NEW LAMBDA BY BISECTION METHOD -------------
FLAMBDA = 0.5D0*( BOTTOLMT + UPPERLMT )
C---------- PREDICTION OF NEW LAMBDA BY NEWTON-RAPHSON METHOD ----------
IF ( DCSN1 .AND. DCSN3 ) THEN
FLAMBDA = FLNEWTON
WRITE (1,*) MODE, "NEWTON-PAPHSON ",FLAMBDA*ENTRYMX
ELSE
WRITE (1,*) MODE, "BISECTION ",FLAMBDA*ENTRYMX
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
DENOMI = DABS ( DMAX1(BOTTOLMT,UPPERLMT) )
IF ( DENOMI .EQ. 0.D0 ) DENOMI = BANDMAX/ENTRYMX
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)/DENOMI .GT. EPS
DCSN3 = .FALSE.
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/DENOMI) .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*ENTRYMX
IF ( DCSN3 ) EIGEN(MODE) = FLNEWTON*ENTRYMX
C----------------- END OF DO MODE --------------------
END DO
DO I = 1 , NLANCZOS
ALPHA(I) = ALPHA(I)*ENTRYMX
BETA (I) = BETA (I)*ENTRYMX
END DO
RETURN
END