--- /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