%PROCESS MARGINS(1,100); /* THE CALCULATION OF SPECIFIED EIGENVECTORS BY INVERSE ITERATION */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp.423-427. */ /* p. 423 */ /* In two places, IF INT(I) THEN .. was changed to IF INT(I) ^= '0'B THEN .. */ /* 22/3/2007 */ TRISTURM: /* p. 423 */ PROC(N,LB,UB,MACHEPS,EPS1,C,B,BETA,M,ROOT,VEC,COUNT,FAIL) OPTIONS (REORDER); DCL (N,M,I,P,Q,R,M1,M2,J,K,S,ITS,GROUP,UBC,LBC) BIN FIXED(31,0), COUNT(*) BIN FIXED(15,0), (LB,UB,MACHEPS,EPS1,NORM,X0,X1,XU,U,V,BI,EPS2,EPS3,EPS4, C(*),B(*),BETA(*),ROOT(*),VEC(*,*)) BIN FLOAT(21), INT(P:Q) BIT(8) CONTROLLED, (X,WU)(M1:M2) BIN FLOAT(21) CONTROLLED,FAIL LABEL, (D,E,F,Y,Z)(P:Q) BIN FLOAT(21) CONTROLLED; LBC=LBOUND(C,1);UBC=LBC+N-1; STURMCNT:PROC(P,Q,D,E,F,LAMBDA) RETURNS(BIN FIXED(31,0)); DCL (P,Q) BIN FIXED(31,0), LAMBDA BIN FLOAT(21), (D,E,F)(*) BIN FLOAT(21), X BIN FLOAT(21), (COUNT,I) BIN FIXED(31,0); COUNT=0; X=1.; DO I=P TO Q; IF X^=0. THEN X=F(I)/X; ELSE X=ABS(E(I)/MACHEPS); X=D(I)-LAMBDA-X; IF X<0. THEN COUNT=COUNT+1; END; RETURN(COUNT); END STURMCNT; IF EPS1<=0. THEN EPS1=(ABS(LB)+ABS(UB))*MACHEPS; /* LOOK FOR SMALL SUBDIAG ENTRIES */ DO I=LBC+1 TO UBC; IF ABS(B(I))<=MACHEPS*(ABS(C(I))+ABS(C(I-1))) THEN B(I),BETA(I)=0E0; END; R=M; M=STURMCNT(LBC,UBC,C,B,BETA,UB)-STURMCNT(LBC,UBC,C,B,BETA,LB); IF M>R THEN GO TO FAIL; Q=LBC-1; R=LBOUND(ROOT,1); NEXTP: P=Q+1; DO Q=P TO UBC-1; IF ABS(B(Q+1))=0E0 THEN GO TO SUB; END; Q=UBC; SUB: IF P=Q THEN DO; IF LB<=C(P)&C(P)2E0*MACHEPS*(ABS(XU)+ABS(X0))+EPS1); S=STURMCNT(P,Q,C,B,BETA,X1); IF S=ABS(U)); IF INT(I) ^= '0'B THEN DO; Y(I),XU=U/BI; D(I-1)=BI; E(I-1)=C(I)-X1; IF I^=Q THEN F(I-1)=B(I+1); ELSE F(I-1)=0E0; U=V-XU*E(I-1); V=-XU*F(I-1); END; ELSE DO; Y(I),XU=BI/U; D(I-1)=U; E(I-1)=V; F(I-1)=0E0; U=C(I)-X1-XU*V; IF I^=Q THEN V=B(I+1); END; END /* I LOOP */; IF U^=0E0 THEN D(Q)=U; ELSE D(Q)=EPS3; E(Q),F(Q)=0E0; /* BACKSUBSTITUTION*/ NEWZ: DO I=Q BY -1 TO P; Z(I)=(Z(I)-U*E(I)-V*F(I))/D(I); V=U; U=Z(I); END /*BACKSUB*/; DO J=R-GROUP TO R-1; /* ORTHOGONALISE WRT PREVIOUS MEMBERS OF GROUP */ XU=0E0; DO I=P TO Q; XU=XU+Z(I)*VEC(I,J); END; DO I=P TO Q; Z(I)=Z(I)-XU*VEC(I,J); END; END /*ORTHOGONALISE */; NORM=0E0; DO I=P TO Q; NORM=NORM+ABS(Z(I)); END; /* FORWARD SUBSTITUTION */ IF NORM<1E0 THEN DO; IF ITS=5 THEN DO; COUNT(R)=6; GO TO FAIL; END; IF NORM=0E0 THEN DO; Z(S)=EPS4; IF S^=Q THEN S=S+1; ELSE S=P; END /* NULL VECTOR */; ELSE DO; XU=EPS4/NORM; DO I=P TO Q; Z(I)=Z(I)*XU; END; END; DO I=P+1 TO Q; IF INT(I) ^= '0'B THEN DO; U=Z(I-1); Z(I-1)=Z(I); Z(I)=U-Y(I)*Z(I); END; ELSE Z(I)=Z(I)-Y(I)*Z(I-1); END; ITS=ITS+1; GO TO NEWZ; END /*FORWARDSUB */; /*NORMALISE SO THAT SUM OF SQRS IS 1 AND EXPAND TO FULL ORDER */ U=0E0; DO I=P TO Q; U=U+Z(I)**2; END; XU=1E0/SQRT(U); DO I=P TO Q; VEC(I,R)=Z(I)*XU; END; DO I=P-1 BY -1 TO LBC, Q+1 TO UBC; VEC(I,R)=0E0; END; COUNT(R)=ITS; R=R+1; X0=X1; END /*K-TH VECTOR */; END /* M1<=M2 */; END /*P^=Q*/; IF Q