code 34140;
comment MCA 2300;
procedure TFMSYMTRI2(A, N, D, B, BB, EM); value N;integer N;
array A, B, BB, D, EM;
begin integer I, J, R, R1;
real W, X, A1, B0, BB0, D0, MACHTOL, NORM;
NORM:= 0;
for J:= 1 step 1 until N do
begin W:= 0;
for I:= 1 step 1 until J do W:= ABS(A[I,J]) + W;
for I:= J + 1 step 1 until N do W:= ABS(A[J,I]) +
W; if W > NORM then NORM:= W
end;
MACHTOL:= EM[0] * NORM; EM[1]:= NORM; R:= N;
for R1:= N - 1 step -1 until 1 do
begin D[R]:= A[R,R]; X:= TAMMAT(1, R - 2, R, R, A, A);
A1:= A[R1,R]; if SQRT(X) <= MACHTOL then
begin B0:= B[R1]:= A1; BB[R1]:= B0 * B0;A[R,R]:= 1 end
else
begin BB0:= BB[R1]:= A1 * A1 + X;
B0:= if A1 > 0 then -SQRT(BB0) else SQRT(BB0);
A1:= A[R1,R]:= A1 - B0; W:= A[R,R]:= 1 / (A1 * B0);
for J:= 1 step 1 until R1 do B[J]:= (TAMMAT(1,
J, J, R, A, A) + MATMAT(J + 1, R1, J, R, A, A)) * W;
ELMVECCOL(1, R1, R, B, A, TAMVEC(1, R1, R, A, B) *
W * .5); for J:= 1 step 1 until R1 do
begin ELMCOL(1, J, J, R, A, A, B[J]);
ELMCOLVEC(1, J, J, A, B, A[J,R])
end; B[R1]:= B0
end; R:= R1
end;
D[1]:= A[1,1]; A[1,1]:= 1; B[N]:= BB[N]:= 0
end TFMSYMTRI2
eop