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