      SUBROUTINE JACOBI(A,N,NP,D,V,IFAIL,NROT)                          JACO001
C ----------------------------------------------------------------------JACO002
C                                                                       JACO003
C     THIS SUBROUTINE COMPUTES ALL EIGENVALUES AND VECTORS OF A REAL    JACO004
C     SYMMETRIC MATRIX A USING THE JACOBI METHOD (Taken from Numerical  JACO005
C     Recipes - Press et al.)                                           JACO006
C                                                                       JACO007
C     PARAMETERS:                                                       JACO008
C                                                                       JACO009
C     A     = REAL SYMMETRIC MATRIX, WHICH IS OF SIZE N BY N, STORED IN JACO010
C             A PHYSICAL NP BY NP ARRAY. ON EXIT, ELEMENTS OF A ABOVE   JACO011
C             THE DIAGONAL ARE DESTROYED.                               JACO012
C     N     = ON ENTRY, ORDER OF MATRIX A FOR WHICH COMPUTATIONS ARE    JACO013
C             DONE. UNCHANGED ON EXIT.                                  JACO014
C     NP    = PHYSICAL ARRAY SUBSCRIPTS OF MATRIX A. UNCHANGED ON EXIT. JACO015
C     D     = REAL ARRAY OF PHYSICAL DIMENSION NP WHICH, ON EXIT,       JACO016
C             CONTAINS EIGENVALUES OF A IN ITS FIRST N ELEMENTS.        JACO017
C     V     = REAL MATRIX WITH THE SAME LOGICAL AND PHYSICAL DIMENSIONS JACO018
C             AS A WHOSE COLUMNS CONTAIN, ON EXIT, THE NORMALIZED       JACO019
C             EIGENVECTORS OF A.                                        JACO020
C     IFAIL = FAILURE INDICATOR. THIS IS SET TO 0 FOR NORMAL RETURN, ANDJACO021
C             TO 1 WHEN MORE THAN 50 * N(N-1)/2 ROTATIONS ARE NEEDED.   JACO022
C     NROT  = COUNTER WHICH, ON EXIT, CONTAINS THE NUMBER OF ROTATIONS. JACO023
C                                                                       JACO024
C     DURING COMPUTATION:                                               JACO025
C                                                                       JACO026
C     Z    = VECTOR WHICH ACCUMULATES TERMS OF FORM TA(PQ)              JACO027
C     B    = VECTOR WHICH CONTAINS CURRENT EIGENVALUES OF A             JACO028
C                                                                       JACO029
C     OPERATION COUNT:                                                  JACO030
C                                                                       JACO031
C     TYPICAL MATRICES REQUIRE 6 TO 10 SWEEPS OF N*(N-1)/2 JACOBI       JACO032
C     ROTATIONS TO ACHIEVE CONVERGENCE, OR 3*N**2 TO 5*N**2             JACO033
C     ROTATIONS. EACH ROTATION REQUIRES OF ORDER 6*N OPERATIONS, SO     JACO034
C     THE TOTAL LABOR IS OF THE ORDER 18*N**3 TO 30*N**3.               JACO035
C                                                                       JACO036
C ----------------------------------------------------------------------JACO037
      IMPLICIT REAL    (A-H,O-Z)                                        JACO038
      IMPLICIT INTEGER (I-N)                                            JACO039
      PARAMETER        (EMIN=1.0E-8)                                    JACO040
      PARAMETER        (NMAX=100)                                       JACO041
      LOGICAL          IFAIL                                            JACO042
      REAL             A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)          JACO043
      IFAIL=.FALSE.                                                     JACO044
                                                                        JACO045
C                                                                       JACO046
C   JACOBI METHOD NOT MOST EFFICIENT METHOD FOR MATRICES OF ORDER       JACO047
C   GREATER THAN 10 ; MAXIMUM ALLOWED PRESET BY NMAX                    JACO048
C                                                                       JACO049
      IF (N.GT.NMAX) THEN                                               JACO050
          WRITE(6,10)                                                   JACO051
10        FORMAT(/' *** MAXIMUM NUMBER OF ELEMENTS EXCEEDED IN ',       JACO052
     &    'SUBROUTINE JACOBI; CHANGE NMAX - RUN TERMINATED')            JACO053
          STOP                                                          JACO054
      ENDIF                                                             JACO055
C                                                                       JACO056
C   INITIALISE V TO THE IDENTITY MATRIX                                 JACO057
C                                                                       JACO058
      DO 12 IP=1,N                                                      JACO059
          DO 11 IQ=1,N                                                  JACO060
              V(IP,IQ)=0.0                                              JACO061
11        CONTINUE                                                      JACO062
          V(IP,IP)=1.0                                                  JACO063
12    CONTINUE                                                          JACO064
C                                                                       JACO065
C   INITIALISE B AND D TO THE DIAGONAL OF A AND Z TO 0.0                JACO066
C                                                                       JACO067
      DO 13 IP=1,N                                                      JACO068
          B(IP)=A(IP,IP)                                                JACO069
          D(IP)=B(IP)                                                   JACO070
          Z(IP)=0.0                                                     JACO071
13    CONTINUE                                                          JACO072
C                                                                       JACO073
C   START OF SWEEPS                                                     JACO074
C                                                                       JACO075
      NROT=0                                                            JACO076
      DO 24 I=1,50                                                      JACO077
C                                                                       JACO078
C---SUM OFF-DIAGONAL ELEMENTS OF MATRIX A                               JACO079
C                                                                       JACO080
          SM=0.0                                                        JACO081
          DO 15 IP=1,N-1                                                JACO082
              DO 14 IQ=IP+1,N                                           JACO083
                  SM=SM+ABS(A(IP,IQ))                                   JACO084
 14           CONTINUE                                                  JACO085
 15       CONTINUE                                                      JACO086
C                                                                       JACO087
C---NORMAL RETURN ON QUADRATIC CONVERGENCE TO MACHINE UNDERFLOW         JACO088
C                                                                       JACO089
          IF (SM.EQ.0.0) RETURN                                         JACO090
C                                                                       JACO091
C---SET TRESHOLD                                                        JACO092
C                                                                       JACO093
C---ON THE FIRST THREE SWEEPS                                           JACO094
          IF (I.LT.4) THEN                                              JACO095
              TRESH=0.2*SM/N**2                                         JACO096
          ELSE                                                          JACO097
                                                                        JACO098
C---THEREAFTER                                                          JACO099
              TRESH=0.0                                                 JACO100
          ENDIF                                                         JACO101
                                                                        JACO102
C---PERFORM PQ ROTATIONS                                                JACO103
C                                                                       JACO104
          DO 22 IP=1,N-1                                                JACO105
              DO 21 IQ=IP+1,N                                           JACO106
                  G=100.0*ABS(A(IP,IQ))                                 JACO107
                                                                        JACO108
C---AFTER 4 SWEEPS, SKIP ROTATION IF THE OFF-DIAGONAL ELEMENT IS SMALL  JACO109
                  IF ( (I.GT.4) .AND. (ABS(D(IP))+G.EQ.ABS(D(IP)))      JACO110
     &                .AND. (ABS(D(IQ))+G.EQ.ABS(D(IQ))) ) THEN         JACO111
                      A(IP,IQ)=0.0                                      JACO112
                                                                        JACO113
C---OTHERWISE ROTATE DEPENDENT ON THE TRESHOLD                          JACO114
                  ELSE IF (ABS(A(IP,IQ)).GT.TRESH) THEN                 JACO115
                      H=D(IQ)-D(IP)                                     JACO116
                      IF (ABS(H).LT.EMIN) H = 0.0                       JACO117
                      IF (ABS(H)+G.EQ.ABS(H)) THEN                      JACO118
                                                                        JACO119
C   T=1/(2THETA)                                                        JACO120
                          T=A(IP,IQ)/H                                  JACO121
                      ELSE                                              JACO122
                                                                        JACO123
C   T=SIGN(THETA) / (ABS(THETA)+ SQRT(THETA**2 + 1))                    JACO124
                          THETA=0.5*H/A(IP,IQ)                          JACO125
                          T=1.0/(ABS(THETA)+SQRT(1.0+THETA**2))         JACO126
                          IF (THETA.LT.0.0) T=-T                        JACO127
                      ENDIF                                             JACO128
                      C=1.0/SQRT(1+T**2)                                JACO129
                      IF (ABS(C).LT.EMIN) C = 0.0                       JACO130
                      S=T*C                                             JACO131
                      IF (ABS(S).LT.EMIN) S = 0.0                       JACO132
                      TAU=S/(1.0+C)                                     JACO133
                      IF (ABS(TAU).LT.EMIN) TAU = 0.0                   JACO134
                      H=T*A(IP,IQ)                                      JACO135
                      IF (ABS(H).LT.EMIN) H = 0.0                       JACO136
                      Z(IP)=Z(IP)-H                                     JACO137
                      IF (ABS(Z(IP)).LT.EMIN) Z(IP) = 0.0               JACO138
                      Z(IQ)=Z(IQ)+H                                     JACO139
                      IF (ABS(Z(IQ)).LT.EMIN) Z(IQ) = 0.0               JACO140
                      D(IP)=D(IP)-H                                     JACO141
                      IF (ABS(D(IP)).LT.EMIN) D(IP) = 0.0               JACO142
                      D(IQ)=D(IQ)+H                                     JACO143
                      IF (ABS(D(IQ)).LT.EMIN) D(IQ) = 0.0               JACO144
                      A(IP,IQ)=0.0                                      JACO145
                                                                        JACO146
C   CASE OF ROTATIONS 1.LE.J.LT.P                                       JACO147
                      DO 16 J=1,IP-1                                    JACO148
                          G=A(J,IP)                                     JACO149
                          IF (ABS(G).LT.EMIN) G = 0.0                   JACO150
                          H=A(J,IQ)                                     JACO151
                          IF (ABS(H).LT.EMIN) H = 0.0                   JACO152
                          A(J,IP)=G-S*(H+G*TAU)                         JACO153
                          IF (ABS(A(J,IP)).LT.EMIN) A(J,IP) = 0.0       JACO154
                          A(J,IQ)=H+S*(G-H*TAU)                         JACO155
                          IF (ABS(A(J,IQ)).LT.EMIN) A(J,IQ) = 0.0       JACO156
16                    CONTINUE                                          JACO157
                                                                        JACO158
C   CASE OF ROTATIONS P.LT.J.LT.Q                                       JACO159
                      DO 17 J=IP+1,IQ-1                                 JACO160
                          G=A(IP,J)                                     JACO161
                          H=A(J,IQ)                                     JACO162
                          A(IP,J)=G-S*(H+G*TAU)                         JACO163
                          IF (ABS(A(IP,J)).LT.EMIN) A(IP,J) = 0.0       JACO164
                          A(J,IQ)=H+S*(G-H*TAU)                         JACO165
                          IF (ABS(A(J,IQ)).LT.EMIN) A(J,IQ) = 0.0       JACO166
17                    CONTINUE                                          JACO167
                                                                        JACO168
C   CASE OF ROTATIONS Q.LT.J.LE.N                                       JACO169
                      DO 18 J=IQ+1,N                                    JACO170
                          G=A(IP,J)                                     JACO171
                          H=A(IQ,J)                                     JACO172
                          A(IP,J)=G-S*(H+G*TAU)                         JACO173
                          IF (ABS(A(IP,J)).LT.EMIN) A(IP,J) = 0.0       JACO174
                          A(IQ,J)=H+S*(G-H*TAU)                         JACO175
                          IF (ABS(A(IQ,J)).LT.EMIN) A(IQ,J) = 0.0       JACO176
18                    CONTINUE                                          JACO177
                      DO 19 J=1,N                                       JACO178
                          G=V(J,IP)                                     JACO179
                          H=V(J,IQ)                                     JACO180
                          V(J,IP)=G-S*(H+G*TAU)                         JACO181
                          IF (ABS(V(J,IP)).LT.EMIN) V(J,IP) = 0.0       JACO182
                          V(J,IQ)=H+S*(G-H*TAU)                         JACO183
                          IF (ABS(V(J,IQ)).LT.EMIN) V(J,IQ) = 0.0       JACO184
19                    CONTINUE                                          JACO185
                      NROT=NROT+1                                       JACO186
                  ENDIF                                                 JACO187
21            CONTINUE                                                  JACO188
22        CONTINUE                                                      JACO189
                                                                        JACO190
C---END OF PQ ROTATIONS                                                 JACO191
C                                                                       JACO192
                                                                        JACO193
C---UPDATE D WITH THE SUM OF TA(PQ) AND REINITIALIZE Z TO 0.0           JACO194
C                                                                       JACO195
          DO 23 IP=1,N                                                  JACO196
              B(IP)=B(IP)+Z(IP)                                         JACO197
              D(IP)=B(IP)                                               JACO198
              Z(IP)=0.0                                                 JACO199
23        CONTINUE                                                      JACO200
24    CONTINUE                                                          JACO201
                                                                        JACO202
C   END OF SWEEPS                                                       JACO203
C                                                                       JACO204
C---ABNORMAL RETURN; 50 SWEEPS IN JACOBI SHOULD NEVER HAPPEN            JACO205
C                                                                       JACO206
      IFAIL=.TRUE.                                                      JACO207
      RETURN                                                            JACO208
      END                                                               JACO209
