C C .................................................................. C C SUBROUTINE ATEIG C C PURPOSE C COMPUTE THE EIGENVALUES OF A REAL ALMOST TRIANGULAR MATRIX C C USAGE C CALL ATEIG(M,A,RR,RI,IANA,IA) C C DESCRIPTION OF THE PARAMETERS C M ORDER OF THE MATRIX C A THE INPUT MATRIX, M BY M C RR VECTOR CONTAINING THE REAL PARTS OF THE EIGENVALUES C ON RETURN C RI VECTOR CONTAINING THE IMAGINARY PARTS OF THE EIGEN- C VALUES ON RETURN C IANA VECTOR WHOSE DIMENSION MUST BE GREATER THAN OR EQUAL C TO M, CONTAINING ON RETURN INDICATIONS ABOUT THE WAY C THE EIGENVALUES APPEARED (SEE MATH. DESCRIPTION) C IA SIZE OF THE FIRST DIMENSION ASSIGNED TO THE ARRAY A C IN THE CALLING PROGRAM WHEN THE MATRIX IS IN DOUBLE C SUBSCRIPTED DATA STORAGE MODE. C IA=M WHEN THE MATRIX IS IN SSP VECTOR STORAGE MODE. C C REMARKS C THE ORIGINAL MATRIX IS DESTROYED C THE DIMENSION OF RR AND RI MUST BE GREATER OR EQUAL TO M C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C QR DOUBLE ITERATION C C REFERENCES C J.G.F. FRANCIS - THE QR TRANSFORMATION---THE COMPUTER C JOURNAL, VOL. 4, NO. 3, OCTOBER 1961, VOL. 4, NO. 4, JANUARY C 1962. J. H. WILKINSON - THE ALGEBRAIC EIGENVALUE PROBLEM - C CLARENDON PRESS, OXFORD, 1965. C C .................................................................. C SUBROUTINE ATEIG(M,A,RR,RI,IANA,IA) DIMENSION A(1),RR(1),RI(1),PRR(2),PRI(2),IANA(1) INTEGER P,P1,Q C E7=1.0E-8 E6=1.0E-6 E10=1.0E-10 DELTA=0.5 MAXIT=30 C C INITIALIZATION C N=M 20 N1=N-1 IN=N1*IA NN=IN+N IF(N1) 30,1300,30 30 NP=N+1 C C ITERATION COUNTER C IT=0 C C ROOTS OF THE 2ND ORDER MAIN SUBMATRIX AT THE PREVIOUS C ITERATION C DO 40 I=1,2 PRR(I)=0.0 40 PRI(I)=0.0 C C LAST TWO SUBDIAGONAL ELEMENTS AT THE PREVIOUS ITERATION C PAN=0.0 PAN1=0.0 C C ORIGIN SHIFT C R=0.0 S=0.0 C C ROOTS OF THE LOWER MAIN 2 BY 2 SUBMATRIX C N2=N1-1 IN1=IN-IA NN1=IN1+N N1N=IN+N1 N1N1=IN1+N1 60 T=A(N1N1)-A(NN) U=T*T V=4.0*A(N1N)*A(NN1) IF(ABS(V)-U*E7) 100,100,65 65 T=U+V IF(ABS(T)-AMAX1(U,ABS(V))*E6) 67,67,68 67 T=0.0 68 U=(A(N1N1)+A(NN))/2.0 V=SQRT(ABS(T))/2.0 IF(T)140,70,70 70 IF(U) 80,75,75 75 RR(N1)=U+V RR(N)=U-V GO TO 130 80 RR(N1)=U-V RR(N)=U+V GO TO 130 100 IF(T)120,110,110 110 RR(N1)=A(N1N1) RR(N)=A(NN) GO TO 130 120 RR(N1)=A(NN) RR(N)=A(N1N1) 130 RI(N)=0.0 RI(N1)=0.0 GO TO 160 140 RR(N1)=U RR(N)=U RI(N1)=V RI(N)=-V 160 IF(N2)1280,1280,180 C C TESTS OF CONVERGENCE C 180 N1N2=N1N1-IA RMOD=RR(N1)*RR(N1)+RI(N1)*RI(N1) EPS=E10*SQRT(RMOD) IF(ABS(A(N1N2))-EPS)1280,1280,240 240 IF(ABS(A(NN1))-E10*ABS(A(NN))) 1300,1300,250 250 IF(ABS(PAN1-A(N1N2))-ABS(A(N1N2))*E6) 1240,1240,260 260 IF(ABS(PAN-A(NN1))-ABS(A(NN1))*E6)1240,1240,300 300 IF(IT-MAXIT) 320,1240,1240 C C COMPUTE THE SHIFT C 320 J=1 DO 360 I=1,2 K=NP-I IF(ABS(RR(K)-PRR(I))+ABS(RI(K)-PRI(I))-DELTA*(ABS(RR(K)) 1 +ABS(RI(K)))) 340,360,360 340 J=J+I 360 CONTINUE GO TO (440,460,460,480),J 440 R=0.0 S=0.0 GO TO 500 460 J=N+2-J R=RR(J)*RR(J) S=RR(J)+RR(J) GO TO 500 480 R=RR(N)*RR(N1)-RI(N)*RI(N1) S=RR(N)+RR(N1) C C SAVE THE LAST TWO SUBDIAGONAL TERMS AND THE ROOTS OF THE C SUBMATRIX BEFORE ITERATION C 500 PAN=A(NN1) PAN1=A(N1N2) DO 520 I=1,2 K=NP-I PRR(I)=RR(K) 520 PRI(I)=RI(K) C C SEARCH FOR A PARTITION OF THE MATRIX, DEFINED BY P AND Q C P=N2 IF (N-3)600,600,525 525 IPI=N1N2 DO 580 J=2,N2 IPI=IPI-IA-1 IF(ABS(A(IPI))-EPS) 600,600,530 530 IPIP=IPI+IA IPIP2=IPIP+IA D=A(IPIP)*(A(IPIP)-S)+A(IPIP2)*A(IPIP+1)+R IF(D)540,560,540 540 IF(ABS(A(IPI)*A(IPIP+1))*(ABS(A(IPIP)+A(IPIP2+1)-S)+ABS(A(IPIP2+2) 1 )) -ABS(D)*EPS) 620,620,560 560 P=N1-J 580 CONTINUE 600 Q=P GO TO 680 620 P1=P-1 Q=P1 IF (P1-1) 680,680,650 650 DO 660 I=2, P1 IPI=IPI-IA-1 IF(ABS(A(IPI))-EPS)680,680,660 660 Q=Q-1 C C QR DOUBLE ITERATION C 680 II=(P-1)*IA+P DO 1220 I=P,N1 II1=II-IA IIP=II+IA IF(I-P)720,700,720 700 IPI=II+1 IPIP=IIP+1 C C INITIALIZATION OF THE TRANSFORMATION C G1=A(II)*(A(II)-S)+A(IIP)*A(IPI)+R G2=A(IPI)*(A(IPIP)+A(II)-S) G3=A(IPI)*A(IPIP+1) A(IPI+1)=0.0 GO TO 780 720 G1=A(II1) G2=A(II1+1) IF(I-N2)740,740,760 740 G3=A(II1+2) GO TO 780 760 G3=0.0 780 CAP=SQRT(G1*G1+G2*G2+G3*G3) IF(CAP)800,860,800 800 IF(G1)820,840,840 820 CAP=-CAP 840 T=G1+CAP PSI1=G2/T PSI2=G3/T ALPHA=2.0/(1.0+PSI1*PSI1+PSI2*PSI2) GO TO 880 860 ALPHA=2.0 PSI1=0.0 PSI2=0.0 880 IF(I-Q)900,960,900 900 IF(I-P)920,940,920 920 A(II1)=-CAP GO TO 960 940 A(II1)=-A(II1) C C ROW OPERATION C 960 IJ=II DO 1040 J=I,N T=PSI1*A(IJ+1) IF(I-N1)980,1000,1000 980 IP2J=IJ+2 T=T+PSI2*A(IP2J) 1000 ETA=ALPHA*(T+A(IJ)) A(IJ)=A(IJ)-ETA A(IJ+1)=A(IJ+1)-PSI1*ETA IF(I-N1)1020,1040,1040 1020 A(IP2J)=A(IP2J)-PSI2*ETA 1040 IJ=IJ+IA C C COLUMN OPERATION C IF(I-N1)1080,1060,1060 1060 K=N GO TO 1100 1080 K=I+2 1100 IP=IIP-I DO 1180 J=Q,K JIP=IP+J JI=JIP-IA T=PSI1*A(JIP) IF(I-N1)1120,1140,1140 1120 JIP2=JIP+IA T=T+PSI2*A(JIP2) 1140 ETA=ALPHA*(T+A(JI)) A(JI)=A(JI)-ETA A(JIP)=A(JIP)-ETA*PSI1 IF(I-N1)1160,1180,1180 1160 A(JIP2)=A(JIP2)-ETA*PSI2 1180 CONTINUE IF(I-N2)1200,1220,1220 1200 JI=II+3 JIP=JI+IA JIP2=JIP+IA ETA=ALPHA*PSI2*A(JIP2) A(JI)=-ETA A(JIP)=-ETA*PSI1 A(JIP2)=A(JIP2)-ETA*PSI2 1220 II=IIP+1 IT=IT+1 GO TO 60 C C END OF ITERATION C 1240 IF(ABS(A(NN1))-ABS(A(N1N2))) 1300,1280,1280 C C TWO EIGENVALUES HAVE BEEN FOUND C 1280 IANA(N)=0 IANA(N1)=2 N=N2 IF(N2)1400,1400,20 C C ONE EIGENVALUE HAS BEEN FOUND C 1300 RR(N)=A(NN) RI(N)=0.0 IANA(N)=1 IF(N1)1400,1400,1320 1320 N=N1 GO TO 20 1400 RETURN END