code 34161;
comment MCA 2321;
integer procedure QRISYMTRI(A, N, D, B, BB, EM); value N;
integer N; array A, D, B, BB, EM;
begin integer I, J, J1, K, M, M1, COUNT, MAX;
real BBMAX, R, S, SIN, T, C, COS, OLDCOS, G, P, W, TOL, TOL2,
LAMBDA, DK1, A0, A1;
TOL:= EM[2] * EM[1]; TOL2:= TOL * TOL; COUNT:= 0; BBMAX:= 0;
MAX:= EM[4]; M:= N;
IN: K:= M; M1:= M - 1;
NEXT: K:= K - 1; if K > 0 then
begin if BB[K] >= TOL2 then goto NEXT;
if BB[K] > BBMAX then BBMAX:= BB[K]
end;
if K = M1 then M:= M1 else
begin
T:= D[M] - D[M1]; R:= BB[M1];
if ABS(T) < TOL then S:= SQRT(R) else
begin W:= 2 / T; S:= W * R / (SQRT(W * W * R + 1) + 1)
end; if K = M - 2 then
begin D[M]:= D[M] + S; D[M1]:= D[M1] - S;
T:= - S / B[M1]; R:= SQRT(T * T + 1); COS:= 1 / R;
SIN:= T / R; ROTCOL(1,N,M1,M,A,COS,SIN); M:= M - 2
end
else
begin COUNT:= COUNT + 1;
if COUNT > MAX then goto END;
LAMBDA:= D[M] + S; if ABS(T) < TOL then
begin W:= D[M1] - S;
if ABS(W) < ABS(LAMBDA) then LAMBDA:= W
end;
K:= K + 1; T:= D[K] - LAMBDA; COS:= 1; W:= B[K];
P:= SQRT(T * T + W * W); J1:= K;
for J:= K + 1 step 1 until M do
begin OLDCOS:= COS; COS:= T / P; SIN:= W / P;
DK1:= D[J] - LAMBDA; T:= OLDCOS * T;
D[J1]:= (T + DK1) * SIN * SIN + LAMBDA + T;
T:= COS * DK1 - SIN * W * OLDCOS; W:= B[J];
P:= SQRT(T * T + W * W); G:= B[J1]:= SIN * P;
BB[J1]:= G * G; ROTCOL(1, N, J1, J, A, COS, SIN);
J1:= J
end;
D[M]:= COS * T + LAMBDA; if T < 0 then B[M1]:= - G
end QRSTEP
end;
if M > 0 then goto IN;
END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRISYMTRI:= M
end QRISYMTRI;
eop