Article 13495 of alt.sources:
Archive-name: invert.pas
Submitted-By: Russell_Schulz@locutus.ofB.ORG (Russell Schulz)

domagoj.karnicnik@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_Schulz@locutus.ofB.ORG  Shad 86c


