Matrix Library [Matrix nversion Problem][Gauss-Jordan]

[KKR_WE_RULE]

I've made a library for Delphi to work with matrices.
Everything in the library works perfectly as intended save the Matrix Inversion Procedure.

I've used the Gauss-Jordan method of Matrix Inversion & everything works excellent for square matrices upto 7x7.

From 8x8 onwards, I get a 'Invalid Floating Point' Error in this line :

tCopy.Matrix[i,j] := tCopy.Matrix[i,j] / tmp;

in the
try .. except block.

I've tried to debug this, but I cant seem to find the reason of it.

If you guys need to the entire library, I'll provide it.

I hope to get some help on this issue.

Best Regards
KKR

type
  TMatrix = record
     Row : Integer;
     Col : Integer;
     Matrix : Array OF Array OF Extended;
  End;

Procedure MatrixInverse(TargetMatrix : TMatrix; var Inverse : TMatrix; var Singular : Boolean);
var
IdentityMatrix, tCopy : TMatrix;
i, j, k : integer;
tmp, det : Extended;
Begin
  MatrixDeterminant(TargetMatrix, det);
  if (det = 0) then
  Begin
    Singular := True;
    CreateNullMatrix(TargetMatrix.Row, TargetMatrix.Col, IdentityMatrix);
    CopyMatrix(IdentityMatrix, Inverse);
    exit;
  End;
  CopyMatrix(TargetMatrix, tCopy);
  if (TargetMatrix.Row <> TargetMatrix.Col) then
  Begin
    Exit;
  End;
  InitMatrix(TargetMatrix.Row, TargetMatrix.Col, Inverse);
  CreateIdentityMatrix(TargetMatrix.Row, TargetMatrix.Col, IdentityMatrix);
  for i := 0 to TargetMatrix.Row - 1 do
   Begin
     tmp := tCopy.Matrix[i,i];
     for j := 0 to TargetMatrix.Row - 1 do
      Begin
      try
         tCopy.Matrix[i,j] := tCopy.Matrix[i,j] / tmp;
         IdentityMatrix.Matrix[i,j] := (IdentityMatrix.Matrix[i,j] / tmp);
       except
         on E:Exception do
         Begin
           Singular := True;
           CreateNullMatrix(TargetMatrix.Row, TargetMatrix.Col, IdentityMatrix);
           CopyMatrix(IdentityMatrix, Inverse);
           exit;
         End;
       End;
      End;
     for k := 0 to TargetMatrix.Row - 1 do
      Begin
        tmp := tCopy.Matrix[k, i];
        for j := 0 to TargetMatrix.Row - 1 do
         Begin
           if (i = k) then
           Begin
             break;
           End
           else
           Begin
             tCopy.Matrix[k,j] := (tCopy.Matrix[k, j] - (tCopy.matrix[i, j] * tmp));
             IdentityMatrix.Matrix[k, j] := (IdentityMatrix.Matrix[k, j] - (IdentityMatrix.Matrix[i, j] * tmp));
           End;
         End;
      End;
   End;
   DestroyMatrix(tCopy);
   CopyMatrix(IdentityMatrix, Inverse);
   DestroyMatrix(IdentityMatrix);
   Singular := False;
End;

[Mikav6]

>From 8x8 onwards, I get a 'Invalid Floating Point' Error in this line :
>

tCopy.Matrix[i,j] := tCopy.Matrix[i,j] / tmp;

Possible reason could be that tmp is zero or near it, so that result of division by tmp gives out a numer over 1.18 x 10^4932( maximum extended range).

You could change the exception handling to log more info (This works with Delphi. You might need to use different functions more suitable for your dev.envinronment):

except
on E:Exception do
Begin
  Showmessage('Error : '+e.message +'.Tmp was : '+formatFloat('',tmp)+'.Matrix['+inttostr(i)+','+inttostr(j)+']:'+formatFloat('',Matrix[i,j]));

[KKR_WE_RULE]

I use Delphi alright :)
Yes, you are correct about that 'tmp = 0' thing.
But if I add a condition like : if (tmp <> 0) then ... the inverse results for matrices is wrong & on most occasion it hangs!

[Mikav6]

It's been 15 years since I've used matrix calculus but a short read on
http://mathworld.wolfram.com/Gauss-JordanElimination.html seems to indicate
that inversion of n dimensional matrix might not be so straightforward.