code 34190; comment MCA 2420; integer procedure COMVALQRI(A, N, EM, RE, IM); value N; integer N; array A, EM, RE, IM; begin integer I, J, P, Q, MAX, COUNT, N1, P1, P2, IMIN1, I1, I2, I3; real DISC, SIGMA, RHO, G1, G2, G3, PSI1, PSI2, AA, E, K, S, NORM, MACHTOL2, TOL, W; boolean B; NORM:= EM[1]; MACHTOL2:= (EM[0] * NORM) ** 2; TOL:= EM[2] * NORM; MAX:= EM[4]; COUNT:= 0; W:= 0; IN: for I:= N, I - 1 while (if I >= 1 then ABS(A[I + 1,I]) > TOL else false ) do Q:= I; if Q > 1 then begin if ABS(A[Q,Q - 1]) > W then W:= ABS(A[Q,Q - 1]) end; if Q >= N - 1 then begin N1:= N - 1; if Q = N then begin RE[N]:= A[N,N]; IM[N]:= 0; N:= N1 end else begin SIGMA:= A[N,N] - A[N1,N1]; RHO:= -A[N,N1] * A[N1,N]; DISC:= SIGMA ** 2 - 4 * RHO; if DISC > 0 then begin DISC:= SQRT(DISC); S:= -2 * RHO / (SIGMA + (if SIGMA >= 0 then DISC else -DISC)); RE[N]:= A[N,N] + S; RE[N1]:= A[N1,N1] - S; IM[N]:= IM[N1]:= 0 end else begin RE[N]:= RE[N1]:= (A[N1,N1] + A[N,N]) / 2; IM[N1]:= SQRT( -DISC) / 2; IM[N]:= -IM[N1] end; N:= N - 2 end end else begin COUNT:= COUNT + 1; if COUNT > MAX then goto OUT; N1:= N - 1; SIGMA:= A[N,N] + A[N1,N1] + SQRT(ABS(A[N1,N - 2] * A[N,N1]) * EM[0]); RHO:= A[N,N] * A[N1,N1] - A[N,N1] * A[N1,N]; for I:= N - 1, I - 1 while (if I - 1 >= Q then ABS(A[I,I - 1] * A[I1,I] * (ABS(A[I,I] + A[I1,I1] - SIGMA) + ABS(A[I + 2,I1]))) > ABS(A[I,I] * ((A[I,I] - SIGMA) + A[I,I1] * A[I1,I] + RHO)) * TOL else false ) do P1:= I1:= I; P:= P1 - 1; P2:= P + 2; for I:= P step 1 until N - 1 do begin IMIN1:= I - 1; I1:= I + 1; I2:= I + 2; if I = P then begin G1:= A[P,P] * (A[P,P] - SIGMA) + A[P,P1] * A[P1,P] + RHO; G2:= A[P1,P] * (A[P,P] + A[P1,P1] - SIGMA); if P1 <= N1 then begin G3:= A[P1,P] * A[P2,P1]; A[P2,P]:= 0 end else G3:= 0 end else begin G1:= A[I,IMIN1]; G2:= A[I1,IMIN1]; G3:= if I2 <= N then A[I2,IMIN1] else 0 end; K:= if G1 >= 0 then SQRT(G1 ** 2 + G2 ** 2 + G3 ** 2) else -SQRT(G1 ** 2 + G2 ** 2 + G3 ** 2); B:= ABS(K) > MACHTOL2; AA:= if B then G1 / K + 1 else 2; PSI1:= if B then G2 / (G1 + K) else 0; PSI2:= if B then G3 / (G1 + K) else 0; if I ^= Q then A[I,IMIN1]:= if I = P then -A[I,IMIN1] else -K; for J:= I step 1 until N do begin E:= AA * (A[I,J] + PSI1 * A[I1,J] + (if I2 <= N then PSI2 * A[I2,J] else 0)); A[I,J]:= A[I,J] - E; A[I1,J]:= A[I1,J] - PSI1 * E; if I2 <= N then A[I2,J]:= A[I2,J] - PSI2 * E end; for J:= Q step 1 until (if I2 <= N then I2 else N) do begin E:= AA * (A[J,I] + PSI1 * A[J,I1] + (if I2 <= N then PSI2 * A[J,I2] else 0)); A[J,I]:= A[J,I] - E; A[J,I1]:= A[J,I1] - PSI1 * E; if I2 <= N then A[J,I2]:= A[J,I2] - PSI2 * E end; if I2 <= N1 then begin I3:= I + 3; E:= AA * PSI2 * A[I3,I2]; A[I3,I]:= -E; A[I3,I1]:= -PSI1 * E; A[I3,I2]:= A[I3,I2] - PSI2 * E end end end; if N > 0 then goto IN; OUT: EM[3]:= W; EM[5]:= COUNT; COMVALQRI:= N end COMVALQRI eop