%BEGIN %CONSTINTEGER MSIZE=4 %ROUTINESPEC TEST MATRIX(%LONGREALARRAYNAME A, %INTEGER N) %ROUTINESPEC HOUSEHOLDER(%LONGREALARRAYNAME A, W, %C %INTEGER N, K) %INTEGER I %LONGREALARRAY A(1:MSIZE,1:MSIZE) %LONGREALARRAY W(1:MSIZE) TEST MATRIX(A,MSIZE) PRINTSTRING(" TEST MATRIX EIGENVALUES ") HOUSEHOLDER(A,W,MSIZE,1) %CYCLE I=1,1,MSIZE PRINT FL(W(I),8) NEWLINE %REPEAT %STOP %LONGREALFN SQRT(%LONGREAL ARG) %INTEGER I,J %LONGREAL Y,OLD %LONGREALARRAY A(0:15) ! %IF ARG<0.0 %THEN MLIBERR(50) %AND ARG=-ARG %IF ARG<0.0 %THEN PRINTSTRING("SQRT ARG NEG ") %AND %RESULT=-1 %IF ARG=0 %THEN %RESULT=0.0 %CYCLE I=0,1,15 Y=ARG/2**I %IF Y<1 %THEN ->OUT %REPEAT OUT: %IF Y>1 %THEN I=17 %AND Y=Y/4 %ELSE I=I+1 %AND Y=Y/2 %IF Y>0.34375 %THEN A(0)=Y/2 +0.4368 %ELSE A(0)=Y/2+0.381 %CYCLE J=0,1,14 OLD=A(J) A(J+1)=(OLD+Y/OLD)/2 %IF MOD(A(J+1)-OLD)<0.000000000000005 %THEN %EXIT %REPEAT Y=A(J+1)*2**(I//2) %IF I&1#0 %THEN Y=Y*1.4142136 %RESULT=Y %END %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) %IF I=1 %THEN ->LAB1 %CYCLE J=1,1,I-1 T=J A(I,J)=-D*F*T A(J,I)=A(I,J) %REPEAT LAB1: %REPEAT %END %ROUTINE HOUSEHOLDER(%LONGREALARRAYNAME A, W, %INTEGER N, K) %ROUTINESPEC HOUSEHOLDER TRIDIAGONALISATION( %C %LONGREALARRAYNAME A, B, C, %INTEGER N) %ROUTINESPEC TRIDIBISECTION( %C %LONGREALARRAYNAME C, B, W, %INTEGER N, %C %LONGREALNAME NORM) %ROUTINESPEC TRIDIINVERSE ITERATION( %C %LONGREALARRAYNAME C, B, W, Z, %INTEGER N, %C %LONGREAL NORM) %ROUTINESPEC BACKTRANSFORMATION( %C %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 %THEN %RETURN 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( %C %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 ->LAB1 BI=SQRT(SIGMA) ->LAB2 LAB1: BI=-SQRT(SIGMA) LAB2: B(I-1)=BI %IF BI=0 %THEN ->LAB3 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 %IF J=1 %THEN ->LAB10 %CYCLE K=J-1,-1,1 BJ=BJ+A(J,K)*A(I,K) %REPEAT LAB10: 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 LAB3: %REPEAT %CYCLE I=N,-1,1 C(I)=A(I,I) %REPEAT B(1)=A(2,1) B(N)=0 %END %ROUTINE TRIDIINVERSE ITERATION( %C %LONGREALARRAYNAME C, B, W, Z, %INTEGER N, %C %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) %THEN ->LAB 1 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 ->LAB2 LAB1: 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 LAB2: %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 ->LAB3 X(I)=U/EPS ->LAB4 LAB3: X(I)=U/P(I) LAB4: 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 ->LAB5 U=X(I-1) X(I-1)=X(I) X(I)=U-M(I)*X(I-1) ->LAB6 LAB5: X(I)=X(I)-M(I)*X(I-1) LAB6: %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 ->LAB7 X(I)=U/EPS ->LAB8 LAB7: X(I)=U/P(I) LAB8: 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, %C %INTEGER N, %LONGREALNAME NORM) %ROUTINESPEC STURMS SEQUENCE %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 ->LAB 1 P(I+1)=B(I)*B(I) ->LAB2 LAB1: P(I+1)=NORM*NORM/TEN**23 LAB2: %REPEAT %CYCLE K=1,1,N G=NORM; H=-NORM %CYCLE J=1,1,39 LAMBDA=(G+H)/2 STURMS SEQUENCE %IF A1>=K %THEN ->LAB3 G=LAMBDA ->LAB4 LAB3: H=LAMBDA LAB4: %REPEAT W(K)=(G+H)/2 %REPEAT %RETURN %ROUTINE STURMS SEQUENCE 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) %C %THEN A1=A1+1 %REPEAT %IF Q1=0 %AND P1>0 %THEN A1=A1-1 %END %END %ROUTINE BACKTRANSFORMATION( %C %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 %THEN ->LAB1 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 LAB1: %REPEAT %REPEAT %END %END %ENDOFPROGRAM