/* EIGENVECTORS OF REAL AND COMPLEX MATRICES */ /* From J. H. Wilkinson & C. Reinsch, "Handbook for Automatic Computation", Vol. II, "Linear Algebra", Springer-Verlag, New York, 1971, pp. 376-388. */ /* P. 381 */ DIRTRANS:PROC(N,LOW,UPP,INT,H,V) OPTIONS (REORDER); /*PAGE 381*/ DCL (N,LOW,UPP,I,J,M,LBN,UBN) BIN FIXED(31,0), INT(*) BIN FIXED(15,0), (V,H)(*,*) BIN FLOAT(21), X BIN FLOAT (21); LBN=LBOUND(H,1); UBN=LBN+N-1; DO I=LBN TO UBN; V(I,I)=1E0; DO J=LBN TO I-1,I+1 TO UBN; V(I,J)=0E0; END; END; DO I=UPP TO LOW+1 BY -1; DO J=LOW+1 TO I-1; V(I,J)=H(I,J-1); END; M=INT(I); IF M^=I THEN DO J=LOW+1 TO UPP; X =V(M,J); V(M,J)=V(I,J); V(I,J)=X; END; END; END DIRTRANS; /* P. 382 */ ORTRANS:PROC(N,LOW,UPP,H,D,V) OPTIONS (REORDER); /*PAGE 382*/ /*00000010*/ DCL (N,LOW,UPP,I,J,K,M,LBN,UBN) BIN FIXED(31,0), /*00000020*/ X BIN FLOAT(53), /*00000030*/ Y BIN FLOAT(21), /*00000040*/ (H,V)(*,*) FLOAT BIN (21), /*00000050*/ D(*) FLOAT BIN(21); /*00000060*/ LBN=LBOUND(H,1); UBN=LBN+N-1; /*00000070*/ DO I=LBN TO UBN; /*00000080*/ V(I,I)=1E0; /*00000090*/ DO J=LBN TO I-1,I+1 TO UBN; /*00000100*/ V(I,J)=0E0; /*00000110*/ END; /*00000120*/ END; /*00000130*/ DO K=UPP-2 TO LOW BY -1; /*00000140*/ M=K+1; Y=H(M,K); /*00000150*/ IF Y=0E0 THEN GO TO SKIP; /*00000160*/ Y=Y*D(M); /*00000170*/ DO I=K+2 TO UPP; /*00000180*/ D(I)=H(I,K); /*00000190*/ END; /*00000200*/ DO J=M TO UPP; /*00000210*/ X=0E0; /*00000220*/ DO I=M TO UPP; /*00000230*/ X=X+D(I)*V(I,J); /*00000240*/ END; /*00000250*/ X=X/Y; /*00000260*/ DO I=M TO UPP; /*00000270*/ V(I,J)=V(I,J)+X*D(I); /*00000280*/ END; /*00000290*/ END; /*00000300*/ SKIP: /*00000310*/ END; /*00000320*/ END ORTRANS; /*00000330*/ /* P. 383 */ HQR2:PROC(N,LOW,UPP,MACHEPS,H,VECS,W,CNT,FAIL) OPTIONS (REORDER); /*PAGE 383*/ DCL (N,LOW,UPP,I,J,K,L,M,NA,ITS,EN,LBN,UBN) BIN FIXED(31), /*00000020*/ (NORM,P,R) BIN FLOAT(53),(V,TR) BIN FLOAT (21) COMPLEX, /*00000030*/ (Q, S,T,WT,X,Y,Z,RA,SA,U,MACHEPS) BIN FLOAT(21), /*00000040*/ NOTLAST BIT(1), /*00000050*/ CNT(*) BIN FIXED(15,0), /*00000060*/ FAIL LABEL, /*00000070*/ H(*,*) BIN FLOAT(21), /*00000080*/ W(*) BIN FLOAT(21) COMPLEX, /*00000090*/ VECS(*,*) BIN FLOAT(21); /*00000100*/ LBN=LBOUND(H,1); UBN=LBN+N-1; /*00000110*/ DO I=LBN TO LOW-1,UPP+1 TO UBN; /*00000120*/ W(I)=H(I,I); CNT(I)=0; /*00000130*/ END; /*00000140*/ EN=UPP; T=0E0; /*00000150*/ NEXTW: /*00000160*/ IF EN0E0 THEN DO; /*00001110*/ /*REAL PAIR */ /*00001120*/ IF P<0E0 THEN Z=P-Z; ELSE Z=P+Z; /*00001130*/ W(NA)=X+Z; W(EN),S=X-WT/Z; /*00001140*/ X=H(EN,NA); R=SQRT(X**2+Z**2); P=X/R; Q=Z/R; /*00001150*/ DO J=NA TO UBN; /*00001160*/ Z=H(NA,J); H(NA,J)=Q*Z+P*H(EN,J); /*00001170*/ H(EN,J)=Q*H(EN,J)-P*Z; /*00001180*/ END;/*END ROW MODIFICATION */ /*00001190*/ DO I=LBN TO EN; /*00001200*/ Z=H(I,NA); H(I,NA)=Q*Z+P*H(I,EN); /*00001210*/ H(I,EN)=Q*H(I,EN)-P*Z; /*00001220*/ END;/*END COLUMN MODIFICATION */ /*00001230*/ DO I=LOW TO UPP; /*00001240*/ Z=VECS(I,NA); VECS(I,NA)=Q*Z+P*VECS(I,EN); /*00001250*/ VECS(I,EN)=Q*VECS(I,EN)-P*Z; /*00001260*/ END;/*END ACCUMULATE */ /*00001270*/ END;/*END PAIR OF REAL ROOTS */ /*00001280*/ ELSE DO; /*00001290*/ /*COMPLEX PAIRS*/ /*00001300*/ REAL(W(NA)),REAL(W(EN))=X+P; /*00001310*/ IMAG(W(NA))=Z; IMAG(W(EN))=-Z; /*00001320*/ END; /*00001330*/ /*END 2 ROOTS FOUND */ /*00001340*/ EN=EN-2; GO TO NEXTW; /*00001350*/ /*ALL ROOTS FOUND, NOW BACKSUBSTITUTE */ /*00001360*/ FIN: /*00001370*/ NORM=0E0; K=LBN; /*00001380*/ DO I=LBN TO UBN; /*00001390*/ DO J=K TO UBN; /*00001400*/ NORM=NORM+ABS(H(I,J)); /*00001410*/ END; /*00001420*/ K=I; /*00001430*/ END; /*00001440*/ /*BACKSUBSTITUTION */ /*00001450*/ DO EN=UBN TO LBN BY -1; /*00001460*/ P=REAL(W(EN)); Q=IMAG(W(EN)); NA=EN-1; /*00001470*/ IF Q=0E0 THEN DO; /*00001480*/ /* REAL VECTOR */ /*00001490*/ M=EN; H(EN,EN)=1E0; /*00001500*/ DO I=NA TO 1 BY -1; /*00001510*/ WT=H(I,I)-P; R=H(I,EN); /*00001520*/ DO J=M TO NA; /*00001530*/ R=R+H(I,J)*H(J,EN); /*00001540*/ END; /*00001550*/ IF IMAG(W(I))<0E0 THEN DO; /*00001560*/ Z=WT;S=R; /*00001570*/ END; /*00001580*/ ELSE DO; /*00001590*/ M=I; /*00001600*/ IF IMAG(W(I))=0E0 THEN DO; /*00001610*/ IF WT^=0E0 THEN /*00001620*/ U=WT; /*00001630*/ ELSE /*00001640*/ U=MACHEPS* /*00001650*/ NORM; /*00001660*/ H(I,EN)=-R/U; /*00001670*/ END; /*00001680*/ ELSE DO; /*00001690*/ /* Solve |wt x||h(i,en) |=|-r| /*00001700*/ | y z||h(i+1,en)| |-s| */ /*00001710*/ X=H(I,I+1); Y=H(I+1,I); /*00001720*/ Q=(REAL(W(I))-P)**2+IMAG(W(I))**2; /*00001730*/ H(I,EN),T=(X*S-Z*R)/Q; /*00001740*/ IF ABS(X)>ABS(Z) THEN U=(-R-WT*T)/X; /*00001750*/ ELSE U=(-S-Y*T)/Z; /*00001760*/ H(I+1,EN)=U; /*00001770*/ END; /*00001780*/ END; /*00001790*/ END; /*00001800*/ END; /*00001810*/ ELSE /*00001820*/ IF Q<0E0 THEN DO; /*00001830*/ /*Complex vector associated with /*00001840*/ lambda=p-i*q */ /*00001850*/ M=NA; /*00001860*/ IF ABS(H(EN,NA))>ABS(H(NA,EN)) THEN /*00001870*/ DO; /*00001880*/ H(NA,NA)=-(H(EN,EN)-P)/H(EN,NA); /*00001890*/ H(NA,EN)=-Q/H(EN,NA); /*00001900*/ END; /*00001910*/ ELSE /*00001920*/ DO; /*00001930*/ TR=-H(NA,EN)/COMPLEX(H(NA,NA)-P,Q); /*00001940*/ H(NA,NA)=REAL(TR); /*00001950*/ H(NA,EN)=IMAG(TR); /*00001960*/ END; /*00001970*/ H(EN,NA)=1E0; H(EN,EN)=0E0; /*00001980*/ DO I=NA-1 TO LBN BY -1; /*00001990*/ WT=H(I,I)-P; RA=H(I,EN); SA=0E0; /*00002000*/ DO J=M TO NA; /*00002010*/ RA=RA+H(I,J)*H(J,NA); /*00002020*/ SA=SA+H(I,J)*H(J,EN); /*00002030*/ END; /*00002040*/ IF IMAG(W(I))<0E0 THEN DO; /*00002050*/ Z=WT;R=RA;S=SA; /*00002060*/ END; /*00002070*/ ELSE DO; /*00002080*/ M=I; /*00002090*/ IF IMAG(W(I))=0E0 /*00002100*/ THEN DO; /*00002110*/ TR=COMPLEX(-RA,-SA)/ /*00002120*/ COMPLEX(WT,Q); /*00002130*/ H(I,NA)=REAL(TR); /*00002140*/ H(I,EN)=IMAG(TR); /*00002150*/ END; /*00002160*/ ELSE DO; /*00002170*/ /*solve complex equations /*00002180*/ |w+qi x ||h(i,na)+h(i,en)i |=|-ra-sa i| /*00002190*/ | y z+qi||h(i+1,na)+h(i+1,en)i| |-r-si | */ /*00002200*/ X=H(I,I+1); Y=H(I+1,I); /*00002210*/ REAL(V)=(REAL(W(I))-P)**2+IMAG(W(I))**2-Q**2; /*00002220*/ IMAG(V)=(REAL(W(I))-P)*2E0*Q; /*00002230*/ IF V=0E0+0E0I THEN V=MACHEPS*NORM*(ABS(WT)+ABS(Q)+ /*00002240*/ ABS(X)+ABS(Y)+ABS(Z)); /*00002250*/ TR=COMPLEX(X*R-Z*RA+Q*SA,X*S-Z*SA-Q*RA)/V; /*00002260*/ H(I,NA)=REAL(TR); H(I,EN)=IMAG(TR); /*00002270*/ IF ABS(X)>ABS(Z)+ABS(Q) /*00002280*/ THEN DO; /*00002290*/ H(I+1,NA)=(-RA-WT*H(I,NA)+Q*H(I,EN))/X; /*00002300*/ H(I+1,EN)=(-SA-WT*H(I,EN)-Q*H(I,NA))/X; /*00002310*/ END; /*00002320*/ ELSE DO; /*00002330*/ TR=COMPLEX(-R-Y*H(I,NA),-S-Y*H(I,EN))/ /*00002340*/ COMPLEX(Z,Q); /*00002350*/ H(I+1,NA)=REAL(TR); H(I+1,EN)=IMAG(TR); /*00002360*/ END; /*00002370*/ END; /*00002380*/ END; /*00002390*/ END; /*00002400*/ END;/*END COMPLEX VECTOR */ /*00002410*/ END;/*END BACKSUBSTITUTION */ /*00002420*/ /* VECTORS OF ISOLATED ROOTS */ /*00002430*/ DO I=LBN TO LOW-1,UPP+1 TO UBN; /*00002440*/ DO J=I+1 TO UBN; /*00002450*/ VECS(I,J)=H(I,J); /*00002460*/ END; /*00002470*/ END; /*00002480*/ /* Multiply by transformation matrix to give vectors of /*00002490*/ original full matrix */ /*00002500*/ DO J=UBN TO LOW BY -1; /*00002510*/ M=MIN(J,UPP); L=J-1; /*00002520*/ IF IMAG(W(J))<0E0 THEN DO; /*00002530*/ DO I=LOW TO UPP; /*00002540*/ P,R=0E0; /*00002550*/ DO K= LOW TO M; /*00002560*/ P=P+VECS(I,K)*H(K,L); /*00002570*/ R=R+VECS(I,K)*H(K,J); /*00002580*/ END; /*00002590*/ VECS(I,L)=P; VECS(I,J)=R; /*00002600*/ END; /*00002610*/ END; /*00002620*/ ELSE DO; /*00002630*/ IF IMAG(W(J))=0E0 /*00002640*/ THEN DO I=LOW TO UPP; /*00002650*/ R=0E0; /*00002660*/ DO K=LOW TO M; /*00002670*/ R=R+VECS(I,K)*H(K,J); /*00002680*/ END; /*00002690*/ VECS(I,J)=R; /*00002700*/ END; /*00002710*/ END; /*00002720*/ END; /*00002730*/ END HQR2; /*00002740*/