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



Wed, 08 Dec 1999 03:00:00 GMT  
 
 [ 4 post ] 

 Relevant Pages 

1. Matrix inversion algorithm

2. linear equations solving / matrix inversion with ill-conditioned matrix?

3. Algorithm for filling a matrix

4. Word & phrase db inversion

5. TP 6.0 algorithms to compute matrixes' determinant and inverse needed

6. Source Code For Matrix Screen Saver In Pascal?

7. matrix problem

8. bvp303d.zip Vector and matrix library for TP7, OptiCode

9. converting byte matrix to graphics file

10. matrix problem

11. inverse matrix

12. printing using DOT MATRIX printer

 

 
Powered by phpBB® Forum Software