code 34191;
comment MCA 2421;
procedure COMVECHES(A, N, LAMBDA, MU, EM, U, V);
value N, LAMBDA, MU;
integer N; real LAMBDA, MU; array A, EM, U, V;
begin integer I, I1, J, COUNT, MAX;
real AA, BB, D, M, R, S, W, X, Y, NORM, MACHTOL, TOL;
array G, F[1:N];
boolean array P[1:N];
NORM:= EM[1]; MACHTOL:= EM[0] * NORM; TOL:= EM[6] * NORM;
MAX:= EM[8];
for I:= 2 step 1 until N do
begin F[I - 1]:= A[I,I - 1]; A[I,1]:= 0 end;
AA:= A[1,1] - LAMBDA; BB:= -MU;
for I:= 1 step 1 until N - 1 do
begin I1:= I + 1; M:= F[I];
if ABS(M) < MACHTOL then M:= MACHTOL;
A[I,I]:= M; D:= AA ** 2 + BB ** 2; P[I]:= ABS(M) < SQRT(D);
if P[I] then
begin comment A[I,J] * FACTOR AND A[I1,J] - A[I,J];
F[I]:= R:= M * AA / D; G[I]:= S:= -M * BB / D;
W:= A[I1,I]; X:= A[I,I1]; A[I1,I]:= Y:= X * S + W * R;
A[I,I1]:= X:= X * R - W * S;
AA:= A[I1,I1] - LAMBDA - X; BB:= -(MU + Y);
for J:= I + 2 step 1 until N do
begin W:= A[J,I]; X:= A[I,J];
A[J,I]:= Y:= X * S + W * R;
A[I,J]:= X:= X * R - W * S; A[J,I1]:= -Y;
A[I1,J]:= A[I1,J] - X
end
end
else
begin comment INTERCHANGE A[I1,J] AND
A[I,J] - A[I1,J] * FACTOR;
F[I]:= R:= AA / M; G[I]:= S:= BB / M;
W:= A[I1,I1] - LAMBDA; AA:= A[I,I1] - R * W - S * MU;
A[I,I1]:= W; BB:= A[I1,I] - S * W + R * MU;
A[I1,I]:= -MU;
for J:= I + 2 step 1 until N do
begin W:= A[I1,J]; A[I1,J]:= A[I,J] - R * W;
A[I,J]:= W;
A[J,I1]:= A[J,I] - S * W; A[J,I]:= 0
end
end
end
P[N]:= true; D:= AA ** 2 + BB ** 2; if D < MACHTOL ** 2
then begin AA:= MACHTOL; BB:= 0; D:= MACHTOL ** 2 end;
A[N,N]:= D; F[N]:= AA; G[N]:= -BB;
for I:= 1 step 1 until N do
begin U[I]:= 1; V[I]:= 0 end;
COUNT:= 0;
FORWARD: if COUNT > MAX then goto OUTM;
for I:= 1 step 1 until N do
begin if P[I] then
begin W:= V[I]; V[I]:= G[I] * U[I] + F[I] * W;
U[I]:= F[I] * U[I] - G[I] * W; if I < N then
begin V[I + 1]:= V[I + 1] - V[I];
U[I + 1]:= U[I + 1] - U[I]
end
end
else
begin AA:= U[I + 1]; BB:= V[I + 1];
U[I + 1]:= U[I] - (F[I] * AA - G[I] * BB); U[I]:= AA;
V[I + 1]:= V[I] - (G[I] * AA + F[I] * BB); V[I]:= BB
end
end FORWARD;
BACKWARD: for I:= N step -1 until 1 do
begin I1:= I + 1;
U[I]:= (U[I] - MATVEC(I1, N, I, A, U) + (if P[I] then
TAMVEC(I1, N, I, A, V) else A[I1,I] * V[I1])) / A[I,I];
V[I]:= (V[I] - MATVEC(I1, N, I, A, V) - (if P[I] then
TAMVEC(I1, N, I, A, U) else A[I1,I] * U[I1])) / A[I,I]
end BACKWARD;
NORMALISE: W:= 1 / SQRT(VECVEC(1, N, 0, U, U) +
VECVEC(1, N, 0, V, V));
for J:= 1 step 1 until N do
begin U[J]:= U[J] * W; V[J]:= V[J] * W end;
COUNT:= COUNT + 1; if W > TOL then goto FORWARD;
OUTM: EM[7]:= W; EM[9]:= COUNT
end COMVECHES;
eop