%BEGIN %CONSTINTEGER NMAX=20 %ROUTINESPEC TEST MATRIX(%LONGREALARRAYNAME A, %INTEGER N) %ROUTINESPEC HOUSEHOLDER(%LONGREALARRAYNAME A,W, %INTEGER N,K) %INTEGER I %LONGREALARRAY A(1:NMAX,1:NMAX) %LONGREALARRAY W(1:NMAX) TEST MATRIX(A,NMAX) PRINTSTRING(" TEST MATRIX EIGENVALUES ") HOUSEHOLDER(A,W,NMAX,1) %CYCLE I=1,1,NMAX PRINT FL(W(I),8) NEWLINE %REPEAT %STOP %ROUTINE TEST MATRIX(%LONGREALARRAYNAME A, %INTEGER N) %INTEGER I,J %LONGREAL T,C,D,F C=N*(N+1)*(2*N-5)/6 D=1/C A(N,N)=-D %CYCLE I=1,1,N-1 F=I A(I,N)=D*F A(N,I)=A(I,N) A(I,I)=D*(C-F**2) %CYCLE J=1,1,I-1 T=J A(I,J)=-D*F*T A(J,I)=A(I,J) %REPEAT %REPEAT ! %CYCLE I=1,1,N ! %CYCLE J=1,1,N ! PRINTFL(A(I,J),4) ! %REPEAT ! NEWLINE ! %REPEAT %END %ROUTINE HOUSEHOLDER(%LONGREALARRAYNAME A,W, %INTEGER N,K) %ROUTINESPEC HOUSEHOLDER TRIDIAGONALISATION(%LONGREALARRAYNAME A,B,C, %INTEGER N) %ROUTINESPEC TRIDIBISECTION(%LONGREALARRAYNAME C,B,W, %INTEGER N, %LONGREALNAME NORM) %ROUTINESPEC TRIDIINVERSE ITERATION(%LONGREALARRAYNAME C,B,W,Z, %INTEGER N, %LONGREAL NORM) %ROUTINESPEC BACKTRANSFORMATION(%LONGREALARRAYNAME A,B,Z, %INTEGER N) %LONGREALARRAY Z(1:N,1:N) %LONGREALARRAY B,C(1:N) %INTEGER I,J %LONGREAL NORM %LONGREAL TEN TEN=10 HOUSEHOLDER TRIDIAGONALISATION(A,B,C,N) TRIDIBISECTION(C,B,W,N,NORM) %IF K=2 %THENRETURN ! %CYCLE I=1,1,N ! PRINTFL(B(I),6) ! PRINTFL(C(I),6) ! PRINTFL(W(I),6) ! NEWLINE ! %REPEAT TRIDIINVERSE ITERATION(C,B,W,Z,N,NORM) BACKTRANSFORMATION(A,B,Z,N) %CYCLE I=1,1,N %CYCLE J=1,1,N A(I,J)=Z(I,J) %REPEAT %REPEAT %RETURN %ROUTINE HOUSEHOLDER TRIDIAGONALISATION(%LONGREALARRAYNAME A,B,C, %INTEGER N) %INTEGER I,J,K %LONGREAL AI,SIGMA,H,BJ,BIGK,BI %LONGREALARRAY Q(1:N-1) %CYCLE I=N,-1,3 SIGMA=0 %CYCLE K=1,1,I-1 SIGMA=SIGMA+A(I,K)**2 %REPEAT AI=A(I,I-1) %IF AI>=0 %THEN BI=-SQRT(SIGMA) %ELSE BI=SQRT(SIGMA) B(I-1)=BI %IF BI=0 %THENCONTINUE H=SIGMA-AI*BI A(I,I-1)=AI-BI %CYCLE J=I-1,-1,1 BJ=0 %CYCLE K=I-1,-1,J BJ=BJ+A(K,J)*A(I,K) %REPEAT %CYCLE K=J-1,-1,1 BJ=BJ+A(J,K)*A(I,K) %REPEAT Q(J)=BJ/H %REPEAT BIGK=0 %CYCLE J=I-1,-1,1 BIGK=BIGK+A(I,J)*Q(J) %REPEAT BIGK=BIGK/(2*H) %CYCLE J=I-1,-1,1 Q(J)=Q(J)-BIGK*A(I,J) %REPEAT %CYCLE J=I-1,-1,1 %CYCLE K=J,-1,1 A(J,K)=A(J,K)-A(I,J)*Q(K)-A(I,K)*Q(J) %REPEAT %REPEAT %REPEAT %CYCLE I=N,-1,1 C(I)=A(I,I) %REPEAT B(1)=A(2,1) B(N)=0 %END %ROUTINE TRIDIINVERSE ITERATION(%LONGREALARRAYNAME C,B,W,Z, %INTEGER N, %LONGREAL NORM) %INTEGER I,J %LONGREAL BI,BI1,Z1,LAMBDA,U,S,V,H,EPS,ETA %LONGREALARRAY M,P,Q,R,INT(1:N) %LONGREALARRAY X(1:N+2) LAMBDA=NORM EPS=NORM/(TEN**11) %CYCLE J=1,1,N LAMBDA=LAMBDA-EPS %IF W(J)MOD(BI) %THENSTART M(I+1)=BI/U P(I)=U; Q(I)=V; R(I)=0 U=C(I+1)-LAMBDA-M(I+1)*V V=BI1; INT(I+1)=-1 %FINISHELSESTART M(I+1)=U/BI %IF M(I+1)=0 %AND BI<=EPS %THEN M(I+1)=1 P(I)=BI; Q(I)=C(I+1)-LAMBDA R(I)=BI1 U=V-M(I+1)*Q(I) V=-M(I+1)*R(I) INT(I+1)=1 %FINISH %REPEAT P(N)=U; Q(N)=0; R(N)=0 X(N+1)=0; X(N+2)=0; H=0; ETA=1/N %CYCLE I=N,-1,1 U=ETA-Q(I)*X(I+1)-R(I)*X(I+2) %IF P(I)\=0 %THEN X(I)=U/P(I) %ELSE X(I)=U/EPS H=H+MOD(X(I)) %REPEAT H=1/H %CYCLE I=1,1,N X(I)=X(I)*H %REPEAT %CYCLE I=2,1,N %IF INT(I)<=0 %THEN X(I)=X(I)-M(I)*X(I-1) %ELSESTART U=X(I-1) X(I-1)=X(I) X(I)=U-M(I)*X(I-1) %FINISH %REPEAT H=0 %CYCLE I=N,-1,1 U=X(I)-Q(I)*X(I+1)-R(I)*X(I+2) %IF P(I)=0 %THEN X(I)=U/EPS %ELSE X(I)=U/P(I) H=H+X(I)*X(I) %REPEAT H=1/SQRT(H) %CYCLE I=1,1,N Z(J,I)=X(I)*H %REPEAT %REPEAT %END %ROUTINE TRIDIBISECTION(%LONGREALARRAYNAME C,B,W, %INTEGER N, %LONGREALNAME NORM) %LONGREAL L,G,H,LAMBDA,P1,Q1,Y %INTEGER I,J,K,A1,A2 %LONGREALARRAY P(1:N) P(1)=0 NORM=MOD(C(1))+MOD(B(1)) %CYCLE I=2,1,N L=MOD(B(I-1))+MOD(C(I))+MOD(B(I)) %IF L>NORM %THEN NORM=1 %REPEAT %CYCLE I=1,1,N-1 %IF B(I)=0 %THEN P(I+1)=NORM*NORM/TEN**23 %ELSE P(I+1)=B(I)*B(I) %REPEAT %CYCLE K=1,1,N G=NORM; H=-NORM %CYCLE J=1,1,39 LAMBDA=(G+H)/2 P1=0; Q1=1; A1=0 %CYCLE I=1,1,N Y=(C(I)-LAMBDA)*Q1-P(I)*P1 P1=Q1; Q1=Y %IF (P1>=0 %AND Q1>=0) %OR (P1<0 %AND Q1<0) %THEN A1=A1+1 %REPEAT %IF Q1=0 %AND P1>0 %THEN A1=A1-1 %IF A1>=K %THEN H=LAMBDA %ELSE G=LAMBDA %REPEAT W(K)=(G+H)/2 %REPEAT %RETURN %END %ROUTINE BACKTRANSFORMATION(%LONGREALARRAYNAME A,B,Z, %INTEGER N) %INTEGER I,J,K %LONGREAL S %CYCLE J=1,1,N %CYCLE K=3,1,N %IF B(K-1)=0 %THENCONTINUE S=0 %CYCLE I=1,1,K-1 S=S+A(K,I)*Z(J,I) %REPEAT S=S/(B(K-1)*A(K,K-1)) %CYCLE I=1,1,K-1 Z(J,I)=Z(J,I)+S*A(K,I) %REPEAT %REPEAT %REPEAT %END %END %ENDOFPROGRAM