PROGRAM PRINCIPAL_STRESS_SINGLE_TRY
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION TAUJ(3), VECJ(3,3)
DIMENSION TAUP(3), VECP(3,3)
DIMENSION TAUN(3), VECN(3,3)
DIMENSION TAUC(3), VECC(3,3)
DIMENSION TAUB(3), VECB(3,3)
C
EPS = 1.D-14
C
TXX = 1.04179775917145D-10
TYY = 2.97717406283481D-12
TZZ = -4.26645385687152D-11
TXY = -7.22008219763553D-11
TXZ = -0.775248853114106
TYZ = 2.51532128459075D-11
C
C TXX = 3.D0
C TYY = 3.D0
C TZZ = 3.D0
C TXY = 0.D0
C TXZ = 0.D0
C TYZ = 0.D0
C
C----- PRINCIPAL STRESSES TAU'S ----------
CALL JACOBI ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICJ, TAUJ, VECJ )
CALL POWER ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICP, TAUP, VECP )
CALL NEWTON ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICN, TAUN, VECN )
CALL CARDANO ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICC, TAUC, VECC )
CALL BSECTION ( EPS,TXX,TYY,TZZ,TXY,TXZ,TYZ, ICB, TAUB, VECB )
C
OPEN (1,FILE='PRINCIPAL-STRESSES-SINGLE.OUT',STATUS='UNKNOWN')
C----- FILE--------VECX(J,MODE)
WRITE(1,*) 'METHOD T1 T2 T3 A1X A1Y A1Z A2X A2Y A2Z A3X A3Y A3Z',
* ' ITERATION'
WRITE(1,*)'CARDANO ',(TAUC(J),J=1,3),((VECC(J,M),J=1,3),M=1,3),ICC
WRITE(1,*) 'NEWTON ',(TAUN(J),J=1,3),((VECN(J,M),J=1,3),M=1,3),ICN
WRITE(1,*) 'JACOBI ',(TAUJ(J),J=1,3),((VECJ(J,M),J=1,3),M=1,3),ICJ
WRITE(1,*) 'POWER ' ,(TAUP(J),J=1,3),((VECP(J,M),J=1,3),M=1,3),ICP
WRITE(1,*) 'BISECTION ',(TAUB(J),J=1,3),((VECB(J,M),J=1,3),M=1,3),
* ICB
STOP' NORMAL TERMINATION'
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
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
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