/* THE JACOBI METHOD FOR REAL SYMMETRIC MATRICES */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 204-206. */ JACOBI:PROC(N,EIVEC,A,D,V,ROT) OPTIONS (REORDER); /*PAGE 204*/ DCL (N,ROT,P,Q,I,J,UBN,LBN) BIN FIXED (31,0), EIVEC BIT(1), (A(*,*), D(*), V(*,*), (B,Z)(LBN:UBN) CTL) BIN FLOAT(21), (SM,C,S,T,H,G,TAU,THETA,TRESH) FLOAT BIN(21); LBN=LBOUND(A,1); UBN=LBN+N-1; ALLOC B,Z; PROGRAM: IF EIVEC THEN DO P=LBN TO UBN; V(P,P)=1E0; DO Q=LBN TO P-1,P+1 TO UBN; V(P,Q)=0E0; END; END; DO P=LBN TO UBN; B(P),D(P)=A(P,P); Z(P)=0E0; END; ROT=0; DO I=1 TO 50; SWP: SM=0E0; DO P=LBN TO UBN-1; DO Q=P+1 TO UBN; SM=SM+ABS(A(P,Q)); END; END; IF SM=0E0 THEN GO TO OUT; IF I<4 THEN TRESH=.2E0*SM/N**2; ELSE TRESH=0E0; DO P=LBN TO UBN-1; DO Q=P+1 TO UBN; G=1E2*ABS(A(P,Q)); IF I>4&(ABS(D(P))+G=ABS(D(P))) &(ABS(D(Q))+G=ABS(D(Q))) THEN A(P,Q)=0E0; ELSE IF ABS(A(P,Q))>TRESH THEN ROTATE: DO; H=D(Q)-D(P); IF ABS(H)+G=ABS(H) THEN T=A(P,Q)/H; ELSE DO; THETA=.5E0*H/A(P,Q); T=1E0/(ABS(THETA)+ SQRT(1E0+THETA**2)); IF THETA<0E0 THEN T=-T; END; /*END COMPUTING TAN OF ROTATION ANGLE*/ C=1E0/SQRT(1E0+T**2); S=T*C; TAU=S/(1E0+C); H=T*A(P,Q); Z(P)=Z(P)-H; Z(Q)=Z(Q)+H; D(P)=D(P)-H; D(Q)=D(Q)+H; A(P,Q)=0E0; DO J=LBN TO P-1; G=A(J,P); H=A(J,Q); A(J,P)=G-S*(H+G*TAU); A(J,Q)=H+S*(G-H*TAU); END;/*END OF CASE 1<=J