Matrix inversion algorithm

[Domagoj Karnicnik]

Hello everybody!

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 ]

[Timo Salmi]

In article <33a2df01...@news.tel.hr>,
Domagoj Karnicnik <domagoj....@po.tel.hr> wrote:
: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

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

[Mark Vaughan]

In article <33a2df01...@news.tel.hr>, domagoj....@po.tel.hr (Domagoj Karnicnik) wrote:
]-Hello everybody!
]-
]-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.
]-

try the Numerical Methods in Pascal web page at

ftp://cloudlab.larc.nasa.gov/pnumsrc.htm

Mark Vaughan

[Russell Schulz]

Archive-name: invert.pas
Submitted-By: Russell...@locutus.ofB.ORG (Russell Schulz)

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