PROGRAM PRINCIPAL_STRESS
C=======================================================================
C   COMPUTATION OF PRINCIPAL STESSES USING various METHODS
C                           EIJI FUKUMORI
C                         30 DECEMBER 2012
C=======================================================================
      IMPLICIT REAL*8 ( A-H , O-Z )
      INCLUDE 'PARAM.DAT'
CCCCC      PARAMETER ( MXE=40000,MXN=53000,MXB=5000, ND=8 )
C=======================================================================
C                                ARRAYS
      CHARACTER INPFILE*14, STRESSF*12
      DIMENSION GLBTAU(7,MXN)
C=======================================================================
      EPS = 1.D-14
C=======================================================================
      INPFILE = 'NONE'
      IF ( ND .EQ.  8 ) THEN
      INPFILE = 'STATIC3D08.DAT'
      STRESSF = 'STRESS08.OUT'
      END IF
      IF ( ND .EQ. 27 ) THEN
      INPFILE = 'STATIC3D27.DAT'
      STRESSF = 'STRESS27.OUT'
      END IF
      IF ( INPFILE .EQ. 'NONE' ) STOP 'NO INPUT FILE AT MAIN'
C=======================================================================
      WRITE(*,*) 'PRINCIPAL STESSES COMPUTATION PROGRAM'
      WRITE (*,*)' READING IN DATA FILES'
      CALL INPUT ( INPFILE,STRESSF,MXN,NNODE,GLBTAU )
C=======================================================================
      WRITE (*,*) 'FILEMK '
      CALL FILEMK ( EPS,ND,MXN, NNODE, GLBTAU )
      STOP' NORMAL TERMINATION'
      END
C
C
      SUBROUTINE FILEMK ( EPS,ND,MXN,NNODE,GLBTAU )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION  GLBTAU(7,MXN)
      DIMENSION TAUJ(3), VECJ(3,3), TAUP(3), VECP(3,3)
      DIMENSION TAUN(3), VECN(3,3), TAUC(3), VECC(3,3)
      DIMENSION TAUB(3), VECB(3,3)
      CHARACTER OUTFILE*15
      LOGICAL YES
C========> GLBTAU(1,I)=TAUXX = MU*(2DU/DX)+LAMBDA*DIVV
C========> GLBTAU(2,I)=TAUYY = MU*(2DV/DY)+LAMBDA*DIVV
C========> GLBTAU(3,I)=TAUZZ = MU*(2DW/DZ)+LAMBDA*DIVV
C========> GLBTAU(4,I)=TAUXY = MU*(DU/DY+DU/DX)
C========> GLBTAU(5,I)=TAUXZ = MU*(DU/DZ+DW/DX)
C========> GLBTAU(6,I)=TAUYZ = MU*(DV/DZ+DW/DY)
C========> GLBTAU(7,I)= DIVV
C--------0---------0---------0---------0---------0---------0---------0--
C
      OUTFILE =  'NONE'
      IF ( ND .EQ.  8 ) THEN
      OUTFILE =  'PRINCIPAL08.OUT'
      END IF
      IF ( ND .EQ. 27 ) THEN
      OUTFILE =  'PRINCIPAL27.OUT'
      END IF
      IF ( OUTFILE .EQ. 'NONE' ) STOP 'NO FILE ASIGNED'
C
      OPEN ( 1, FILE=OUTFILE, STATUS='UNKNOWN' )
      WRITE (1,*) '# TXX TYY TZZ TXY TXZ TYZ ',
     *            ' ITJ T1J T2J T3J ',
     * ' A1XJ A1YJ A1ZJ A2XJ A2YJ A2ZJ A3XJ A3YJ A3ZJ ',
     *            ' ITP T1P T2P T3P ',
     * ' A1XP A1YP A1ZP A2XP A2YP A2ZP A3XP A3YP A3ZP ',
     *            ' ITN T1N T2N T3N ',
     * ' A1XN A1YN A1ZN A2XN A2YN A2ZN A3XN A3YN A3ZN ',
     *            ' ITC T1C T2C T3C ',
     * ' A1XC A1YC A1ZC A2XC A2YC A2ZC A3XC A3YC A3ZC ',
     *            ' ITB T1B T2B T3B ',
     * ' A1XB A1YB A1ZB A2XB A2YB A2ZB A3XB A3YB A3ZB '
C
      DO I = 1 , NNODE
C----- PRINCIPAL STRESSES TAU'S ----------
      CALL JACOBI (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
     *  GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICJ, TAUJ, VECJ )
C
      CALL POWER (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
     *  GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICP, TAUP, VECP )
C
      CALL NEWTON (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
     *  GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICN, TAUN, VECN )
C
      CALL CARDANO (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
     *  GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICC, TAUC, VECC )
C
      CALL BSECTION (EPS,GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
     *  GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), ICB, TAUB, VECB )
C
C----- FILE--------VECX(J,MODE)
      WRITE (1,*)
     * I, (GLBTAU(J,I),J=1,6),
     * ICJ, (TAUJ(J),J=1,3),((VECJ(J,MODE),J=1,3),MODE=1,3),
     * ICP, (TAUP(J),J=1,3),((VECP(J,MODE),J=1,3),MODE=1,3),
     * ICN, (TAUN(J),J=1,3),((VECN(J,MODE),J=1,3),MODE=1,3),
     * ICC, (TAUC(J),J=1,3),((VECC(J,MODE),J=1,3),MODE=1,3),
     * ICB, (TAUB(J),J=1,3),((VECB(J,MODE),J=1,3),MODE=1,3)
      END DO
      CLOSE (1)
      RETURN
      END

C
C
      SUBROUTINE INPUT ( INPFILE,STRESSF,MXN,NNODE,GLBTAU )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION GLBTAU(7,MXN)
      CHARACTER*14 INPFILE*14, STRESSF*12
      LOGICAL YES
C========> FILENAME INPFILE
      OPEN ( 1, FILE=INPFILE, STATUS='UNKNOWN' )
C========> PARAMETERS
      READ (1,*) YOUNG,  POISSON
C========> ELEMENTS
      READ (1,*) NE
      DO I = 1 , NE
      READ (1,*)
      END DO
CC========>
      READ (1,*) NNODE
      CLOSE (1)
C========> FILENAME STRESSF
      OPEN ( 1, FILE=STRESSF, STATUS='UNKNOWN' )
      READ (1,*)
      DO I = 1 , NNODE
      READ (1,*) NODE,(GLBTAU(J,I),J=1,7)
      END DO
      CLOSE (1)
      RETURN
      END
C
C
      SUBROUTINE JACOBI (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C=======================================================================
C                 PRINCIPAL STRESSES USING JACOB METHOD
C=======================================================================
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER ( MAXTRTN=1000 )
      DIMENSION A(3,3),U(3,3), TAU(3), UTEMP(3)
C---------------------- RESETTING EIGEN VECTORS ------------------------
      IC = 0
      NNODE = 3
      DO I=1,NNODE
      DO J=1,NNODE
        U(I,J) = 0.0D0
      END DO 
        U(I,I) = 1.0D0
      END DO
C------------ SUBSTITUTION OF STRESSES INTO [A] MATRIX -----------------
      A(1,1) = TXX
      A(2,2) = TYY
      A(3,3) = TZZ
      A(1,2) = TXY
      A(1,3) = TXZ
      A(2,3) = TYZ
      A(2,1) = A(1,2)
      A(3,1) = A(1,3)
      A(3,2) = A(2,3)
C
      ITERATION = 0
      AMXOFF = EPS
      DMX = EPS
C=======================================================================
      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)-------------------------
      IC = ITERATION
      TAU(1) = A(1,1)
      TAU(2) = A(2,2)
      TAU(3) = A(3,3)
C-------- MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER -----------
      IF ( TAU(2) .GT. TAU(1) ) THEN
          TEMP =   TAU(1)
          TAU(1) = TAU(2)
          TAU(2) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,1)
          U  (J,1) = U  (J,2)
          U  (J,2) = UTEMP(J)
          END DO
      END IF
      IF ( TAU(3) .GT. TAU(2) ) THEN
          TEMP =   TAU(2)
          TAU(2) = TAU(3)
          TAU(3) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,2)
          U  (J,2) = U  (J,3)
          U  (J,3) = UTEMP(J)
          END DO
      END IF
      IF ( TAU(2) .GT. TAU(1) ) THEN
          TEMP =   TAU(1)
          TAU(1) = TAU(2)
          TAU(2) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,1)
          U  (J,1) = U  (J,2)
          U  (J,2) = UTEMP(J)
          END DO
      END IF
      RETURN
      END
C
C
      SUBROUTINE NEWTON (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU,U )
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION TAU(3),DRV(3,3),DI(3,3),F(3),H(3,3),G(3,3),U(3,3)
      DIMENSION TEMPU(3)
      DATA MAXITA / 100 /
C-------------------------- STRESS INVARIANTS --------------------------
      A = TXX+TYY+TZZ
      B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
      C = TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
     * TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY
C---------------- ENERGY BY COMPUTAED GIVEN STRESSES -------------------
      ENERGY0 = (TXX-TYY)**2+(TYY-TZZ)**2+(TZZ-TXX)**2+
     * 6.D0*(TXY**2+TYZ**2+TXZ**2)
C------------------------ RESET EIGEN VECTORS --------------------------
      DO I = 1 , 3
      DO J = 1 , 3
      U(I,J) = 0.D0
      END DO
      U(I,I) = 1.D0
      END DO
C-------------------------- CASE OF A=B=C=0 ----------------------------
      IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) THEN
      IC = 0
      TAU(1) = 0.D0
      TAU(2) = 0.D0
      TAU(3) = 0.D0
      RETURN
      END IF
C--------------------------- CASE OF B=C=0 -----------------------------
      IC = -1
      IF ( B .EQ. 0.D0 .AND. C .EQ. 0.D0 ) THEN
      TAU(2) = 0.D0
      IF ( A .GT. 0.D0 ) THEN
      TAU(1) = A
      TAU(3) = 0.D0
      ELSE
      TAU(1) = 0.D0
      TAU(3) = A
      END IF
      RETURN
      END IF
C---------------------------- CASE OF C=0 ------------------------------
      IC = -2
      IF ( C .EQ. 0.D0 ) THEN
      TAU(1) = (A+DSQRT(A*A-4.D0*B))/2.D0
      TAU(2) = 0.D0
      TAU(3) = (A-DSQRT(A*A-4.D0*B))/2.D0
      IF ( TAU(3) .GT. 0.D0 ) THEN
      TAU(2) = TAU(3)
      TAU(3) = 0.D0
      END IF
      END IF
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
      IC = -3
      P = -(-A)**2/3.D0 + B
      Q = 2.D0*(-A)**3/27.D0 -(-A)*B/3.D0 + (-C)
      IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
      TAU(1) = TXX
      TAU(2) = TYY
      TAU(3) = TZZ
      RETURN
      END IF
C=================== NEWTON-RAPHSON METHOD ITERATION ===================
C--------------------- IC = NUMBER OF ITERATIONS -----------------------
C---------------------- INITIAL GUESSES FOR TAU'S ----------------------
      TAU(1) = DSQRT(2.D0) + A/3.D0
      TAU(2) = DSQRT(3.D0) + B/TAU(1)
      TAU(3) = DATAN(1.D0) + C/(TAU(1)*TAU(2))
C
      DO I = 1 , MAXITA
      IC = I
      F(1) = TAU(1)+TAU(2)+TAU(3)-A
      F(2) = TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1)-B
      F(3) = TAU(1)*TAU(2)*TAU(3)-C
C--------- DERIVATIVE MATRIX [DRV] ------
      DRV(1,1) = 1.D0
      DRV(1,2) = 1.D0
      DRV(1,3) = 1.D0
      DRV(2,1) = TAU(2)+TAU(3)
      DRV(2,2) = TAU(1)+TAU(3)
      DRV(2,3) = TAU(1)+TAU(2)
      DRV(3,1) = TAU(2)*TAU(3)
      DRV(3,2) = TAU(1)*TAU(3)
      DRV(3,3) = TAU(1)*TAU(2)
C------- DETERMINANT OF [DRV]
      DETJ = DRV(1,1) * ( DRV(2,2)*DRV(3,3)-DRV(2,3)*DRV(3,2) )
     *     - DRV(1,2) * ( DRV(2,1)*DRV(3,3)-DRV(2,3)*DRV(3,1) )
     *     + DRV(1,3) * ( DRV(2,1)*DRV(3,2)-DRV(2,2)*DRV(3,1) )
C------- [DI] = INVERSE OF THE MATRIX [DRV]
      DI(1,1) = (DRV(2,2)*DRV(3,3) - DRV(2,3)*DRV(3,2))/DETJ
      DI(2,1) = (DRV(2,3)*DRV(3,1) - DRV(2,1)*DRV(3,3))/DETJ
      DI(3,1) = (DRV(2,1)*DRV(3,2) - DRV(2,2)*DRV(3,1))/DETJ
      DI(1,2) = (DRV(1,3)*DRV(3,2) - DRV(1,2)*DRV(3,3))/DETJ
      DI(2,2) = (DRV(1,1)*DRV(3,3) - DRV(1,3)*DRV(3,1))/DETJ
      DI(3,2) = (DRV(1,2)*DRV(3,1) - DRV(1,1)*DRV(3,2))/DETJ
      DI(1,3) = (DRV(1,2)*DRV(2,3) - DRV(1,3)*DRV(2,2))/DETJ
      DI(2,3) = (DRV(1,3)*DRV(2,1) - DRV(1,1)*DRV(2,3))/DETJ
      DI(3,3) = (DRV(1,1)*DRV(2,2) - DRV(1,2)*DRV(2,1))/DETJ
C-------- NEW VALUES FOR TAU'S
      TAU(1) = TAU(1) - ( DI(1,1)*F(1)+DI(1,2)*F(2)+DI(1,3)*F(3) )
      TAU(2) = TAU(2) - ( DI(2,1)*F(1)+DI(2,2)*F(2)+DI(2,3)*F(3) )
      TAU(3) = TAU(3) - ( DI(3,1)*F(1)+DI(3,2)*F(2)+DI(3,3)*F(3) )
C--------- ENERGY BY COMPUTAED PRINCIPAL STRESSES ------
      ENERGY1 = 2.D0*(TAU(1)+TAU(2)+TAU(3))**2-
     *      6.D0*(TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1))
C--------- DECISION --------
      IF ( DABS(ENERGY1-ENERGY0)/ENERGY0 .LE. EPS ) EXIT
      END DO
C--------------------------- NO CONVERGENCE ----------------------------
      IF ( I .GE. MAXITA ) THEN
      IC = -4
      END IF
C===================== EIGEN VECTOR COMPUTATION ========================
      CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICV, TAU, U )
      IF ( IC .GT. 0 ) IC = IC + ICV
C============ MAKING TAU(i) AND U(i,j) DESCENDING ORDER ================
      CALL ORDERING ( TAU, U )
      RETURN
      END
C
C
      SUBROUTINE CARDANO (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C=========== SOLUTION OF CUBIC EQUATION BY CARDANO METHOD ==============
C-------------------- X**3 + Ax**2 + Bx +C = 0 -------------------------
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION TAU(3), U(3,3)
      COMPLEX*16 T1, T2, ALPHA, X1, X2, X3, BETA
C------------------------ STRESS INVARIANTS ----------------------------
      A = -(TXX+TYY+TZZ)
      B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
      C = -(TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
     * TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY)
C------------------------- INITIAL SETTING -----------------------------
      DO I = 1 , 3
      DO J = 1 , 3
      U(I,J) = 0.D0
      END DO
      U(I,I) = 1.D0
      END DO
C------------------------- CASE OF A=B=C=0 -----------------------------
      IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) THEN
      IC = 0
      TAU(1) = 0.D0
      TAU(2) = 0.D0
      TAU(3) = 0.D0
      RETURN
      END IF
C------------------------ VARIABLES P AND Q ----------------------------
      P = -A**2/3.D0 + B
      Q = 2.D0*A**3/27.D0 -A*B/3.D0 + C
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
      IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
      IC = -3
      TAU(1) = TXX
      TAU(2) = TYY
      TAU(3) = TZZ
      RETURN
      END IF
C------------------- COMPUTATION OF QUBIC ROOTS ------------------------
      DISCR = Q**2 + (4.D0/27.D0)*P**3
      BETA = (-DCMPLX(Q) - CDSQRT(DCMPLX(DISCR)))/2.0D0
      ALPHA = BETA**(1.D0/3.D0)
      T1 = DCMPLX(-0.5D0, +0.5D0*DSQRT(3.D0))
      T2 = DCMPLX(-0.5D0, -0.5D0*DSQRT(3.D0))
      X1 = ALPHA   -DCMPLX(P)/(3.D0*ALPHA)   -DCMPLX(A)/3.D0
      X2 = ALPHA*T1-DCMPLX(P)/(3.D0*ALPHA*T1)-DCMPLX(A)/3.D0
      X3 = ALPHA*T2-DCMPLX(P)/(3.D0*ALPHA*T2)-DCMPLX(A)/3.D0
C----------------- ALL PRINCIPAL STRESSES ARE REAL.---------------------
      TAU(1) = DREAL(X1)
      TAU(2) = DREAL(X2)
      TAU(3) = DREAL(X3)
C====================== EIGEN VECTOR COMPUTATION =======================
      CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C========== MAKING TAU(i) ALONG WITH U(i,j) DESCENDING ORDER ===========
      CALL ORDERING ( TAU, U )
      RETURN
      END
C
C
      SUBROUTINE POWER (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ,NITERA,TAU,U)
C---------- COMPUTATION OF PRINCIPAL STRESSES BY POWER METHOD ----------
      IMPLICIT REAL*8 ( A-H , O-Z )
      PARAMETER ( MXITERA=99999 )
      DIMENSION A(3,3), X(3),Y(3),U(3,3),TAU(3), UTEMP(3)
      N = 3
C--------------- CASE OF P=Q=0 MEANING TRIPLE ROOTS---------------------
      AA = -(TXX+TYY+TZZ)
      B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
      C = -(TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
     * TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY)
      P = -AA**2/3.D0 + B
      Q = 2.D0*AA**3/27.D0 -AA*B/3.D0 + C
      IF ( P.EQ.0.D0 .AND. Q.EQ.0.D0 ) THEN
      TAU(1) = TXX
      TAU(2) = TYY
      TAU(3) = TZZ
      DO I=1,3
      DO J=1,3
        U(I,J) = 0.0D0
      END DO 
        U(I,I) = 1.0D0
      END DO
      NITERA = -3
      RETURN
      END IF
C----------------- SUBSTITUTION OF TAU-IJ INTO A(I,J) ------------------
      A(1,1) = TXX
      A(2,2) = TYY
      A(3,3) = TZZ
      A(1,2) = TXY
      A(1,3) = TXZ
      A(2,3) = TYZ
      A(2,1) = A(1,2)
      A(3,1) = A(1,3)
      A(3,2) = A(2,3)
C----------------- CASE OF ALL STRESSES ARE ZERO.  IC=0 ----------------
      SUM = 0.D0
      DO I = 1 , N
      DO J = 1 , N
      SUM = SUM + DABS(A(I,J))
      END DO
      END DO
      IF ( SUM .EQ. 0.D0 ) THEN
         TAU(1) = 0.D0
         TAU(2) = 0.D0
         TAU(3) = 0.D0
         IC = 0
         RETURN
      END IF
C================ ADDING BIAS STRESS INTO NORMAL STRESSES ==============
      TAUMIN =               DABS(A(1,1))+DABS(A(1,2))+DABS(A(1,3))
      TAUMIN = DMAX1 (TAUMIN,DABS(A(2,1))+DABS(A(2,2))+DABS(A(2,3)))
      TAUMIN = DMAX1 (TAUMIN,DABS(A(3,1))+DABS(A(3,2))+DABS(A(3,3)))
      A(1,1) = A(1,1) + 2.D0*TAUMIN
      A(2,2) = A(2,2) + 2.D0*TAUMIN
      A(3,3) = A(3,3) + 2.D0*TAUMIN
C============== PRINCIPAL STRESS COMPUTATION STARTS HERE ===============
      NITERA = 0
      DO MODE = 1 , 3
       IF ( MODE .GT. 1 ) THEN
C-------------------- FOR SECOND AND THIRD MODES -----------------------
          A(1,1) = A(1,1) - TAU(MODE-1)*U(1,MODE-1)*U(1,MODE-1)
          A(1,2) = A(1,2) - TAU(MODE-1)*U(1,MODE-1)*U(2,MODE-1)
          A(1,3) = A(1,3) - TAU(MODE-1)*U(1,MODE-1)*U(3,MODE-1)
          A(2,1) = A(2,1) - TAU(MODE-1)*U(2,MODE-1)*U(1,MODE-1)
          A(2,2) = A(2,2) - TAU(MODE-1)*U(2,MODE-1)*U(2,MODE-1)
          A(2,3) = A(2,3) - TAU(MODE-1)*U(2,MODE-1)*U(3,MODE-1)
          A(3,1) = A(3,1) - TAU(MODE-1)*U(3,MODE-1)*U(1,MODE-1)
          A(3,2) = A(3,2) - TAU(MODE-1)*U(3,MODE-1)*U(2,MODE-1)
          A(3,3) = A(3,3) - TAU(MODE-1)*U(3,MODE-1)*U(3,MODE-1)
       END IF
C----------- INITIAL SETTING FOR ITERATION PROCEDURE -------------------
      U(1,MODE) = 0.D0
      U(2,MODE) = 0.D0
      U(3,MODE) = 0.D0
      U(MODE,MODE) = 1.D0
      DIFF = 1.D0
C--------------- ITERATION FOR MODE=I STARTS HERE ----------------------
      DO WHILE ( DIFF .GT. EPS .AND. NITERA .LE. MXITERA )
      NITERA = NITERA + 1
C-------- COMPUTATION OF [A]{X}={Y}
      Y(1) = A(1,1)*U(1,MODE) + A(1,2)*U(2,MODE) + A(1,3)*U(3,MODE)
      Y(2) = A(2,1)*U(1,MODE) + A(2,2)*U(2,MODE) + A(2,3)*U(3,MODE)
      Y(3) = A(3,1)*U(1,MODE) + A(3,2)*U(2,MODE) + A(3,3)*U(3,MODE)
C-------- COMPUTATION OF PRINCIPAL STRESSES
      TAU(MODE) = DSQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3))
C-------- COMPUTATION OF ERROR
      IF ( TAU(MODE) .NE. 0.D0 ) THEN
         X(1) = Y(1)/TAU(MODE)
         X(2) = Y(2)/TAU(MODE)
         X(3) = Y(3)/TAU(MODE)
         DIFF=(DABS(U(1,MODE)-X(1))+DABS(U(2,MODE)-X(2))+
     *    DABS(U(3,MODE)-X(3)))/3
         U(1,MODE) = X(1)
         U(2,MODE) = X(2)
         U(3,MODE) = X(3)
      ELSE
         DIFF = 0.D0
         U(1,MODE) = 0.D0
         U(2,MODE) = 0.D0
         U(3,MODE) = 0.D0
         U(MODE,MODE) = 1.D0
      END IF
C----------------- END OF ITERATION FOR MODE=I -------------------------
      END DO
C------------------------ END OF DO MODE -------------------------------
      END DO
C============== SUBTRACT BIAS FROM NORMAL STRESSES =====================
      TAU(1) = TAU(1) - 2.D0*TAUMIN
      TAU(2) = TAU(2) - 2.D0*TAUMIN
      TAU(3) = TAU(3) - 2.D0*TAUMIN
C======== MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER ===========
      CALL ORDERING ( TAU, U )
      RETURN
      END
C
C
      SUBROUTINE BSECTION ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC,TAU,U )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION P(3), TAU(3), U(3,3)
      LOGICAL DCSN2, DCSN5
      MAXITA = 1000
      WRITE (2,*) '@BSECTION2'
      WRITE (2,*) 'MODE BISECTION-ITERATION LAMBDA ERROR'
C-------- 
      NEIGEN = 3
C-------- COMPUTATION OF LIMITS WHERE EIGENVALUES EXIST
      TAUMAX =               DABS(TXX)+DABS(TXY)+DABS(TXZ)
      TAUMAX = DMAX1 (TAUMAX,DABS(TXY)+DABS(TYY)+DABS(TYZ))
      TAUMAX = DMAX1 (TAUMAX,DABS(TXZ)+DABS(TYZ)+DABS(TZZ))
C-------- COMPUTATION OF IGENVALUES ----------
      DO MODE = 1, NEIGEN
       ITA = 0
       BOTTOLMT =-TAUMAX*1.1D0
       UPPERLMT = TAUMAX*1.2D0
       FMAG = 1.D0
C------------ INITIAL SETTING FOR CONVERGENCE PARAMETERS ---------------
       DCSN2 = .TRUE.
       DCSN5 = .TRUE.
C*** BEGIN  OF DO WHILE ***
      DO WHILE ( DCSN2 .AND. DCSN5 )
       DCSN2 = DABS(UPPERLMT-BOTTOLMT)/FMAG .GT. EPS 
       ITA = ITA + 1
       DCSN5 = ITA .LT. MAXITA
       FLMD = 0.5D0*(BOTTOLMT+UPPERLMT)
C------- FORMING  STURM SEQUENCES
C------- P(0), P(1), P(2), P(3).  NOTE: P(0)=1.
       P(1) =  TXX-FLMD
       P(2) = (TXX-FLMD)*(TYY-FLMD)-TXY**2
       POSI = (TXX-FLMD)*(TYY-FLMD)*(TZZ-FLMD)+2.D0*TXY*TYZ*TXZ
       ANEG = (TYY-FLMD)*TXZ**2+(TXX-FLMD)*TYZ**2+(TZZ-FLMD)*TXY**2
       P(3) = POSI - ANEG
C------- COUNT NUMBER OF SIGN-CHANGES IN P(I)
       KOUNT = 0
       PRODUCT = 1.D0
       DO J = 1 , NEIGEN
       IF ( PRODUCT*P(J) .LT. 0.D0 ) THEN
       KOUNT = KOUNT + 1
       PRODUCT = - PRODUCT
       END IF
       END DO
C------- NARROW DOWN LIMTS WHERE EIGENVALUE OF MODE EXISTS
       IF ( KOUNT .LT. MODE ) THEN
           BOTTOLMT = FLMD
       ELSE
           UPPERLMT = FLMD
       END IF
       FMAG = DMAX1 (DABS(UPPERLMT), DABS(BOTTOLMT))
       IF ( FMAG .EQ. 0.D0 ) FMAG = 1.D0
      END DO
C*** END  OF DO WHILE ***
        TAU(MODE) = FLMD
        WRITE (2,*) MODE,ITA,FLMD,DABS(UPPERLMT-BOTTOLMT)
      END DO
C====================== EIGEN VECTOR COMPUTATION =======================
      CALL EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
      IC = IC + ITA
C========== MAKING TAU(i) ALONG WITH U(i,j) DESCENDING ORDER ===========
      CALL ORDERING ( TAU, U )
      RETURN
      END
C
C
      SUBROUTINE EVECTOR (EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU, U )
C========= SOLUTION OF EIGENVECTORS BY INVERSE POWER METHOD ============
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION TAU(3), F(3,3), G(3,3), U(3,3), TEMPU(3)
      DATA MAXITA / 100 /
C
      DO MODE = 1 , 3
      F(1,1) = TXX-TAU(MODE)*(1.D0-EPS)
      F(2,2) = TYY-TAU(MODE)*(1.D0-EPS)
      F(3,3) = TZZ-TAU(MODE)*(1.D0-EPS)
      F(1,2) = TXY
      F(1,3) = TXZ
      F(2,3) = TYZ
      F(2,1) = F(1,2)
      F(3,1) = F(1,3)
      F(3,2) = F(2,3)
C----------- DETERMINANT OF [F]=[TAU(I,J)]-TAUI(1-EPS)[I] --------------
      POSI = F(1,1)*F(2,2)*F(3,3)+F(2,1)*F(3,2)*F(1,3)+
     *       F(3,1)*F(1,2)*F(2,3)
      ANEG = F(1,3)*F(2,2)*F(3,1)+F(1,1)*F(2,3)*F(3,2)+
     *       F(2,1)*F(1,2)*F(3,3)
      DET = POSI - ANEG
C------- 1ST ROW OF INVERSE OF [F]
      G(1,1) =  (F(2,2)*F(3,3)-F(2,3)*F(3,2))/DET
      G(1,2) = -(F(2,1)*F(3,3)-F(2,3)*F(3,1))/DET
      G(1,3) =  (F(2,1)*F(3,2)-F(2,2)*F(3,1))/DET
C------- 2ND ROW OF INVERSE OF [F]
      G(2,1) = -(F(1,2)*F(3,3)-F(1,3)*F(3,2))/DET
      G(2,2) =  (F(1,1)*F(3,3)-F(1,3)*F(3,1))/DET
      G(2,3) = -(F(1,1)*F(3,2)-F(1,2)*F(3,1))/DET
C------- 3RD ROW OF INVERSE OF [F]
      G(3,1) =  (F(1,2)*F(2,3)-F(1,3)*F(2,2))/DET
      G(3,2) = -(F(1,1)*F(2,3)-F(1,3)*F(2,1))/DET
      G(3,3) =  (F(1,1)*F(2,2)-F(1,2)*F(2,1))/DET
C------- INVERSE OF [F] TIMES TAU*EPS
      DO I = 1 , 3
      DO J = 1 , 3
      G(I,J) = TAU(MODE)*EPS*G(I,J)
      END DO
      END DO
C------------- INITIAL SETTING OF EIGEN VECTOR U(J,MODE) ---------------
      U(1,MODE) = 0.D0
      U(2,MODE) = 0.D0
      U(3,MODE) = 0.D0
      U(MODE,MODE) = 1.D0
C============================= ITERATION ===============================
      VECL = 9999.D0
      DO ITERATIN = 1 , MAXITA
       TEMPU(1) = G(1,1)*U(1,MODE)+G(1,2)*U(2,MODE)+G(1,3)*U(3,MODE)
       TEMPU(2) = G(2,1)*U(1,MODE)+G(2,2)*U(2,MODE)+G(2,3)*U(3,MODE)
       TEMPU(3) = G(3,1)*U(1,MODE)+G(3,2)*U(2,MODE)+G(3,3)*U(3,MODE)
C--------------------------- VECTOR LENGTH -----------------------------
       VECLENGH = DSQRT(TEMPU(1)**2+TEMPU(2)**2+TEMPU(3)**2)
C--------------------------- SUBSTITUTION ------------------------------
       U(1,MODE) = TEMPU(1)/VECLENGH
       U(2,MODE) = TEMPU(2)/VECLENGH
       U(3,MODE) = TEMPU(3)/VECLENGH
C------------------------- CHECK CONVERGENCE ---------------------------
       DELTA = DABS(VECL-VECLENGH)/DMAX1(VECL,VECLENGH)
       IF ( DELTA .LE. EPS ) EXIT
       VECL = VECLENGH
      END DO
C------------------------ BELOW END OF DO MODE -------------------------
      END DO
      IC = ITERATIN
      RETURN
      END
C
C
      SUBROUTINE ORDERING ( TAU, U )
C======== MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER ===========
      IMPLICIT REAL*8 ( A-H , O-Z )
      DIMENSION U(3,3),TAU(3), UTEMP(3)
      IF ( TAU(2) .GT. TAU(1) ) THEN
          TEMP =   TAU(1)
          TAU(1) = TAU(2)
          TAU(2) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,1)
          U  (J,1) = U  (J,2)
          U  (J,2) = UTEMP(J)
          END DO
      END IF
      IF ( TAU(3) .GT. TAU(2) ) THEN
          TEMP =   TAU(2)
          TAU(2) = TAU(3)
          TAU(3) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,2)
          U  (J,2) = U  (J,3)
          U  (J,3) = UTEMP(J)
          END DO
      END IF
      IF ( TAU(2) .GT. TAU(1) ) THEN
          TEMP =   TAU(1)
          TAU(1) = TAU(2)
          TAU(2) = TEMP
          DO J = 1 , 3
          UTEMP(J) = U  (J,1)
          U  (J,1) = U  (J,2)
          U  (J,2) = UTEMP(J)
          END DO
      END IF
      RETURN
      END