program gaus; { -> 75 }
{ pascal program to perform simultaneous solution by Gaussian elimination }
{ procedure GAUSS is included }

const maxr = 8;
maxc = 8;

type ary = array[1..maxr] of real;
arys = array[1..maxc] of real;
ary2s = array[1..maxr,1..maxc] of real;

var y : arys;
coef : arys;
a : ary2s;
n,m : integer;
first,
error : boolean;

external procedure cls;

procedure get_data(var a: ary2s;
var y: arys;
var n,m: integer);

{ get values for n and arrays a,y }

var i,j : integer;

begin
writeln;
repeat
write('How many equations? ');
if not first then cls else first:=false;
m:=n
until n if n>1 then
begin
for i:=1 to n do
begin
writeln('Equation',i:3);
for j:=1 to n do
begin
write(j:3,':');
end;
write(',C:');
end;
writeln;
for i:=1 to n do
begin
for j:=1 to m do
write(a[i,j]:7:4);
writeln(':',y[i]:7:4)
end;
writeln
end { if n>1 }
end; { procedure get_data}

procedure write_data;
{ print out the answeres }

var i : integer;

begin
for i:=1 to m do
write(coef[i]:9:5);
writeln
end; { write_data }

procedure gauss
(a : ary2s;
y : arys;
var coef : arys;
ncol : integer;
var error : boolean);

{ matrix solution by Gaussian Elimination }

var
b : ary2s; { work array, nrow,ncol }
w : arys; { work array, ncol long }
i,j,i1,k,
l,n : integer;
hold,sum,
t,ab,big: real;

begin
error:=false;
n:=ncol;
for i:=1 to n do
begin { copy to work arrays }
for j:=1 to n do
b[i,j]:=a[i,j];
w[i]:=y[i]
end;
for i:=1 to n-1 do
begin
big:=abs(b[i,i]);
l:=i;
i1:=i+1;
for j:=i1 to n do
begin { search for largest element }
ab:=abs(b[j,i]);
if ab>big then
begin
big:=ab;
l:=j
end
end;
if big=0.0 then error:= true
else
begin
if l<>i then
begin
{ interchange rows to put largest element on diagonal }
for j:=1 to n do
begin
hold:=b[l,j];
b[l,j]:=b[i,j];
b[i,j]:=hold
end;
hold:=w[l];
w[l]:=w[i];
w[i]:=hold
end; { if l<>i }
for j:=i1 to n do
begin
t:=b[j,i]/b[i,i];
for k:=i1 to n do
b[j,k]:=b[j,k]-t*b[i,k];
w[j]:=w[j]-t*w[i]
end { j-loop }
end { if big }
end; { i-loop }
if b[n,n]=0.0 then error:=true
else
begin
coef[n]:=w[n]/b[n,n];
i:=n-1;
{ back substitution }
repeat
sum:=0.0;
for j:=i+1 to n do
sum:=sum+b[i,j]*coef[j];
coef[i]:=(w[i]-sum)/b[i,i];
i:=i-1
until i=0
end; { if b[n,n]=0 }
if error then writeln(chr(7),'ERROR: Matrix is singular')
end; { GAUSS }

begin { MAIN }
first:=true;
cls;
writeln;
writeln('Simultaneous solution by Gauss elimination');
repeat
get_data(a,y,n,m);
if n>1 then
begin
gauss(a,y,coef,n,error);
if not error then write_data
end
until n<2
end.
