+++ /dev/null
- SUBROUTINE DJACOB(N,NMAX,MAXJAC,E,A,C,AII)
- IMPLICIT REAL*8 (A-H,O-Z)
-C THE JACOBI DIAGONALIZATION PROCEDURE
- COMMON INP,IOUT,IPN
- DIMENSION A(NMAX,N),C(NMAX,N),AII(150),AJJ(150)
- SIN45 = .70710678
- COS45 = .70710678
- S45SQ = 0.50
- C45SQ = 0.50
-C UNIT EIGENVECTOR MATRIX
- DO 70 I = 1,N
- DO 7 J = I,N
- A(J,I)=A(I,J)
- C(I,J) = 0.0
- 7 C(J,I) = 0.0
- 70 C(I,I) = 1.0
-C DETERMINATION OF SEARCH ARGUMENT, TEST
- AMAX = 0.0
- DO 1 I = 1,N
- DO 1 J = 1,I
- TEMPA=DABS(A(I,J))
- IF (AMAX-TEMPA) 2,1,1
- 2 AMAX = TEMPA
- 1 CONTINUE
- TEST = AMAX*E
-C SEARCH FOR LARGEST OFF DIAGONAL ELEMENT
- DO 72 IJAC=1,MAXJAC
- AIJMAX = 0.0
- DO 3 I = 2,N
- LIM = I-1
- DO 3 J = 1,LIM
- TAIJ=DABS(A(I,J))
- IF (AIJMAX-TAIJ) 4,3,3
- 4 AIJMAX = TAIJ
- IPIV = I
- JPIV = J
- 3 CONTINUE
- IF(AIJMAX-TEST)300,300,5
-C PARAMETERS FOR ROTATION
- 5 TAII = A(IPIV,IPIV)
- TAJJ = A(JPIV,JPIV)
- TAIJ = A(IPIV,JPIV)
- TMT = TAII-TAJJ
- IF(DABS(TMT/TAIJ)-1.0D-12) 60,60,6
- 60 IF(TAIJ) 10,10,11
- 6 ZAMMA=TAIJ/(2.0*TMT)
- 90 IF(DABS(ZAMMA)-0.38268)8,8,9
- 9 IF(ZAMMA)10,10,11
- 10 SINT = -SIN45
- GO TO 12
- 11 SINT = SIN45
- 12 COST = COS45
- SINSQ = S45SQ
- COSSQ = C45SQ
- GO TO 120
- 8 GAMSQ=ZAMMA*ZAMMA
- SINT=2.0*ZAMMA/(1.0+GAMSQ)
- COST = (1.0-GAMSQ)/(1.0+GAMSQ)
- SINSQ=SINT*SINT
- COSSQ=COST*COST
-C ROTATION
- 120 DO 13 K = 1,N
- TAIK = A(IPIV,K)
- TAJK = A(JPIV,K)
- A(IPIV,K) = TAIK*COST+TAJK*SINT
- A(JPIV,K) = TAJK*COST-TAIK*SINT
- TCIK = C(IPIV,K)
- TCJK = C(JPIV,K)
- C(IPIV,K) = TCIK*COST+TCJK*SINT
- 13 C(JPIV,K) = TCJK*COST-TCIK*SINT
- A(IPIV,IPIV) = TAII*COSSQ+TAJJ*SINSQ+2.0*TAIJ*SINT*COST
- A(JPIV,JPIV) = TAII*SINSQ+TAJJ*COSSQ-2.0*TAIJ*SINT*COST
- A(IPIV,JPIV) = TAIJ*(COSSQ-SINSQ)-SINT*COST*TMT
- A(JPIV,IPIV) = A(IPIV,JPIV)
- DO 30 K = 1,N
- A(K,IPIV) = A(IPIV,K)
- 30 A(K,JPIV) = A(JPIV,K)
- 72 CONTINUE
- WRITE (IOUT,1000) AIJMAX
- 1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',
- 1 'MENT = ',1PE14.7)
-C ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER
- 300 DO 14 I=1,N
- 14 AJJ(I)=A(I,I)
- LT=N+1
- DO15 L=1,N
- LT=LT-1
- AIIMIN=1.0E+30
- DO16 I=1,N
- IF(AJJ(I)-AIIMIN)17,16,16
- 17 AIIMIN=AJJ(I)
- IT=I
- 16 CONTINUE
- IN=L
- AII(IN)=AIIMIN
- AJJ(IT)=1.0E+30
- DO15 K=1,N
- 15 A(IN,K)=C(IT,K)
- DO 18 I=1,N
- IF(A(I,1))19,22,22
- 19 T=-1.0
- GO TO 91
- 22 T=1.0
- 91 DO 18 J=1,N
- 18 C(J,I)=T*A(I,J)
- RETURN
- END