Category : Files from Magazines
Archive   : MATRIX2.ZIP
Filename : MATRIX2.PAS

{®RM80¯}
program invertmatrix;
{$N+}

uses crt, dos;

{see Byte,Apr,1986,®MDUL¯Inverting Large Matrices®MDNM¯}

type
stringtype = string[25];
matrix = array[1..32,1..32] of double;
var
im,m1,m2,e : matrix; {matrices, product, and inverse matrix}
i,j,nr,nc,repets,numiter:integer; {max values of indices}
checksolution, insolution, outsolution:matrix;
ch:char;
thld,err:double;
io,disk:boolean;
infile,outfile:text;
infilestr,outfilestr:string;


procedure matrixmultiply (var p,a,b:matrix;nr1,nc1nr2,nc2:integer);
var
i,j,jj :integer;
temp:double;
begin {procedure matrixmultiply}
for i:=1 to nr1 do begin
for j:=1 to nc2 do begin
temp:=0;
for jj:=1 to nc1nr2 do temp:=temp+a[i,jj]*b[jj,j];
p[i,j]:=temp;
end {for j};
end {for i};
end; {procedure matrixmultiply}

procedure matrixadd(var s,a,b:matrix; nc,nr:integer);
var
i,j:integer;
begin
for i:=1 to nr do begin
for j:=1 to nc do s[i,j]:=a[i,j]+b[i,j];
end; {for i}
end; {procedure}

procedure matrixsubtract(var d,p,n:matrix; nr,nc:integer);
var
i,j:integer;
begin
for i:=1 to nr do begin
for j:=1 to nc do d[i,j]:=p[i,j]-n[i,j]; {n=neg; p=pos}
end; {for i}
end; {procedure}

procedure matrixtranspose(var t,m:matrix; nr,nc:integer);
var
i,j:integer;
begin
for i:=1 to nr do begin
for j:=1 to nc do t[i,j]:=m[j,i];
end; {for i}
end; {procedure}

procedure rowsum(var sum:double; var m:matrix; r,nc:integer);
var
j:integer;
temp:double;
begin
temp:=0;
for j:=1 to nc do temp:=temp+abs(m[r,j]);
sum:=temp;
end; {procedure}

procedure colsum(var sum:double; var m:matrix; c,nr:integer);
var
i:integer;
temp:double;
begin
temp:=0;
for i:=1 to nr do temp:=temp+abs(m[i,c]);
sum:=temp;
end; {procedure}

function maxerr(var m:matrix; nc,nr:integer; limit:double):double;
var
temp:double;
i,j:integer;
begin {function maxerr}
temp:=0;
i:=0;
repeat
i:=i+1;
j:=0;
repeat
j:=j+1;
if abs(m[i,j]) > temp then temp:=abs(m[i,j]);
until (temp > limit) or (j >= nc);
until (temp > limit) or (i >= nr);
{limit stops the search when it's uninteresting}
maxerr:=temp;
end; {function maxerr}

procedure initialize(var bb,m:matrix; nr,nc:integer);

{
Find t1=maximum sum of row element magnitudes over all rows.
Find t2=maximum sum of column element magnitudes over all columns.
t=1/(t1*t2).
bb=m(transpose)*t.
}
var
t1,t2,t,temp:double;
i,j:integer;
begin {procedure initialize}
t1:=0;
t2:=0;
for i:=1 to nr do begin
rowsum(temp,m,i,nc);
if temp > t1 then t1:=temp;
end; {for i}
for j:=1 to nc do begin
colsum(temp,m,j,nr);
if temp > t2 then t2:=temp;
end; {for j}
t:=1/(t1*t2);
matrixtranspose(bb,m,nr,nc);
for i:=1 to nr do begin
for j:=1 to nc do bb[i,j]:=bb[i,j]*t;
end; {for i}
end; {procedure initialize}


procedure iterate(var bb,m,e:matrix; nr,nc:integer);
{
Let E(n)=I-BB(n-1)*M (starting with B(0) from initialization).
Let BB(n)=(I+E(n))*BB(n-1) (Newton's iteration).
}
var
tm:matrix;
i,j,n:integer;
begin
matrixmultiply(tm,bb,m,nr,nc,nr);
matrixsubtract(e,im,tm,nr,nc); {Gives error,E}
matrixadd(tm,im,e,nr,nc);
matrixmultiply(bb,tm,bb,nr,nc,nr); {update BB, the estimate of the inverse}
end; {procedure iterate}


procedure enter;
var
i,j:integer;
cch:char;
begin
clrscr;
disk:=false;
io:=true;
cch:=#0;
gotoxy(5,2);
writeln
('PROGRAM TO SOLVE SIMULTANEOUS EQUATIONS USING PAN & REIF MATRIX INVERSION');
gotoxy(15,4);
writeln
('See Byte, April, 1986');
gotoxy(5,6);
write('Enter number of equations and variables (max of 32): ');
readln(nr);
nc:=nr;
gotoxy(5,8);
writeln
('[A] Enter/read data manually; or [B] read from/write to disk? Enter A or B');
cch:=upcase(readkey);
clrscr;

if cch = 'B' then begin
disk:=true;
io:=true;
repeat
gotoxy(5,5); write('Enter input path\filname: ');clreol;
readln(infilestr);
assign(infile,infilestr);
{$I-}
reset(infile);
{$I+}
if ioresult <> 0 then begin
io:=false;
gotoxy(1,25);
clreol;
write(#7,'File not found. Press Q to quit or C to continue.');
cch:=Upcase(readkey);
gotoxy(1,25);clreol;
if cch = 'Q' then halt;
end {if ioresult <>0}
else io:=true;
until io = true;
repeat
gotoxy(5,8); write('Enter output path\filename: ');clreol;
readln(outfilestr);
assign(outfile,outfilestr);
{$I-}
rewrite(outfile);
{$I+}
if ioresult <> 0 then begin
io:=false;
gotoxy(1,25);
clreol;
write(#7,'File not found. Press Q to quit or C to continue.');
cch:=Upcase(readkey);
gotoxy(1,25);clreol;
if cch = 'Q' then halt;
end {if ioresult <>0}
else io:=true;
until io = true;
end;{if cch = B}

if disk = false then begin
clrscr;
for i:= 1 to nr do begin
writeln
('Type equation #',i,' coeff.`s and solution separated by 1 space; then ENTER.');
for j:=1 to nc do read(m1[i,j]);
readln(insolution[i,1]);
end;{for i}
end;{if disk = false}

if disk = true then begin
clrscr;
writeln('INPUT EQUATIONS');
for i:=1 to nr do begin
for j:=1 to nc do begin
{$I-}
read(infile,m1[i,j]);
{$I+}
if ioresult <> 0 then begin
clrscr;
gotoxy(1,25);
write('File not found or bad file format');
halt;
end;{if ioresult <> 0}
write(m1[i,j]:10:2);
end;{for j}
{$I-}
readln(infile,insolution[i,1]);
{$I+}
if ioresult <> 0 then begin
clrscr;
gotoxy(1,25);
write('File not found or bad file format');
halt;
end;{if ioresult <> 0}
writeln(insolution[i,1]:10:2);
end;{for i}
if wherex > 1 then gotoxy(1,wherey+1);
write('Press a key.');
cch:=readkey;cch:=#0;
end;{if disk = true}
end; {procedure enter}


begin {program invmatrix.pas}
for i:=1 to 32 do begin
for j:=1 to 32 do begin
im[i,j]:=0.0;
if i=j then im[i,j]:=1.0;
m1[i,j]:=0.0;
m2[i,j]:=0.0;
e[i,j]:=0.0;
insolution[i,j]:=0.0;
outsolution[i,j]:=0.0;
checksolution[i,j]:=0.0;
end;{for j}
end;{for i}
ch:=#0;
thld:=1.0e-8;
repeat
clrscr;
if ch = 'R' then begin
ch:=#0;
gotoxy(5,8);clreol;
write
('[A] Iterate same solution; or [B] enter new equations? Enter A or B');
ch:=upcase(readkey);
if ch = 'A' then begin
clrscr;
writeln('INPUT EQUATIONS');
for i:= 1 to nr do begin
gotoxy(1,wherey);
for j:=1 to nc do write(m1[i,j]:10:2);
writeln(insolution[i,1]:10:2);
end;{for i};
if wherex > 1 then gotoxy(1,wherey+1);
writeln('Press a key.');
ch:=readkey;
ch:='A';
end;{if ch = A}
end;{if ch = R}
if ch <> 'A' then begin
ch:=#0;
numiter:=0;
enter;
gotoxy(5,25);
clreol;
write('...................Initializing.........................');
initialize(m2,m1,nr,nc);
end; {if ch <> A}
gotoxy(1,25);
clreol;
write('....................Computing...........................');
repets:=20;
if ch = 'A' then repets:=5;
{On first pass do 20 iterations; on successive passes do 5.}

for i:= 1 to repets do iterate(m2,m1,e,nr,nc);
{initializes [e] to (hopefully) valid 1st value.}

numiter:=numiter+repets;


repeat
numiter:=numiter+1;
err:=maxerr(e,nr,nc,1.0); {ensures maxerr is the largest in [e].}
iterate(m2,m1,e,nr,nc);
until maxerr(e,nr,nc,1.0) >= err;
{This code keeps the iteration going as long it helps.}

if err > thld then begin
gotoxy(1,25); clreol;
write
('Significant solution error!',#7,' Max. error in I-M*MI = ',err:14:4);
write('. Press a key.');
ch:=readkey;ch:=#0;
end;{if err > thld}
matrixmultiply(outsolution,m2,insolution,nc,nc,1);
clrscr;
writeln('INVERSE MATRIX');
for i:=1 to nr do begin
for j:=1 to nc do write(m2[i,j]:10:4);
writeln;
end;
writeln('No. iterations = ',numiter:4);
writeln('Max. error in I-M*MI = ',err);
for i:=1 to nc do begin
write(' X',i,' = ',outsolution[i,1]:10:4);
if (i mod 3) = 0 then writeln;
end;{for i}
if wherey = 25 then writeln else
if wherex > 1 then gotoxy(1,wherey+1);
write('Press any key');
ch:=readkey;ch:=#0;
writeln('....For a check the input solution vector computes as...');
matrixmultiply(checksolution,m1,outsolution,nc,nc,1);
for i:=1 to nc do begin
write(checksolution[i,1]:11:4,' ');
if (i mod 4) = 0 then writeln;
end;{for i}

if disk = true then begin
{$I-}
writeln(outfile,'INVERTED MATRIX');
{$I+}
if ioresult <> 0 then begin
clrscr;
gotoxy(1,25);
write('File not found.');
halt;
end;{if ioresult <> 0}
for i:=1 to nr do begin
for j:=1 to nc do write(outfile,m2[i,j]);
writeln(outfile);
end; {for i}
writeln(outfile,'OUTPUT SOLUTION VECTOR');
for i:=1 to nr do write(outfile,outsolution[i,1]);
writeln(outfile,'No. iterations = ',numiter:4);
writeln(outfile,'Max. error in I-M*MI = ',err);
writeln(outfile,'CHECK SOLUTION OF INPUT SOLUTION VECTOR');
for i:=1 to nr do write(outfile,checksolution[i,1]);
writeln(outfile);
end;{if disk = true}

gotoxy(1,25);clreol;
write(' ....Press Q to quit; R to repeat.............');
ch:=upcase(readkey);
until ch = 'Q';
clrscr;
if disk = true then begin
close(infile);
close(outfile);
end;
end.