Back to home

Calculating Matrix Determinant Using Gauss-Jordan Elimination





This can be considered as the continuation from
the first one. This time, the Gauss-Jordan elimination
is used. However, the full-fledged version of
the elimination is not needed.

This elimination comes into play due to the fact
that if the elements in lower triangle of the
matrix are all zeroes, the determinant is simply
the product of the diagonals. This also applies
in the upper triangle zeroes. It means:

a ? ? ?        a 0 0 0
0 b ? ?   or   ? b 0 0  's determinant = a*b*c*d
0 0 c ?        ? ? c 0
0 0 0 d        ? ? ? d

? = don't care.

So, the elimination works to shape the matrix into
one of these forms. My solution aims for the lower
triangle zeroes. Full Gauss-Jordan elimination
requires both upper and lower triangles be zeroes.

I don't want to explain this in details as you can
easily find it in good algebra books. I just give
you a debug procedure "printarr", so that you can
inject a line at mark {1} in order to learn how the
algorithm works. If you'd like to have a more
detailed look, put it at mark {2} instead.

Try it using this matrix:
 1  0 -1  0
 0  3  0  2
-1  0  2  1
 0  5  0  7

and then this one to understand when to "re-sort":
 0  0 -1  0
 0  3  0  2
-1  0  2  1
 0  5  0  7


Here is my solution:
----------------------------------------------------
type { you can modify this dimension }
  arr = array[1..4, 1..4] of real;

procedure printarr(new_a: arr; dim: integer);
var i, j: integer;
begin
   for i := 1 to dim do
   begin
      for j := 1 to dim do
         write(new_a[i,j]:5:3,' ');
      writeln;
   end;
end;

{ Gauss-Jordan Elimination }
function gauss(a: arr; dim: integer): real;
var
  new_a : arr;
  i, j, k : integer;
  factor, temp, det : real;
begin
   { copy the array }
   for i := 1 to dim do
      for j := 1 to dim do
         new_a[i,j] := a[i,j];

   det := 1.0;
   { do the elimination }
   for i := 1 to dim-1 do
   begin
     { if the main diagonal value is zero }
     { re-sort the array }
     if (new_a[i,i] = 0) then
     begin
       for j := i+1 to dim do
       begin
         if (new_a[j,i] <> 0) then
         begin
           for k := 1 to dim do
           begin
             temp := new_a[i,k];
             new_a[i,k] := new_a[j,k];
             new_a[j,k] := temp;
           end;
           { For Gauss-Jordan Elimination, }
           { if we do a switch, the determinant }
           { switches sign. }
           det := -det;
           break;
         end;
       end;
     end;

     { if after the resorting, the value is still zero }
     { then the determinant is definitely zero }
     if (new_a[i,i] = 0) then
     begin
       gauss := 0; exit;
     end;

     { eliminate the lower rows to achieve triangular zeroes }
     for j := i+1 to dim do
     begin
       if (new_a[j,i] <> 0) then
       begin
         factor := (new_a[j,i] * 1.0) / new_a[i,i];
         for k := i to dim do
         begin
           new_a[j,k] := new_a[j,k] - factor * new_a[i,k];
           {2}
         end;
         {1}
       end;
     end;
   end;

   { calculate the main diagonal }
   for i := 1 to dim do
      det := det * new_a[i,i];

   gauss := det;
end;
----------------------------------------------------

Back to home


Roby Joehanes, © 2001