code 34160;
    integer procedure QRIVALSYMTRI(D, BB, N, EM); value N;
    integer N; array D, BB, EM;
    begin integer I, I1, LOW, OLDLOW, N1, COUNT, MAX;
        real BBTOL, BBMAX, BBI, BBN1, MACHTOL, DN, DELTA, F, NUM,
        SHIFT, G, H, T, P, R, S, C, OLDG;
        BBTOL:= (EM[2] * EM[1]) ** 2; MACHTOL:= EM[0] * EM[1];
        MAX:= EM[4]; BBMAX:= 0; COUNT:= 0; OLDLOW:= N;
        for N1:= N - 1 while N > 0 do 
        begin 
            for I:= N, I - 1 while (if I >= 1 then 
            BB[I] > BBTOL else false ) do LOW:= I;
            if LOW > 1 then begin if BB[LOW-1] > BBMAX then 
            BBMAX:= BB[LOW-1] end;
            if LOW = N then N:= N1 else 
            begin DN:= D[N]; DELTA:= D[N1] - DN;
                BBN1:= BB[N1];
                if ABS(DELTA) < MACHTOL then R:= SQRT(BBN1) else 
                begin 
                    F:= 2 / DELTA; NUM:= BBN1 * F;
                    R:= -NUM / (SQRT(NUM * F + 1) + 1)
                end;
                if LOW = N1 then 
                begin D[N]:= DN + R; D[N1]:= D[N1] - R; N:= N - 2
                end 
                else 
                begin COUNT:= COUNT + 1;
                    if COUNT > MAX then goto END;
                    if LOW < OLDLOW then 
                    begin SHIFT:= 0; OLDLOW:= LOW end 
                    else  SHIFT:= DN + R;
                    H:= D[LOW] - SHIFT;
                    if ABS(H) < MACHTOL then H:= if H <= 0 then 
                    -MACHTOL else MACHTOL;
                    G:= H; T:= G * H;
                    BBI:= BB[LOW]; P:= T + BBI; I1:= LOW;
                    for I:= LOW + 1 step 1 until N do 
                    begin S:= BBI / P; C:= T / P;
                        H:= D[I] - SHIFT - BBI / H;
                        if ABS(H) < MACHTOL then H:= if H <= 0
                        then -MACHTOL else MACHTOL;
                        OLDG:= G; G:= H * C; T:= G * H;
                        D[I1]:= OLDG - G + D[I];
                        BBI:= if I = N then 0 else BB[I];
                        P:= T + BBI; BB[I1]:= S * P; I1:= I
                    end;
                    D[N]:= G + SHIFT
                end QRSTEP
            end 
        end;
     END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRIVALSYMTRI:= N
    end QRIVALSYMTRI

         eop