SUBROUTINE EIGENZ(T,R,VA,A,N,NX,NV) C EIG0010 IMPLICIT REAL *8(A-H,O-Z) EIG0020 DIMENSION T(NX,NX),VA(NX),A(6),R(9) C EIG0050 K=0 EIG0060 DO 10 J=1,N EIG0070 DO 10 I=1,J EIG0080 K=K+1 EIG0090 A(K)=T(I,J) EIG0100 10 CONTINUE EIG0110 NN=N*(N+1)/2 EIG0120 20 IF (NV-1) 30,60,30 EIG0130 30 IQ=-N EIG0140 DO 50 J=1,N EIG0150 IQ=IQ+N EIG0160 DO 50 I=1,N EIG0170 IJ=IQ+I EIG0180 R(IJ)=0.0D0 EIG0190 IF (I-J) 50,40,50 EIG0200 40 R(IJ)=1.0D0 EIG0210 50 CONTINUE EIG0220 C EIG0230 C COMPUTE INITIAL AND FINAL NORMS ANORM AND ANORMX EIG0240 C EIG0250 60 ANORM=0.0D0 EIG0260 DO 80 I=1,N EIG0270 DO 80 J=I,N EIG0280 IF (I-J) 70,80,70 EIG0290 70 IA=I+(J*J-J)/2 EIG0300 ANORM=ANORM+A(IA)*A(IA) EIG0310 80 CONTINUE EIG0320 IF (ANORM) 370,370,90 EIG0330 90 ANORM=DSQRT(2.0D2*ANORM) EIG0340 ANRMX=ANORM*1.0D-06/DFLOAT(N) EIG0350 C EIG0360 C INITIALIZE INDICATORS AND COMPUTE THRESHOULD, THR EIG0370 C EIG0380 IND=0 EIG0390 THR=ANORM EIG0400 100 THR=THR/DFLOAT(N) EIG0410 110 L=1 EIG0420 120 M=L+1 EIG0430 C EIG0440 C COMPUTE SIN AND COS EIG0450 C EIG0460 130 MQ=(M*M-M)/2 EIG0470 LQ=(L*L-L)/2 EIG0480 LM=L+MQ EIG0490 140 IF (DABS(A(LM))-THR) 300,150,150 EIG0500 150 IND=1 EIG0510 LL=L+LQ EIG0520 MM=M+MQ EIG0530 X=0.5D0*(A(LL)-A(MM)) EIG0540 160 Y=-A(LM)/DSQRT(A(LM)*A(LM)+X*X) EIG0550 IF (X) 170,180,180 EIG0560 170 Y=-Y EIG0570 180 SINX=Y/DSQRT(2.0D0*(1.0D0+(DSQRT(1.0D0-Y*Y)))) EIG0580 SINX2=SINX*SINX EIG0590 190 COSX=1.0D0-SINX2 EIG0600 COSX=DSQRT(COSX) EIG0610 COSX2=COSX*COSX EIG0620 SINCS=SINX*COSX EIG0630 C EIG0640 C ROTATE L AND M COLUMNS EIG0650 C EIG0660 ILQ=N*(L-1) EIG0670 IMQ=N*(M-1) EIG0680 DO 290 I=1,N EIG0690 IQ=(I*I-I)/2 EIG0700 IF (I-L) 200,270,200 EIG0710 200 IF (I-M) 210,270,220 EIG0720 210 IM=I+MQ EIG0730 GO TO 230 EIG0740 220 IM=M+IQ EIG0750 230 IF (I-L) 240,250,250 EIG0760 240 IL=I+LQ EIG0770 GO TO 260 EIG0780 250 IL=L+IQ EIG0790 260 X=A(IL)*COSX-A(IM)*SINX EIG0800 A(IM)=A(IL)*SINX+A(IM)*COSX EIG0810 A(IL)=X EIG0820 270 IF (NV-1) 280,290,280 EIG0830 280 ILR=ILQ+I EIG0840 IMR=IMQ+I EIG0850 X=R(ILR)*COSX-R(IMR)*SINX EIG0860 R(IMR)=R(ILR)*SINX+R(IMR)*COSX EIG0870 R(ILR)=X EIG0880 290 CONTINUE EIG0890 X=2.0D0*A(LM)*SINCS EIG0900 Y=A(LL)*COSX2+A(MM)*SINX2-X EIG0910 X=A(LL)*SINX2+A(MM)*COSX2+X EIG0920 A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) EIG0930 A(LL)=Y EIG0940 A(MM)=X EIG0950 C EIG0960 C TESTS FOR COMPLITION EIG0970 C EIG0980 C TEST FOR M LAST COLUMN EIG0990 C EIG1000 300 IF (M-N) 310,320,310 EIG1010 310 M=M+1 EIG1020 GO TO 130 EIG1030 C EIG1040 C TEST FOR L SECOND FROM LAST COLUMN EIG1050 C EIG1060 320 IF (L-(N-1)) 330,340,330 EIG1070 330 L=L+1 EIG1080 GO TO 120 EIG1090 340 IF (IND-1) 360,350,360 EIG1100 350 IND=0 EIG1110 GO TO 110 EIG1120 C EIG1130 C COMPARE THRESHOULD WITH FINAL NORM EIG1140 C EIG1150 360 IF (THR-ANRMX) 370,370,100 EIG1160 C EIG1170 C SORT EIGENVALUES AND EIGENVECTORS EIG1180 C EIG1190 370 IQ=-N EIG1200 DO 410 I=1,N EIG1210 IQ=IQ+N EIG1220 LL=I+(I*I-I)/2 EIG1230 JQ=N*(I-2) EIG1240 DO 410 J=I,N EIG1250 JQ=JQ+N EIG1260 MM=J+(J*J-J)/2 EIG1270 IF (A(LL)-A(MM)) 380,410,410 EIG1280 380 X=A(LL) EIG1290 A(LL)=A(MM) EIG1300 A(MM)=X EIG1310 IF (NV-1) 390,410,390 EIG1320 390 DO 400 K=1,N EIG1330 ILR=IQ+K EIG1340 IMR=JQ+K EIG1350 X=R(ILR) EIG1360 R(ILR)=R(IMR) EIG1370 400 R(IMR)=X EIG1380 410 CONTINUE EIG1390 K=0 EIG1400 DO 420 I=1,N EIG1410 K=K+I EIG1420 420 VA(I)=A(K) EIG1430 RETURN EIG1440 END EIG1450