C #JACOB   V1C   17-FEB-72
C LAST UPDATE 26-APR-73
C
      SUBROUTINE JACOB(A,U,N,M,EPS)
      DIMENSION A(M,M),U(M,M)
C
C *PURPOSE.
C       TO DIAGONALIZE A REAL SYMMETRIC MATRIX.
C
C *PARAMETERS.
C
C  A   - REAL ARRAY OF DIMENSION "M" X "M" WHICH CONTAINS THE
C 	SYMMETRIC REAL "N" X "N" MATRUX TO BE DIAGONALIZED IN
C 	ITS TOP LEFT CORNER.  ONLY THE UPPER TRIANGLE OF "A"
C 	IS IN FACT USED.  ON RETURN THE EIGENVALUES ARE ALONG
C 	THE DIAGONAL IN ALGEBARAICALLY INCREASING ORDER. THE
C 	OFF DIAGONAL ELEMENTS ARE REDUCED TO SMALL NUMBERS.
C  U   - REAL ARRAY OF DIMENSION "M" X "M" IN WHICH THE EIGENVECTORS
C 	ARE RETURNED.  THE K'TH EIGENVECTOR CORRESPONDING TO
C 	THE EIGENVALUE A(K,K) IS IN THE K'TH COLUMN U(*,K).
C  N   - THE LOGICAL SIZE OF THE MATRIX TO BE DIAGONALIZED.
C  M   - THE PHYSICAL SIZE OF THE ARRAYS CONTAINING THE MATRICES.
C  EPS - A TOLERANCE PARAMETER CONTROLLING THE CUTOFF POINT IN
C 	THE JACOBI ITERATIVE PROCESS.  TYPICALLY "EPS" IS
C 	ABOUT 1 OR 2 ORDERS OF MAGNITUDE GREATER THAN THE MACHINE
C 	SIGNIFICANCE.
C
C INPUT PARAMETERS:
C 	A, N, M, EPS
C OUTPUT PARAMETERS:
C 	A, U
C
C *METHOD.
C       THE STANDARD JACOBI ITERATIVE PROCESS IS USED, BUT THE THRESHOLD
C IS REDUCED AFTER EACH ROTATION AND MAINTAINED AT A LITTLE MORE THAN
C THE R.M.S. VALUE OF THE OFF DIAGONAL ELEMENTS.  THIS IS DESIGNED TO
C REDUCE THE NUMBER OF ROTATIONS TO A MINIMUM.
C
C *ACCURACY.
C       DEPENDS ON EPS AND THE MACHINE SIGNIFICANCE,  HOWEVER IF EPS
C IS NOT TOO SMALL THE ROUNDING ERRORS CAN BE KEPT TO ONE ORDER
C OF MAGNITUDE GREATER THEN THE SIGNIFICANCE.
C
C *RESTRICTIONS.
C
C *ERROR CONDITIONS.
C
C *NON STANDARD ROUTINES CALLED.
C
C *TYPICAL TIMES.
C       THE TIME IS ROUGHLY PROPORTIONAL TO N**3.
C
C *ORIGIN.
C 	THIS ROUTINE IS TAKEN FROM ONE BY J.P.M.BAILEY WHO TOOK
C IT FROM ONE BY SOMEONE ELSE.  IT HAS NOW BEEN REFINED TO THE N'TH
C DEGREE AND IS PROBABLY ONE OF THE MOST EFFICIENT
C JACOBI ROUTINES AROUND.  SOME OF THE DEVELOPMENT DONE BY
C M.R.MANNING AND J.R.HULSTON.
C
C *COMMENTS.
C 	FOR MATRICES THAT ARE LARGE (SAY N>10) AND NOT SPARSE
C THERE CAN BE LITTLE JUSTIFICATION FOR USING "JACOB" IN
C PREFERENCE TO "QREGB" OR A SIMILAR ROUTINE.
C
C
C #END.
C
C
      DATA  EPS2 / 0.5E-4/
C  EPS2 IS USED TO SET DSTOP - SEE BELOW.
C #START.
C
C  INITIALIZE U ARRAY AND GET MAX ELEMENT OF A
C
      NM1=N-1
      IF(NM1) 800, 800, 5
5     ROOT2=SQRT(0.5)
      DO 40 J=1,N
      DO 20 I=1,N
20    IF( I .NE. J)  U(I,J)=0.
40    U(J,J)=1.
C
      IF(EPS) 800, 800, 60
60    AVGF=FLOAT(N*NM1)*0.55
      D=0.
      DSTOP=0.
C
C  START A SWEEP
C
100   IFLAG=0
      DO 400 JCOL=2,N
      JCOL1=JCOL-1
      J1=JCOL+1
      DO 400 IROW=1,JCOL1
C
C  CHECK SUM OF SQUARES OF OFF-DIAGONAL ELEMENTS TO SEE IF
C  IT NEEDS RECALCULATING.
      IF( D.GT.DSTOP) GO TO 150
C   RECALCULATE D WHEN IT GETS TOO SMALL BECAUSE OF ROUNDING ERRORS.
C
      D=0.
      DO 130 J=2, N
      JM1= J - 1
      DO 130 I= 1, JM1
      S= A(I, J)
130   D= D + S*S
      DSTOP= EPS2*D
C  AND RESET THRESHOLD NOW.
      THRSH= SQRT(D/AVGF)
C
150   AIJ=A(IROW,JCOL)
C
C  COMPARE OFF DIAGONAL ELEMENT WITH THRESHOLD - NOTE THAT
C  THIS IS CONSTRUCTED SO THERE SHOULD ALWAYS BE AT LEAST
C  ONE OFF DIAGONAL ELEMENT LARGER!!
C
      IF(ABS(AIJ) .LE. THRSH )  GO TO 400
      AJJ=A(JCOL,JCOL)
      AII=A(IROW,IROW)
      V=AJJ-AII
C  CHECK IF CLOSE TO 45 DEGREES.
      IF(ABS(V) .GT. EPS*ABS(AIJ)) GO TO 200
      S=ROOT2
      C=ROOT2
      GO TO 220
C
C  NORMAL CALCULATION OF SIN AND COS
C
200    T = AIJ/V
C
C   CHECK WHETHER ROTATION SIGNIFICANT.
C  THIS IS THE ONE AND ONLY TEST FOR CONVERGENCE!!!!
C
      IF(ABS(T) .LT. EPS) GO TO 400
      S=0.25/SQRT(0.25 + T*T)
      C=SQRT(0.5+S)
      S=(T+T)*S/C
C
C  ROTATION APPLIED TO BOTH SIDES OF A.
C
220    IFLAG = 1
      I1=IROW-1
      IF(I1) 250,250,230
230   DO 240 I=1,I1
      T=A(I,IROW)
      U=A(I,JCOL)
      A(I,IROW)=C*T-S*U
240   A(I,JCOL)=S*T+C*U
250   I1=IROW+1
      IF(I1.GT.JCOL1)  GO TO 270
      DO 260 I=I1,JCOL1
      T=A(I,JCOL)
      U=A(IROW,I)
      A(IROW,I)=C*U-S*T
260   A(I,JCOL)=S*U+C*T
270   IF(J1.GT.N) GO TO 290
      DO 280 J=J1,N
      T=A(IROW,J)
      U=A(JCOL,J)
      A(IROW,J)=C*T-S*U
280   A(JCOL,J)=S*T+C*U
290   T=C*C*AII + S*( S*AJJ - (C+C)*AIJ)
      A(JCOL,JCOL)=AII+AJJ-T
      A(IROW,IROW)=T
      A(IROW,JCOL)=(C*C - S*S)*AIJ - S*C*V
C
C   ROTATION NOW APPLED TO ONE SIDE OF U.
C
      DO 310 I=1,N
      T=U(I,IROW)
      U=U(I,JCOL)
      U(I,IROW)=C*T-U*S
310   U(I,JCOL)=S*T+U*C
C
C  CALCULATE NEW SUM OF SQUARES OF OFF DIAGONALS
C
      D=D- AIJ*AIJ
C  SET NEW THRESHOLD AFTER REMOVING EACH OFF DIAGONAL ELEMENT.
C
      IF( D .GE. 0.)  THRSH= SQRT( D/AVGF )
400   CONTINUE
C
      IF(IFLAG) 100, 500, 100
C
C  CONVERGED NOW ORDER EIGENVALUES
C
500   DO 600 I= 1, NM1
C  FIND THE I'TH SMALLEST EIGENVALUE.
      K= I
      T= A(I,I)
      I1=I + 1
      DO 550 J= I1, N
      IF ( A(J,J) .GE. T)  GO TO 550
      K=J
      T= A(J,J)
550   CONTINUE
      IF ( K .EQ. I) GO TO 600
C  IF NOT IN I'TH COLUMN PUT IT THERE.
      A(K,K)= A(I,I)
      A(I,I)= T
      DO 580 J=1, N
      T= U(J, I)
      U(J, I)= U(J, K)
580   U(J, K)= T
600   CONTINUE
C
800   RETURN
      END
���������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������