%control x'11111111'

%BEGIN
%CONSTINTEGER NMAX=20
%ROUTINESPEC TEST MATRIX(%REALARRAYNAME A, %INTEGER N)
%ROUTINESPEC HOUSEHOLDER(%REALARRAYNAME A,W, %INTEGER N,K)
%INTEGER I
%REALARRAY A(1:NMAX,1:NMAX)
%REALARRAY 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
      pprofile
      printstring("
Test is OK provided all answers are 1.00000 except last 2 which should be
1.63669..@ -2 and -2.49383..@-2
")
      %STOP
!% %REAL %FN SQRT(% %REAL ARG)
!!***********************************************************************
!!*    THIS ROUTINE SCALES WITH A LOOP AND HENCE IS VERY SLOW           *
!!*    MUCHFASTER M-C DEPENDENT SCALING IS POSSIBLE                     *
!!***********************************************************************
!%INTEGER I,J,h,n
!% %REAL Y,OLD,NEW,f,ta
!      %IF ARG<0.0 %THEN ARG=-ARG
!      %IF ARG=0 %THEN %RESULT=0.0
!      i=0
!      y=arg
!      y=y/16 %AND i=1+1 %WHILE y>1.0
!      y=y*16.0 %AND i=i-1 %WHILE y<1/16
!      %IF Y>0.34375 %THEN OLD=Y/2+0.4368 %ELSE OLD=Y/2+0.381
!      %CYCLE J=0,1,14
!         NEW=(OLD+Y/OLD)/2
!         %IF MOD(NEW-OLD)<0.000000000000005 %THEN %EXIT
!         OLD=NEW
!      %REPEAT
!      %IF I>0 %THEN NEW=NEW*4.0**i
!      %IF i<0 %THEN new=new/(4.0**imod(i))
!      %RESULT=NEW
!%END
%ROUTINE TEST MATRIX(%REALARRAYNAME A, %INTEGER N)
%INTEGER I,J
%REAL 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(%REALARRAYNAME A,W, %INTEGER N,K)
%ROUTINESPEC HOUSEHOLDER TRIDIAGONALISATION(%REALARRAYNAME A,B,C,
    %INTEGER N)
%ROUTINESPEC TRIDIBISECTION(%REALARRAYNAME C,B,W, %INTEGER N,
    %REALNAME NORM)
%ROUTINESPEC TRIDIINVERSE ITERATION(%REALARRAYNAME C,B,W,Z, %INTEGER N,
    %REAL NORM)
%ROUTINESPEC BACKTRANSFORMATION(%REALARRAYNAME A,B,Z, %INTEGER N)
%REALARRAY Z(1:N,1:N)
%REALARRAY B,C(1:N)
%INTEGER I,J
%REAL NORM
%CONSTREAL 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(%REALARRAYNAME A,B,C, %INTEGER N)
%INTEGER I,J,K
%REAL AI,SIGMA,H,BJ,BIGK,BI
%REALARRAY 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(%REALARRAYNAME C,B,W,Z, %INTEGER N,
    %REAL NORM)
%INTEGER I,J
%REAL BI,BI1,Z1,LAMBDA,U,S,V,H,EPS,ETA
%REALARRAY M,P,Q,R,INT(1:N)
%REALARRAY 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)**2
         %REPEAT
         H=1/SQRT(H)
         %CYCLE I=1,1,N
            Z(J,I)=X(I)*H
         %REPEAT
      %REPEAT
%END
%ROUTINE TRIDIBISECTION(%REALARRAYNAME C,B,W, %INTEGER N,
    %REALNAME NORM)
%REAL L,G,H,LAMBDA,P1,Q1,Y
%INTEGER I,J,K,A1,A2
%REALARRAY 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)**2
      %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(%REALARRAYNAME A,B,Z, %INTEGER N)
%INTEGER I,J,K
%REAL 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