source/unres/src_MD/djacob.f

   1       SUBROUTINE DJACOB(N,NMAX,MAXJAC,E,A,C,AII)
   2       IMPLICIT REAL*8 (A-H,O-Z)
   3 C     THE JACOBI DIAGONALIZATION PROCEDURE
   4       COMMON INP,IOUT,IPN
   5       DIMENSION A(NMAX,N),C(NMAX,N),AII(150),AJJ(150)
   6       SIN45 = .70710678
   7       COS45 = .70710678
   8       S45SQ = 0.50
   9       C45SQ = 0.50
  10 C     UNIT EIGENVECTOR MATRIX
  11       DO 70 I = 1,N
  12       DO 7 J = I,N
  13       A(J,I)=A(I,J)
  14       C(I,J) = 0.0
  15     7 C(J,I) = 0.0
  16    70  C(I,I) = 1.0
  17 C     DETERMINATION OF SEARCH ARGUMENT, TEST
  18       AMAX = 0.0
  19       DO 1 I = 1,N
  20       DO 1 J = 1,I
  21       TEMPA=DABS(A(I,J))
  22       IF (AMAX-TEMPA) 2,1,1
  23     2 AMAX = TEMPA
  24     1 CONTINUE
  25       TEST = AMAX*E
  26 C     SEARCH FOR LARGEST OFF DIAGONAL ELEMENT
  27       DO 72 IJAC=1,MAXJAC
  28       AIJMAX = 0.0
  29       DO 3 I = 2,N
  30       LIM = I-1
  31       DO 3 J = 1,LIM
  32       TAIJ=DABS(A(I,J))
  33        IF (AIJMAX-TAIJ) 4,3,3
  34     4 AIJMAX = TAIJ
  35       IPIV = I
  36       JPIV = J
  37     3 CONTINUE
  38       IF(AIJMAX-TEST)300,300,5
  39 C     PARAMETERS FOR ROTATION
  40     5 TAII = A(IPIV,IPIV)
  41       TAJJ = A(JPIV,JPIV)
  42       TAIJ = A(IPIV,JPIV)
  43       TMT = TAII-TAJJ
  44       IF(DABS(TMT/TAIJ)-1.0D-12) 60,60,6
  45    60 IF(TAIJ) 10,10,11
  46     6 ZAMMA=TAIJ/(2.0*TMT)
  47    90 IF(DABS(ZAMMA)-0.38268)8,8,9
  48     9 IF(ZAMMA)10,10,11
  49    10 SINT = -SIN45
  50       GO TO 12
  51    11 SINT = SIN45
  52    12 COST = COS45
  53       SINSQ = S45SQ
  54       COSSQ = C45SQ
  55       GO TO 120
  56     8 GAMSQ=ZAMMA*ZAMMA
  57       SINT=2.0*ZAMMA/(1.0+GAMSQ)
  58       COST = (1.0-GAMSQ)/(1.0+GAMSQ)
  59       SINSQ=SINT*SINT
  60       COSSQ=COST*COST
  61 C     ROTATION
  62   120 DO 13 K = 1,N
  63       TAIK = A(IPIV,K)
  64       TAJK = A(JPIV,K)
  65       A(IPIV,K) = TAIK*COST+TAJK*SINT
  66       A(JPIV,K) = TAJK*COST-TAIK*SINT
  67       TCIK = C(IPIV,K)
  68       TCJK = C(JPIV,K)
  69       C(IPIV,K) = TCIK*COST+TCJK*SINT
  70    13 C(JPIV,K) = TCJK*COST-TCIK*SINT
  71       A(IPIV,IPIV) = TAII*COSSQ+TAJJ*SINSQ+2.0*TAIJ*SINT*COST
  72       A(JPIV,JPIV) = TAII*SINSQ+TAJJ*COSSQ-2.0*TAIJ*SINT*COST
  73       A(IPIV,JPIV) = TAIJ*(COSSQ-SINSQ)-SINT*COST*TMT
  74       A(JPIV,IPIV) = A(IPIV,JPIV)
  75       DO 30 K = 1,N
  76       A(K,IPIV) = A(IPIV,K)
  77    30 A(K,JPIV) = A(JPIV,K)
  78    72 CONTINUE
  79       WRITE (IOUT,1000) AIJMAX
  80  1000 FORMAT (/1X,'NONCONVERGENT JACOBI. LARGEST OFF-DIAGONAL ELE',
  81      1 'MENT = ',1PE14.7)
  82 C     ARRANGEMENT OF EIGENVALUES IN ASCENDING ORDER
  83   300 DO 14 I=1,N
  84    14 AJJ(I)=A(I,I)
  85       LT=N+1
  86       DO15 L=1,N
  87       LT=LT-1
  88       AIIMIN=1.0E+30
  89       DO16 I=1,N
  90       IF(AJJ(I)-AIIMIN)17,16,16
  91    17 AIIMIN=AJJ(I)
  92       IT=I
  93    16 CONTINUE
  94       IN=L
  95       AII(IN)=AIIMIN
  96       AJJ(IT)=1.0E+30
  97       DO15 K=1,N
  98    15 A(IN,K)=C(IT,K)
  99       DO 18 I=1,N
 100       IF(A(I,1))19,22,22
 101    19 T=-1.0
 102       GO TO 91
 103    22 T=1.0
 104    91 DO 18 J=1,N
 105    18 C(J,I)=T*A(I,J)
 106       RETURN
 107       END