code 34181;
comment MCA 2411;
procedure REAVECHES(A, N, LAMBDA, EM, V); value N, LAMBDA;
integer N; real LAMBDA; array A, EM, V;
begin integer I, I1, J, COUNT, MAX;
real M, R, NORM, MACHTOL, TOL;
boolean array P[1:N];
NORM:= EM[1]; MACHTOL:= EM[0] * NORM; TOL:= EM[6] * NORM;
MAX:= EM[8]; A[1,1]:= A[1,1] - LAMBDA;
GAUSS: for I:= 1 step 1 until N - 1 do
begin I1:= I + 1; R:= A[I,I]; M:= A[I1,I];
if ABS(M) < MACHTOL then M:= MACHTOL;
P[I]:= ABS(M) <= ABS(R);
if P[I] then
begin A[I1,I]:= M:= M / R;
for J:= I1 step 1 until N do
A[I1,J]:= (if J > I1 then A[I1,J]
else A[I1,J] - LAMBDA) - M * A[I,J]
end
else
begin A[I,I]:= M; A[I1,I]:= M:= R / M;
for J:= I1 step 1 until N do
begin R:= (if J > I1 then A[I1,J] else
A[I1,J] - LAMBDA);
A[I1,J]:= A[I,J] - M * R; A[I,J]:= R
end
end
end GAUSS;
if ABS(A[N,N]) < MACHTOL then A[N,N]:= MACHTOL;
for J:= 1 step 1 until N do V[J]:= 1; COUNT:= 0;
FORWARD: COUNT:= COUNT + 1; if COUNT > MAX then goto OUT;
for I:= 1 step 1 until N - 1 do
begin I1:= I + 1;
if P[I] then V[I1]:= V[I1] - A[I1,I] * V[I] else
begin R:= V[I1]; V[I1]:= V[I] - A[I1,I] * R;
V[I]:=R
end
end FORWARD;
BACKWARD: for I:= N step -1 until 1 do
V[I]:= (V[I] - MATVEC(I + 1, N, I, A, V)) / A[I,I];
R:= 1 / SQRT(VECVEC(1, N, 0, V, V));
for J:= 1 step 1 until N do V[J]:= V[J] * R;
if R > TOL then goto FORWARD;
OUT: EM[7]:= R; EM[9]:= COUNT
end REAVECHES
eop