PROGRAM LANCZOS_PRINCIPLE_BASIC
C=======================================================================
C      FORMATION OF ORTHONORMAL VECTOR AND TRI-DAIGONAL MATRIX BY
C                           LANCZOS PRINCIPAL
C           EIGENVALUES ARE SOLVED BY THE FOLLOWING METHOD:
C                           BISECTION METHOD
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), T(MXN,MXN)
C------- BISECTION METHOD MEMORY
      DIMENSION F(MXN)
C------- COMMON
      DIMENSION EIGEN(MXN)
C=======================================================================
      OPEN (1,FILE='LANCZOS-BASIC1.OUT', STATUS='UNKNOWN')
      WRITE (1,*) '** FORMATION OF [T] BY LANCZOS METHOD **'
      WRITE (1,*) '******* m = n *******'
      WRITE (1,*) 'EIGENVALUES BY BISECTION1'
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=[U(I,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(I,J)] BY LANCZOS PRINCIPAL'
      DO I = 1 , N
      WRITE (1,*) ( U(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-------- 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
C=======================================================================
C------- BISECTION METHOD
      CALL BSECTION1 ( EPS, MXN, N, T, F, EIGEN )
      WRITE (1,*)
      WRITE (1,*) 'EIGENVALUES BY BISECTION METHOD WITH [T]'
      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 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