code 34186;
    comment MCA 2416;
    integer procedure REAQRI(A, N, EM, VAL, VEC); value N;
    integer N; array A, EM, VAL, VEC;
    begin integer M1, I, I1, M, J, Q, MAX, COUNT;
        real W, SHIFT, KAPPA, NU, MU, R, TOL, S, MACHTOL,
        ELMAX, T, DELTA, DET;
        array TF[1:N];

        MACHTOL:= EM[0] * EM[1]; TOL:= EM[1] * EM[2]; MAX:= EM[4];
        COUNT:= 0; ELMAX:= 0; M:= N;
        for I:= 1 step 1 until N do 
        begin VEC[I,I]:= 1;
            for J:= I + 1 step 1 until N do 
            VEC[I,J]:= VEC[J,I]:= 0
        end;
     IN: M1:= M - 1;
        for I:= M, I - 1 while (if I >= 1 then 
        ABS(A[I + 1,I]) > TOL else false ) do Q:= I;
        if Q > 1 then 
        begin if ABS(A[Q,Q - 1]) > ELMAX then 
            ELMAX:= ABS(A[Q, Q - 1])
        end;
        if Q = M then 
        begin VAL[M]:= A[M,M]; M:= M1 end 
        else 
        begin DELTA:= A[M,M] - A[M1,M1]; DET:= A[M,M1] * A[M1,M];
            if ABS(DELTA) < MACHTOL then S:= SQRT(DET) else 
            begin W:= 2 / DELTA; S:= W * W * DET + 1;
                S:= if S <= 0 then -DELTA * .5 else 
                W * DET / (SQRT(S) + 1)
            end;
            if Q = M1 then 
            begin A[M,M]:= VAL[M]:= A[M,M] + S;
                A[Q,Q]:= VAL[Q]:= A[Q,Q] - S;
                T:= if ABS(S) < MACHTOL then 
                (S + DELTA) / A[M,Q] else A[Q,M] / S;
                R:= SQRT(T * T + 1); NU:= 1 / R;
                MU:= -T * NU; A[Q,M]:= A[Q,M] - A[M,Q];
                ROTROW(Q + 2, N, Q, M, A, MU, NU);
                ROTCOL(1, Q - 1, Q, M, A, MU, NU);
                ROTCOL(1, N, Q, M, VEC, MU, NU); M:= M - 2
            end 

            else 
            begin COUNT:= COUNT + 1;
                if COUNT > MAX then goto END;
                SHIFT:= A[M,M] + S; if ABS(DELTA) < TOL then 
                begin W:= A[M1,M1] - S;
                    if ABS(W) < ABS(SHIFT) then SHIFT:= W
                end;
                A[Q,Q]:= A[Q,Q] - SHIFT;
                for I:= Q step 1 until M1 do 
                begin I1:= I + 1; A[I1,I1]:= A[I1,I1] - SHIFT;
                    KAPPA:= SQRT(A[I,I] ** 2 + A[I1,I] ** 2);
                    if I > Q then 
                    begin A[I,I - 1]:= KAPPA * NU;
                        W:= KAPPA * MU
                    end 
                    else W:= KAPPA; MU:= A[I,I] / KAPPA;
                    NU:= A[I1,I] / KAPPA; A[I,I]:= W;
                    ROTROW(I1, N, I, I1, A, MU, NU);
                    ROTCOL(1, I, I, I1, A, MU, NU);
                    A[I,I]:= A[I,I] + SHIFT;
                    ROTCOL(1, N, I, I1, VEC, MU, NU)
                end;
                A[M,M1]:= A[M,M] * NU; A[M,M]:= A[M,M] * MU + SHIFT
            end 
        end;
        if M > 0 then goto IN;
        for J:= N step -1 until 2 do 
        begin TF[J]:= 1; T:= A[J,J];
            for I:= J - 1 step -1 until 1 do 
            begin DELTA:= T - A[I,I];
                TF[I]:= MATVEC(I + 1, J, I, A, TF) /
                (if ABS(DELTA) < MACHTOL then MACHTOL else DELTA)
            end;
            for I:= 1 step 1 until N do 
            VEC[I,J]:= MATVEC(1, J, I, VEC, TF)
        end;
     END: EM[3]:= ELMAX; EM[5]:= COUNT; REAQRI:= M
    end REAQRI

        eop