code 34428; procedure DECSOLTRIPIV(SUB, DIAG, SUPER, N, AUX, B); value N; integer N; array SUB, DIAG, SUPER, AUX, B; begin integer I, I1, N1, N2; real D, R, S, U, T, Q, V, W, NORM, NORM1, NORM2, TOL, BI, BI1, BI2; boolean array PIV[1:N]; TOL:= AUX[2]; D:= DIAG[1]; R:= SUPER[1]; BI:= B[1]; NORM:= NORM2:= ABS(D) + ABS(R); N2:= N - 2; for I:= 1 step 1 until N2 do begin I1:= I + 1; S:= SUB[I]; T:= DIAG[I1]; Q:= SUPER[I1]; BI1:= B[I1]; NORM1:= NORM2; NORM2:= ABS(S) + ABS(T) + ABS(Q); if NORM2 > NORM then NORM:= NORM2; if ABS(D) * NORM2 < ABS(S) * NORM1 then begin if ABS(S) <= TOL * NORM2 then begin AUX[3]:= I - 1; AUX[5]:= S; goto EXIT end; U:= SUPER[I]:= T / S; BI1:= B[I]:= BI1 / S; BI:= BI - BI1 * D; V:= SUB[I]:= Q / S; W:= SUPER[I1]:= -V * D; D:= DIAG[I1]:= R - U * D; R:= W; NORM2:= NORM1; PIV[I]:= true end else begin if ABS(D) <= TOL * NORM1 then begin AUX[3]:= I - 1; AUX[5]:= D; goto EXIT end; U:= SUPER[I]:= R / D; BI:= B[I]:= BI / D; BI:= BI1 - BI * S; D:= DIAG[I1]:= T - U * S; PIV[I]:= false; R:= Q end end; N1:= N - 1; S:= SUB[N1]; T:= DIAG[N]; NORM1:= NORM2; BI1:= B[N]; NORM2:= ABS(S) + ABS(T); if NORM2 > NORM then NORM:= NORM2; if ABS(D) * NORM2 < ABS(S) * NORM1 then begin if ABS(S) <= TOL * NORM2 then begin AUX[3]:= N2; AUX[5]:= S; goto EXIT end; U:= SUPER[N1]:= T / S; BI1:= B[N1]:= BI1 / S; BI:= BI - BI1 * D; D:= R - U * D; NORM2:= NORM1 end else begin if ABS(D) <= TOL * NORM1 then begin AUX[3]:= N2; AUX[5]:= D; goto EXIT end; U:= SUPER[N1]:= R / D; BI:= B[N1]:= BI / D; BI:= BI1 - BI * S; D:= T - U * S end; if ABS(D) <= TOL * NORM2 then begin AUX[3]:= N1; AUX[5]:= D; goto EXIT end; AUX[3]:= N; AUX[5]:= NORM; BI1:= B[N]:= BI / D; BI:= B[N1]:= B[N1] - SUPER[N1] * BI1; for I:= N - 2 step -1 until 1 do begin BI2:= BI1; BI1:= BI; BI:= B[I]:= B[I] - SUPER[I] * BI1 - (if PIV[I] then SUB[I] * BI2 else 0) end; EXIT: end DECSOLTRIPIV; eop