PROGRAM LANCZOS_PRINCIPLE1
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-PRINCIPLE1.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 = 3
C--------- GIVEN [A]
      A(1,1) = 3.D0
      A(1,2) = 2.D0
      A(1,3) = 1.D0
      A(2,1) = 2.D0
      A(2,2) = 2.D0
      A(2,3) = 1.D0
      A(3,1) = 1.D0
      A(3,2) = 1.D0
      A(3,3) = 1.D0
      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
      U(1,1) = 1.D0/DSQRT(2.D0)
      U(2,1) = 0.D0
      U(3,1) = 1.D0/DSQRT(2.D0)
      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