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