code 34161; comment MCA 2321; integer procedure QRISYMTRI(A, N, D, B, BB, EM); value N; integer N; array A, D, B, BB, EM; begin integer I, J, J1, K, M, M1, COUNT, MAX; real BBMAX, R, S, SIN, T, C, COS, OLDCOS, G, P, W, TOL, TOL2, LAMBDA, DK1, A0, A1; TOL:= EM[2] * EM[1]; TOL2:= TOL * TOL; COUNT:= 0; BBMAX:= 0; MAX:= EM[4]; M:= N; IN: K:= M; M1:= M - 1; NEXT: K:= K - 1; if K > 0 then begin if BB[K] >= TOL2 then goto NEXT; if BB[K] > BBMAX then BBMAX:= BB[K] end; if K = M1 then M:= M1 else begin T:= D[M] - D[M1]; R:= BB[M1]; if ABS(T) < TOL then S:= SQRT(R) else begin W:= 2 / T; S:= W * R / (SQRT(W * W * R + 1) + 1) end; if K = M - 2 then begin D[M]:= D[M] + S; D[M1]:= D[M1] - S; T:= - S / B[M1]; R:= SQRT(T * T + 1); COS:= 1 / R; SIN:= T / R; ROTCOL(1,N,M1,M,A,COS,SIN); M:= M - 2 end else begin COUNT:= COUNT + 1; if COUNT > MAX then goto END; LAMBDA:= D[M] + S; if ABS(T) < TOL then begin W:= D[M1] - S; if ABS(W) < ABS(LAMBDA) then LAMBDA:= W end; K:= K + 1; T:= D[K] - LAMBDA; COS:= 1; W:= B[K]; P:= SQRT(T * T + W * W); J1:= K; for J:= K + 1 step 1 until M do begin OLDCOS:= COS; COS:= T / P; SIN:= W / P; DK1:= D[J] - LAMBDA; T:= OLDCOS * T; D[J1]:= (T + DK1) * SIN * SIN + LAMBDA + T; T:= COS * DK1 - SIN * W * OLDCOS; W:= B[J]; P:= SQRT(T * T + W * W); G:= B[J1]:= SIN * P; BB[J1]:= G * G; ROTCOL(1, N, J1, J, A, COS, SIN); J1:= J end; D[M]:= COS * T + LAMBDA; if T < 0 then B[M1]:= - G end QRSTEP end; if M > 0 then goto IN; END: EM[3]:= SQRT(BBMAX); EM[5]:= COUNT; QRISYMTRI:= M end QRISYMTRI; eop