Matrix inversion algorithm
Author Message
Matrix inversion algorithm

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  ]

Wed, 01 Dec 1999 03:00:00 GMT
Matrix inversion algorithm

: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

....................................................................

Moderating at ftp:// & http://garbo.uwasa.fi archives  193.166.120.5
Department of Accounting and Business Finance  ; University of Vaasa

Wed, 01 Dec 1999 03:00:00 GMT
Matrix inversion algorithm

]-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

Fri, 03 Dec 1999 03:00:00 GMT
Matrix inversion algorithm

Archive-name: invert.pas

Quote:

> 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

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;

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
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;

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;

var
acol: integer;

begin

for acol := 1 to numrows do
begin
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);
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
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;
display(matrix);
solve;
writeln('identity=');
display(matrix);
writeln('inverse=');
display(inverse);
end.
--

Wed, 08 Dec 1999 03:00:00 GMT

 Page 1 of 1 [ 4 post ]

Relevant Pages