PROGRAM JACOB35
C==================================================================
C           COMPUTATION OF EIGENVALUES AND EIGENVECTORS USING
C                       JACOB METHOD
C        EQUATION TO BE SOLVED: [A]{X}=LAMBDA{X}
C            CONDITION: [A]= SYMMETRIC AND POSITIVE DIFINITE
C              JULY 18, 2008 EIJI FUKUMORI
C==================================================================
      IMPLICIT REAL*8(A-H,O-Z)
C
      PARAMETER (MXN=160, EPS = 1.D-15, ITAMAX = 50)
      DIMENSION A(MXN,MXN),U(MXN,MXN), FLAMBDA(MXN)
C------ BELOW IS A SET OF TEST DATA.
C    NORNALLY TAKEN FROM FILE
      NNODE = 3
      DATA A(1,1), A(1,2), A(1,3) / 3.D0, 2.D0, 1.D0 /
      DATA A(2,1), A(2,2), A(2,3) / 2.D0, 2.D0, 1.D0 /
      DATA A(3,1), A(3,2), A(3,3) / 1.D0, 1.D0, 1.D0 /
C------- COMPUTE EIGENVALUES AND VECTORS ---------
      MAXTRTN = NNODE*ITAMAX
      CALL JACOBS ( MXN, EPS, MAXTRTN, NNODE, A, U, FLAMBDA )
C------------ PRINT
      OPEN ( 1, FILE='JACOB.OUT', STATUS='UNKNOWN' )
      DO J = 1 , NNODE
      WRITE (1,*) 'EIGENVALUE (',J,') = ',  FLAMBDA(J)
      WRITE(1,*) 'EIGENVECTOR'
      DO I = 1 , NNODE
      WRITE(1,*)  U(I,J)
      END DO
      WRITE (1,*)
      END DO
      CLOSE (1)
      STOP 'NORMAL TERMINATION'
      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