code 34152;
    comment MCA 2312;
    procedure VECSYMTRI(D, B, N, N1, N2, VAL, VEC, EM);
    value N, N1, N2;
    integer N, N1, N2; array D, B, VAL, VEC, EM;
    begin integer I, J, K, COUNT, MAXCOUNT, COUNTLIM, ORTH, IND;
        real BI, BI1, U, W, Y, MI1, LAMBDA, OLDLAMBDA, ORTHEPS,
        VALSPREAD, SPR, RES, MAXRES, OLDRES, NORM, NEWNORM, OLDNORM,
        MACHTOL, VECTOL;
        array M, P, Q, R, X[1:N];
        boolean array INT[1:N];

        NORM:= EM[1]; MACHTOL:= EM[0] * NORM; VALSPREAD:= EM[4] * NORM;
        VECTOL:= EM[6] * NORM; COUNTLIM:= EM[8]; ORTHEPS:= SQRT(EM[0]);
        MAXCOUNT:= IND:= 0; MAXRES:= 0;
        if N1 > 1 then 
        begin ORTH:= EM[5]; OLDLAMBDA:= VAL[N1 - ORTH];
            for K:= N1 - ORTH + 1 step 1 until N1 - 1  do 
            begin LAMBDA:= VAL[K]; SPR:= OLDLAMBDA - LAMBDA;
                if SPR < MACHTOL then LAMBDA:= OLDLAMBDA - MACHTOL;
                OLDLAMBDA:= LAMBDA
            end 
        end else ORTH:= 1;
        for K:= N1 step 1 until N2 do 
        begin LAMBDA:= VAL[K]; if K > 1 then 
            begin SPR:= OLDLAMBDA - LAMBDA;
                if SPR < VALSPREAD then 
                begin if SPR < MACHTOL then 
                    LAMBDA:= OLDLAMBDA - MACHTOL;
                    ORTH:= ORTH +1
                end else ORTH:= 1
            end;
            COUNT:= 0; U:= D[1] - LAMBDA; BI:= W:= B[1];
            if ABS(BI) < MACHTOL then BI:= MACHTOL;
            for I:= 1 step 1 until N - 1  do 
            begin BI1:= B[I + 1];
                if ABS(BI1) < MACHTOL then BI1:= MACHTOL;
                if ABS(BI) >= ABS(U) then 
                begin MI1:= M[I + 1]:= U / BI; P[I]:= BI;
                    Y:= Q[I]:= D[I + 1] - LAMBDA; R[I]:= BI1;
                    U:= W - MI1 * Y; W:= - MI1 * BI1; INT[I]:= true 
                end 
                else 
                begin MI1:= M[I + 1]:= BI / U; P[I]:= U; Q[I]:= W;
                    R[I]:= 0; U:= D[I + 1] - LAMBDA - MI1 * W;W:= BI1;
                    INT[I]:= false 
                end;
                X[I]:= 1; BI:= BI1
            end TRANSFORM
            P[N]:= if ABS(U) < MACHTOL then MACHTOL else U;
            Q[N]:= R[N]:= 0; X[N]:= 1; goto ENTRY;
        ITERATE: W:= X[1];
            for I:= 2 step 1 until N do 
            begin if INT[I - 1] then 
                begin U:= W; W:= X[I - 1]:= X[I] end 
                else U:= X[I]; W:= X[I]:= U - M[I] * W
            end ALTERNATE;
        ENTRY: U:= W:= 0;
            for I:= N step -1 until 1 do 
            begin Y:= U; U:= X[I]:= (X[I] - Q[I] * U - R[I] * W) /
                P[I]; W:= Y
            end NEXT ITERATION;
            NEWNORM:= SQRT(VECVEC(1, N, 0, X, X)); if ORTH > 1then 
            begin OLDNORM:= NEWNORM;
                for J:= K - ORTH + 1 step 1 until K - 1 do 
                ELMVECCOL(1, N, J, X, VEC, -TAMVEC(1, N, J, VEC, X));
                NEWNORM:= SQRT(VECVEC(1, N, 0, X, X));
                if NEWNORM < ORTHEPS * OLDNORM then 
                begin IND:= IND + 1; COUNT:= 1;
                    for I:= 1 step 1 until IND - 1,
                    IND + 1 step 1 until N do X[I]:= 0;
                    X[IND]:= 1; if IND = N then IND:= 0;
                    goto ITERATE
                end NEW START
            end ORTHOGONALISATION;
            RES:= 1 / NEWNORM; if RES > VECTOL or COUNT = 0 then 
            begin COUNT:= COUNT + 1; if COUNT <= COUNTLIM then 
                begin for I:= 1 step 1 until N do 
                    X[I]:= X[I] * RES; goto ITERATE
                end 
            end;
            for I:= 1 step 1 until N do VEC[I,K]:= X[I] * RES;
            if COUNT > MAXCOUNT then MAXCOUNT:= COUNT;
            if RES > MAXRES then MAXRES:= RES; OLDLAMBDA:= LAMBDA
        end;
        EM[5]:= ORTH; EM[7]:= MAXRES; EM[9]:= MAXCOUNT
    end VECSYMTRI

        eop