%endoflist
%external %routine %spec SELECT INPUT(%integer STREAM)
%external %routine %spec SELECT OUTPUT(%integer STREAM)
%external %routine %spec NEWLINE
%external %routine %spec SPACE
%external %routine %spec SKIP SYMBOL
%external %routine %spec READ STRING(%string %name S)
%external %routine %spec NEWLINES(%integer N)
%external %routine %spec SPACES(%integer N)
%external %integer %fn %spec NEXT SYMBOL
%external %routine %spec PRINT SYMBOL(%integer SYMBOL)
%external %routine %spec READ SYMBOL(%name SYMBOL)
%external %routine %spec READ(%name NUMBER)
%external %routine %spec WRITE(%integer VALUE,PLACES)
%external %routine %spec NEWPAGE
%external %integer %fn %spec ADDR(%name VARIABLE)
%external %long %real %fn %spec ARCSIN(%long %real X)
%external %integer %fn %spec INT(%long %real X)
%external %integer %fn %spec INTPT(%long %real X)
%external %long %real %fn %spec FRACPT(%long %real X)
%external %routine %spec PRINT(%long %real NUMBER,%integer BEFORE,AFTER)
%external %routine %spec PRINTFL(%long %real NUMBER,%integer PLACES)
%external %real %map %spec REAL(%integer VAR ADDR)
%external %integer %map %spec INTEGER(%integer VAR ADDR)
%external %long %real %fn %spec MOD(%long %real X)
%external %long %real %fn %spec ARCCOS(%long %real X)
%external %long %real %fn %spec SQRT(%long %real X)
%external %long %real %fn %spec LOG(%long %real X)
%external %long %real %fn %spec SIN(%long %real X)
%external %long %real %fn %spec COS(%long %real X)
%external %long %real %fn %spec TAN(%long %real X)
%external %long %real %fn %spec EXP(%long %real X)
%external %routine %spec CLOSE STREAM(%integer STREAM)
%external %byte %integer %map %spec BYTE INTEGER(%integer VAR ADDR)
%external %integer %fn %spec EVENTINF
%external %long %real %fn %spec RADIUS(%long %real X,Y)
%external %long %real %fn %spec ARCTAN(%long %real X,Y)
%external %byte %integer %map %spec LENGTH(%string %name  S)
%external %routine %spec PRINT STRING(%string(255) STR)
%external %integer %fn %spec NL
%external %long %real %map %spec LONG REAL(%integer VAR ADDR)
%external %routine %spec PRINT CH(%integer CHARACTER)
%external %routine %spec READ CH(%name CHARACTER)
%external %string %map %spec STRING(%integer VAR ADDR)
%external %routine %spec READ ITEM(%string %name ITEM)
%external %string (1) %fn %spec NEXT ITEM
%external %byte %integer %map %spec CHARNO(%string %name STR,%integer CHAR REQD)
%external %string (1) %fn %spec TO STRING(%integer ch)
%external %string (255) %fn %spec SUB STRING(%string(255) STR, %integer FIRST, LAST)
!%external %record (*) %map %spec RECORD(%integer REC ADDR)
!%external %array %map %spec ARRAY(%integer A1ADDR,%array %name FORMAT)
%external %integer %fn %spec SIZEOF(%name X)
%external %integer %fn %spec IMOD(%integer VALUE)
%external %long %real %fn %spec PI
%external %integer %fn %spec EVENTLINE
%external %long %integer %map %spec LONGINTEGER(%integer ADR)
%external %long %long %real %map %spec LONGLONGREAL(%integer ADR)
%external %long %integer %fn %spec LENGTHENI(%integer VAL)
%external %long %long %real %fn %spec LENGTHENR(%long %real VAL)
%external %integer %fn %spec SHORTENI(%long %integer VAL)
%external %long %real %fn %spec SHORTENR(%long %long %real VAL)
%external %integer %fn %spec NEXTCH
%external %half %integer %map %spec HALFINTEGER(%integer ADDR)
%external %routine %spec PPROFILE
%external %long %real %fn %spec FLOAT(%integer VALUE)
%external %long %integer %fn %spec LINT(%long %long %real X)
%external %long %integer %fn %spec LINTPT(%long %long %real X)
%external %short %integer %map %spec SHORTINTEGER(%integer N)
%external %integer %fn %spec TRUNC(%long %real X)
!%endoffile
%list

%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
%CONSTLONGREAL 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)<LAMBDA %THEN LAMBDA=W(J)
         U=C(1)-LAMBDA; V=B(1)
         %IF V=0 %THEN V=EPS
         %CYCLE I=1,1,N-1
            BI=B(I)
            %IF BI=0 %THEN BI=EPS
            BI1=B(I+1)
            %IF BI1=0 %THEN BI1=EPS
            %IF MOD(U)>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