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