I have a little problem. I need to write a piece of code that would do
matrix inversion, but I have no literature where I can find the
algorithm. I remember it's done by writing a matrix with 1s as
diagonal members on the right side, but what follows I don't remember.
I have it written down somewhere in my schoolbooks, but I just can't
find it.
Could anyone send me algorithm or code snippet that solves this
problem?
TIA, Domagoj.
An example:
[ 5 -2 4]
A = [-3 1 -5]
[ 2 -1 3]
[ 5 -2 4 | 1 0 0 ]
-> [-3 1 -5 | 0 1 0 ]
[ 2 -1 3 | 0 0 1]
-> ... ???
Result:
[ 1/2 -1/2 -3/2 ]
Ainv = [ 1/4 -1/4 -13/4 ]
[-1/4 -1/4 1/4 ]
303778 May 2 1991 ftp://garbo.uwasa.fi/pc/turbopas/nrpas13.zip
nrpas13.zip Numerical Recipes Pascal shareware version, 303K!
All the best, Timo
....................................................................
Prof. Timo Salmi Co-moderator of news:comp.archives.msdos.announce
Moderating at ftp:// & http://garbo.uwasa.fi archives 193.166.120.5
Department of Accounting and Business Finance ; University of Vaasa
mailto:t...@uwasa.fi <URL:http://uwasa.fi/~ts> ; FIN-65101, Finland
try the Numerical Methods in Pascal web page at
ftp://cloudlab.larc.nasa.gov/pnumsrc.htm
Mark Vaughan
domagoj....@po.tel.hr (Domagoj Karnicnik) writes:
> I remember it's done by writing a matrix with 1s as diagonal
> members on the right side, but what follows I don't remember. [...]
it's one approach, but there are others. (you can search the web
for `Gaussian elimination', and no doubt find better, more stable
code than this).
you can add a multiple of a row to any other row, or multiply any
row by a constant, or swap any two rows.
I've set followups to c.l.p.misc -- there's nothing especially
Borlandish about this, I think.
program invert; { invert a small matrix }
const
maxsize=20;
epsilon=1e-10;
type
{ for Delphi }
{$ifndef VER90}
{$ifndef VER80}
double=real;
{$endif}
{$endif}
matrixt=array[1..maxsize,1..maxsize] of double;
var
matrix: matrixt;
inverse: matrixt;
numrows: integer;
matrixfn: string;
procedure detab(var astring: string);
var
tempint: integer;
begin
for tempint := 1 to length(astring) do
if astring[tempint]=chr(9) then
astring[tempint] := ' ';
end;
function alltrim(astring: string): string;
var
firstpos: integer;
lastpos: integer;
done: boolean;
begin
firstpos := 1;
lastpos := length(astring);
if astring='' then
alltrim := ''
else
begin
done := false;
while not done do
begin
if firstpos>length(astring) then
done := true
else if (astring[firstpos]<>' ') and (astring[firstpos]<>chr(9)) then
done := true
else
firstpos := firstpos+1;
end;
done := false;
while not done do
begin
if lastpos<1 then
done := true
else if (astring[lastpos]<>' ') and (astring[lastpos]<>chr(9)) then
done := true
else
lastpos := lastpos-1;
end;
if firstpos>length(astring) then
alltrim := ''
else
alltrim := copy(astring,firstpos,lastpos-firstpos+1);
end;
end;
function chopword(var astring: string): string;
var
tempstring: string;
lastpos: integer;
begin
astring := alltrim(astring);
lastpos := pos(' ',astring);
if lastpos=0 then
lastpos := length(astring)
else
lastpos := lastpos-1;
tempstring := copy(astring,1,lastpos);
astring := copy(astring,lastpos+1,255);
chopword := tempstring;
end;
function min(i1,i2: integer): integer;
begin
if i1<i2 then
min := i1
else
min := i2;
end;
procedure usage;
begin
writeln('invert: invert a small matrix');
writeln;
writeln('usage: invert matrix.txt');
writeln('will read in matrix.txt (one row per line) and write its');
writeln('inverse to stdout');
halt(1);
end;
procedure msgusage(amsg: string);
begin
writeln(amsg);
usage;
end;
procedure display(amatrix: matrixt);
var
arow, acol: integer;
begin
for arow := 1 to numrows do
begin
for acol := 1 to numrows do
write(amatrix[arow,acol]:10:2);
if acol=numrows then
writeln
else
write(' ');
end;
end;
function isdiff(a,b: double): boolean;
begin
isdiff := (abs(a-b)>epsilon);
end;
procedure initialize;
begin
if paramcount<>1 then
usage;
matrixfn := paramstr(1);
end;
procedure readin;
var
matrixf: text;
matrixline: string;
numfields: integer;
currfield: integer;
afield: string;
avalue: double;
code: word;
begin
assign(matrixf,matrixfn);
{$I-}
reset(matrixf);
{$I+}
if ioresult<>0 then
msgusage('could not open '+matrixfn);
numrows := 0;
numfields := 0;
while not eof(matrixf) do
begin
readln(matrixf,matrixline);
numrows := numrows+1;
currfield := 1;
matrixline := alltrim(matrixline);
detab(matrixline);
while matrixline<>'' do
begin
afield := chopword(matrixline);
val(afield,avalue,code);
if code<>0 then
begin
close(matrixf);
writeln('"',afield,'" is not a valid number');
msgusage('no matrix found in '+matrixfn);
end;
if (currfield>maxsize) or (numrows>maxsize) then
begin
close(matrixf);
msgusage('recompile with larger MAXSIZE');
end;
matrix[numrows,currfield] := avalue;
inverse[numrows,currfield] := 0;
if numrows=currfield then
inverse[numrows,currfield] := 1;
inc(currfield);
end;
if numrows=1 then
numfields := currfield-1;
if numfields<>currfield-1 then
begin
close(matrixf);
writeln('line ',numrows,' of ',matrixfn,' is inconsistent');
writeln('found field count: ',currfield-1);
writeln('expected field count: ',numfields);
msgusage('no matrix found in '+matrixfn);
end;
end;
close(matrixf);
if numfields<>numrows then
msgusage('not a square matrix in '+matrixfn);
if numfields=0 then
msgusage('no matrix found in '+matrixfn);
writeln('number of rows: ',numrows);
end;
procedure addrows(r1,r2: integer; destrow: integer);
var
acol: integer;
begin
writeln('adding rows ',r1,' and ',r2,' into row ',destrow);
for acol := 1 to numrows do
begin
matrix[destrow,acol] := matrix[r1,acol]+matrix[r2,acol];
inverse[destrow,acol] := inverse[r1,acol]+inverse[r2,acol];
end;
end;
procedure multrow(adouble: double; destrow: integer);
var
acol: integer;
begin
writeln('multiplying row ',destrow,' by ',adouble:10:2);
for acol := 1 to numrows do
begin
matrix[destrow,acol] := matrix[destrow,acol]*adouble;
inverse[destrow,acol] := inverse[destrow,acol]*adouble;
end;
end;
procedure solvecell(arow,acol: integer; desired: double);
var
otherrow: integer;
canuse: boolean;
othercol: integer;
begin
if isdiff(matrix[arow,acol],desired) then
begin
if isdiff(matrix[arow,acol],0) and (desired<>0) then
multrow(1/matrix[arow,acol],arow);
if isdiff(matrix[arow,acol],desired) then
begin
for otherrow := 1 to numrows do
if (otherrow<>arow) and isdiff(matrix[otherrow,acol],0) and
isdiff(matrix[arow,acol],desired) then
begin
canuse := true;
for othercol := 1 to acol-1 do
if isdiff(matrix[otherrow,othercol],0) then
canuse := false;
if canuse then
begin
multrow(-matrix[otherrow,acol]/matrix[arow,acol],arow);
addrows(otherrow,arow,arow);
end;
end;
end;
if isdiff(matrix[arow,acol],desired) then
begin
writeln('can''t solve row=',arow,'; col=',acol,'; value=',
matrix[arow,acol]:10:2,'; desired=',desired:10:2);
writeln('currently, matrix=');
display(matrix);
writeln('and inverse=');
display(inverse);
halt(1);
end;
end;
end;
procedure populatediagonal(arow: integer);
var
otherrow: integer;
begin
writeln('populatediagonal(',arow,');');
if not isdiff(matrix[arow,arow],0) then
begin
for otherrow := 1 to numrows do
if isdiff(matrix[otherrow,arow],0) then
if not isdiff(matrix[arow,arow],0) then
addrows(otherrow,arow,arow);
end;
end;
procedure solveleft(arow: integer);
var
acol: integer;
begin
for acol := 1 to arow-1 do
solvecell(arow,acol,0);
end;
procedure solveright(arow: integer);
var
acol: integer;
begin
for acol := numrows downto arow+1 do
solvecell(arow,acol,0);
end;
procedure solve;
var
arow: integer;
begin
for arow := 1 to numrows do
populatediagonal(arow);
for arow := 1 to numrows do
solveleft(arow);
for arow := 1 to numrows do
solvecell(arow,arow,1);
for arow := numrows downto 1 do
solveright(arow);
for arow := 1 to numrows do
solvecell(arow,arow,1);
end;
begin
initialize;
readin;
display(matrix);
solve;
writeln('identity=');
display(matrix);
writeln('inverse=');
display(inverse);
end.
--
Russell...@locutus.ofB.ORG Shad 86c