real procedure INNERPRODUCT(u, v, k, s, f);
   value s, f; integer k, s, f; real u, v;
     comment;
begin
   real h;
   h ≔ 0; for k ≔ s step 1 until f do h ≔ h + u × v;
   INNERPRODUCT ≔ h
end INNERPRODUCT;
procedure CROUT(A, b, n, y, pivot, INNERPRODUCT);
   value n; array A, b, y; integer n; integer array pivot;
     real procedure INNERPRODUCT;
     comment;
begin
   integer k, i, j, imax, p; real TEMP, quot;
   for k ≔ 1 step 1 until n do
      begin
         TEMP ≔ 0;
         for i ≔ k step 1 until n do
            begin
               A[i, k] ≔ A[i, k] - INNERPRODUCT(A[i,p], A[p,k],
                 p, 1, k - 1);
               if abs(A[i,k]) > TEMP then
               begin
                  TEMP ≔ abs(A[i, k]); imax ≔ i
               end
            end;
         pivot[k] ≔ imax;
         comment;
         if imax ≠ k then begin for j ≔ 1 step 1 until n do
               begin
                  TEMP ≔ A[k,j]; A[k,j] ≔ A[imax,j];
                  A[imax,j] ≔ TEMP
               end;
            TEMP ≔ b[k]; b[k] ≔ b[imax]; b[imax] ≔ TEMP
         end;
         comment;
         if A[k,k] = 0 then goto singular;
         for i ≔ k+1 step 1 until n do
            begin quot ≔ 1·0/A[k,k]; A[i,k] ≔ quot × A[i,k]
            end;
         for j ≔ k + 1 step 1 until n do
              A[k, j] ≔ A[k, j] - INNERPRODUCT(A[k,p],
              A[p,j], p, 1, k - 1);
         b[k] ≔ b[k] - INNERPRODUCT(A[k, p], b[p], p,
           1, k - 1)
      end;
   comment;
   for k ≔ n step -1 until 1 do
        y[k] ≔ (b[k] - INNERPRODUCT(A[k,p], y[p], p, k + 1, n))/A[k, k];
   singular:
end CROUT;

procedure SOLVE(B, c, n, z, pivot, INNERPRODUCT);
   value n; array B, c, z; integer n; integer array pivot;
     real procedure INNERPRODUCT;
     comment;
begin
   integer k, p; real TEMP;
   for k ≔ 1 step 1 until n do
      begin
         TEMP ≔ c[pivot[k]]; c[pivot[k]] ≔ c[k]; c[k] ≔            TEMP; c[k] ≔ c[k] - INNERPRODUCT(B[k, p],
           c[p], p, 1, k - 1)
      end;
   for k ≔ n step -1 until 1 do
        z[k] ≔ (c[k] - INNERPRODUCT(B[k,p], z[p], p,
        k+1, n))/B[k,k]
end;