/* THE IMPLICIT Q R ALGORITHM */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 243-244. */ IMTQL1:PROC(N,MACHEPS,D,E,FAIL) OPTIONS (REORDER); /*PAGE 243*/ /*00000010*/ DCL (N,I,J,K,M,ITS,UBN,LBN) BIN FIXED (31,0), /*00000020*/ FAIL LABEL, /*00000030*/ (MACHEPS,H,C,P,Q,S,T,U) BIN FLOAT(21), /*00000040*/ (D,E)(*) BIN FLOAT(21); /*00000050*/ LBN=LBOUND(D,1); UBN=LBN+N-1; /*00000060*/ DO I=LBN+1 TO UBN; /*00000070*/ E(I-1)=E(I); /*00000080*/ END; /*00000090*/ E(UBN)=0E0; K=UBN-1; /*00000100*/ DO J=LBN TO UBN; /*00000110*/ ITS=0; /*00000120*/ /* LOOK FOR SINGLE SMALL SUBDIAGONAL ELEMENT */ /*00000130*/ TEST: /*00000140*/ DO M=J TO K; /*00000150*/ IF ABS(E(M))<=MACHEPS*(ABS(D(M))+ABS(D(M+1)))THEN GO TO CONT1;/*00000160*/ END; /*00000170*/ M=N; /*00000180*/ CONT1: /*00000190*/ U=D(J); /*00000200*/ IF M^=J THEN DO; /*00000210*/ IF ITS=30 THEN GO TO FAIL; /*00000220*/ ITS=ITS+1; /*00000230*/ /*FORM SHIFT */ /*00000240*/ Q=(D(J+1)-U)/2E0*E(J) ; T=SQRT(1E0+Q**2); /*00000250*/ IF Q<0E0 THEN C=Q-T; ELSE C=Q+T; /*00000260*/ Q=D(M)-U+E(J)/C; U=0E0; S,C=1E0; /*00000270*/ DO I=M-1 TO J BY -1; /*00000280*/ P=S*E(I); H=C*E(I); /*00000290*/ IF ABS(P)>=ABS(Q) THEN DO; /*00000300*/ C=Q/P; /*00000310*/ T=SQRT(C**2+1E0); /*00000320*/ E(I+1)=P*T; /*00000330*/ S=1E0/T; C=C*S; /*00000340*/ END; /*00000350*/ ELSE DO; /*00000360*/ S=P/Q; /*00000370*/ T=SQRT(S**2+1E0); /*00000380*/ E(I+1)=Q*T; /*00000390*/ C=1E0/T; S=S*C; /*00000400*/ END; /*00000410*/ Q=D(I+1)-U; T=(D(I)-Q)*S+2E0*C*H; /*00000420*/ U=S*T; D(I+1)=Q+U; Q=C*T-H; /*00000430*/ END; /*00000440*/ D(J)=D(J)-U; E(J)=Q; E(M)=0E0; GO TO TEST; /*00000450*/ END; /*00000460*/ /*ORDER EIGENVALUE */ /*00000470*/ DO I=J TO LBN+1 BY -1; /*00000480*/ IF U