PROGRAM LANCZOS_PRINCIPLE_BASIC2
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                      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), F(MXN), EIGEN(MXN)
C=======================================================================
      OPEN (1,FILE='LANCZOS-BASIC2.OUT', STATUS='UNKNOWN')
      WRITE (1,*) '** SOLUTION OF LANCZOS METHOD/BISECTION **'
      WRITE (1,*) '******* m = n *******'
      WRITE (1,*) 'EIGENVALUES BY BISECTION2'
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 BSECTION2 ( EPS, MXN, N, ALPHA,BETA, F, EIGEN )
      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 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