Const ndim = 10;
ndimS = 10;
ndimC = 10;
maxnfe = 3000;
type
float = extended;
type floatvector = array [1..ndim] of float;
floatmatrix = array [1..ndim,1..ndim] of float;
rvecn = array [1..6*ndim+3] of float;
ivec5 = array [1..5] of integer;
ivector = array [1..ndim] of integer;
floatvectorS = array [1..ndimS] of float;
floatvectorC = array [1..ndimC] of float;
floatmatrixCS = array [1..ndimS,1..ndimC] of float;
procedure DECOMP(
n : integer;
var A : floatmatrix ;
var cond : float ;
var ipvt : ivector;
var work : floatvector
);
VAR
nm1,i,j,k,kp1,kb,km1,m : integer;
ek,t,anorm,ynorm,znorm : float;
begin
ipvt[n]:=1;
if n = 1 then begin
cond:=1.0;
if a[1,1] <> 0.0 then exit
else begin
cond:=1.0e+32;
exit
end
end;
nm1:=n-1;
anorm:=0.0;
for j:=1 to n do begin
t:=0.0;
for i:=1 to n do t:=t+abs(a[i,j]);
if t > anorm then anorm:=t;
end;
{do 35} for k:=1 to nm1 do begin
kp1:=k + 1;
m:=k;
for i:=kp1 to n do
if abs(a[i,k]) > abs(a[m,k]) then m:=i;
ipvt[k]:=m;
if m <> k then ipvt[n]:=-ipvt[n];
t:=a[m,k];
a[m,k]:=a[k,k];
a[k,k]:=t;
if t <> 0.0 then begin
{20} for i:=kp1 to n do a[i,k]:=-a[i,k]/t;
{do 30} for j:=kp1 to n do begin
t:=a[m,j];
a[m,j]:=a[k,j];
a[k,j]:=t;
if t <> 0.0 then for i:=kp1 to n do
{25} a[i,j]:=a[i,j] + a[i,k]*t;
{30} end;
{35} end;
end;
{do 50}for k:=1 to n do begin
t:=0.0;
if k <> 1 then begin
km1:=k-1;
for i:=1 to km1 do
t:=t + a[i,k] * work[i];
end;
{45} ek:=1.0;
if t < 0.0 then ek:=-1.0;
if a[k,k] = 0.0 then begin
cond:=1.0E+32;
exit
end;
work[k]:=-(ek+t) / a[k,k];
{50} end;
{do 60}for kb:=1 to nm1 do begin
k:=n-kb;
t:=0.0;
kp1:=k + 1;
for i:=kp1 to n do t:=t + a[i,k] * work[k];
work[k]:=t;
m:=ipvt[k];
if m <> k then begin
t:=work[m];
work[m]:=work[k];
work[k]:=t
end;
{60} end;
ynorm:=0.0;
for i:=1 to n do ynorm:=ynorm + abs(work[i]);
solve(n,a,work,ipvt);
znorm:=0.0;
for i:=1 to n do znorm:=znorm + abs(work[i]);
cond:=anorm*znorm/ynorm;
if cond < 1.0 then cond:=1.0;
exit
END;
