Delphi - solving equations

Theory

equations of the type N equations & N unknown can be solved 
with the gaussian reduction algorithm. 

a11*x1+a12*x2+......+a1n*xn=c1
a21*x1+a22*x2+......+a2n*xn=c2
.
.
.
an1*x1+an2*x2+......+ann*xn=cn

Fill them into a matrix as :

|a11 a12 ..... a1n|   |x1|   |c1|
|a21 a22 ..... a2n|   |x2|   |c2|
|.                | * |. | = |. |
|.                |   |. |   |. |
|.                |   |. |   |. |
|an1 an2 ..... ann|   |xn|   |cn|

where the 
 aij are the constant coefficients, 
 xi the unknown, 
 ci again constants

Whole rows, including aij and ci are interchangeable, leave the x.
The row with highest absolute value of a?1 is taken as top row.
by dividing the row with the leading a, the first aij becomes 1.
By adding weighted rows, the other rows get a leading zero.

Now, from the second row down, take the one with the highest absolute
value and take it as second row. repeat the above steps

Until the matrix looks like :

|1 * ..... *|   |x1|   |*|
|0 1 ..... *|   |x2|   |*|
|0 0 1      | * |. | = |*|
|0 0        |   |. |   |*|
|0 0        |   |. |   |*|
|0 0 0 ... 1|   |xn|   |*|

Now as the solution of Xn is known, the lines can be added
from the bottom up in the same weighted manner until the matrix
looks like :

|1 0 0 0 . 0|   |x1|   |*|
|0 1 0 0 . 0|   |x2|   |*|
|0 0 1 0 . 0| * |. | = |*|
|0 0 0 1 0 0|   |. |   |*|
|0 0 0 0 1 0|   |. |   |*|
|0 0 0 0 0 1|   |xn|   |*|


The constant terms * are the solutions. The top * is x1 and so on.

Now the program :

Implementation

 class TLinEq :
 solves linear equations :

 |m11[0,0]   m12[0,1]   m13[0,2]  m14[0,3]   .. m1N[0,N-1]   | | x1 |    |c1[0]  |
 |m21[1,0]   m22[1,1]   m23[1,2]  m24[1,3]   .. m2N[1,N-1]   | | x2 |    |c2[1]  |
 |m31[2,0]   m32[2,1]   m33[2,2]  m34[2,3]   .. m3N[2,N-1]   | | x3 |    |c3[2]  |
 |m41[3,0]   m42[3,1]   m43[3,2]  m44[3,3]   .. m4N[3,N-1]   | | x4 |  = |c4[3]  |
 |..         ..         ..        ..         .. ..           | | .. |    | ..    |
 |mN1[N-1,0] mN2[N-1,1] mN3[0,2]  mN4[N-1,3] .. mNN[N-1,N-1] | | xN |    |cN[N-1]|

 the matix is the property M, the constants are property C.
 note the indices are 0 based !!
 solve returns true, when the solution was possible.

 eg :
 var eq:TLinEq;

 eq:=TLinEq.create(2);
 eq.M[0,0]:=2; eq.M[0,1]:=1; eq.C[0]:=3;
 eq.M[1,0]:=1; eq.M[1,1]:=-1;eq.C[1]:=0;

 if eq.solve then
  begin
   result[1]:=eq.C[0];
   result[2]:=eq.C[1];
  end;


type
  TEVect=class(TObject)
  private
   fp:pointer;
   fsize:integer;
   function getitem(i:integer):extended;
   procedure setitem(i:integer;z:extended);
  public
   constructor create(size:integer);
   constructor createnull(size:integer);
   destructor destroy;     override;
   procedure clear;
   procedure mult(z:extended);  //multiply with z
   procedure sub(v:TEVect);     // this:=this-v
   function copy:TEVect;
   procedure submul(z:extended;v:TEVect);     // this:=this-z*v
  published
   property Q[index:integer]:extended read getitem write setitem; default;
  end;

  TLinEq = class(TObject)
  private
    { Private declarations }
   vec:TList;
   equs:integer;
   function getMitem(index1,index2:integer):extended;
   procedure setMitem(index1,index2:integer;z:extended);
   function getCitem(index:integer):extended;
   procedure setCitem(index:integer;z:extended);
  protected
    { Protected declarations }
   function greatestpivot(n:integer):integer;
   procedure swaplines(a,b:integer);
  public
    { Public declarations }
   constructor create(dim:integer);
   destructor destroy; override;
   function solve:boolean;
  published
    { Published declarations }
   property M[index1,index2:integer]: extended read getMitem write setMitem;
   property C[index:integer]: extended read getCitem write setCitem;
  end;
......................................................................
type MyVectException = class(Exception);
type extptr=^extended;
.......................................................................
function TEVect.getitem(i:integer):extended;
var p:extptr;
    k:integer absolute p;
begin
 if (i < 0)or(i >=fsize) then 
  raise(MyVectException.create('LinEq:index out of vector boundary'));
 p:=fp;
 k:=k+i*sizeof(extended);
 result:=p^;
end;
procedure TEVect.setitem(i:integer;z:extended);
var p:extptr;
    k:integer absolute p;
begin
 if (i < 0)or(i >=fsize) then 
  raise(MyVectException.create('LinEq:index out of vector boundary'));
 p:=fp;
 k:=k+i*sizeof(extended);
 p^:=z;
end;
constructor TEVect.create(size:integer);
begin
 inherited create;
 fsize:=size;
 Getmem(fp,fsize*sizeof(extended));
end;
constructor TEVect.createnull(size:integer);
begin
 inherited create;
 fsize:=size;
 Getmem(fp,fsize*sizeof(extended));
 FillChar(fp^,fsize*sizeof(extended),0); // clear
end;
destructor TEVect.destroy;
begin
 Freemem(fp,fsize*sizeof(extended));
 inherited destroy;
end;
procedure TEVect.clear;
begin
 FillChar(fp^,fsize*sizeof(extended),0);
end;
procedure TEVect.mult(z:extended);  //multiply with z
var i:integer;
//    p:extptr;
//    k:integer absolute p;
begin
 for i:=0 to fsize-1 do Q[i]:=Q[i]*z;
end;
procedure TEVect.sub(v:TEVect);     // this:=this-v
var i:integer;
begin
 for i:=0 to fsize-1 do Q[i]:=Q[i]-v[i];
end;
function TEVect.copy:TEVect;
begin
 result:=TEVect.create(fsize);
 move(fp^,result.fp^,fsize*sizeof(extended));
end;
procedure TEVect.submul(z:extended;v:TEVect);     // this:=this-z*v
var i:integer;
begin
 for i:=0 to fsize-1 do Q[i]:=Q[i]-z*v[i];
end;
......................................................................
function TLinEq.getMitem(index1,index2:integer):extended;
begin
 result:=TEVect(vec.items[index1])[index2];
end;
procedure TLinEq.setMitem(index1,index2:integer;z:extended);
begin
 TEVect(vec.items[index1])[index2]:=z;
end;
function TLinEq.getCitem(index:integer):extended;
begin
 result:=TEVect(vec.items[index])[equs];
end;
procedure TLinEq.setCitem(index:integer;z:extended);
begin
 TEVect(vec.items[index])[equs]:=z;
end;
function TLinEq.greatestpivot(n:integer):integer;
var i,j:integer;
    t,max:extended;
begin
 j:=0;max:=0;
 for i:=n to equs-1 do begin
   t:=abs(M[i,n]);
   if t > max then begin
     max:=t;
     j:=i;
    end;
  end;
 result:=j;
end;
procedure TLinEq.swaplines(a,b:integer);
var v:TEVect;
begin
 v:=vec.items[a];
 vec.items[a]:=vec.items[b];
 vec.items[b]:=v;
end;
constructor TLinEq.create(dim:integer);
var i:integer;
    v:TEVect;
begin
 inherited create;
 vec:=TList.create;
 for i:=1 to dim do begin
   v:=TEVect.create(dim+1);
   v.clear;
   vec.add(v);
  end;
 equs:=dim;
end;
destructor TLinEq.destroy;
var i:integer;
    v:TEVect;
begin
 for i:=0 to equs-1 do begin
  v:=vec[i];
  v.destroy;
 end;
 vec.destroy;
 inherited destroy;
end;
//
Jean Bernard found an error that arised from copy/pasting :
the current line is the one with the greatest pivot.
//
function TLinEq.solve:boolean;
var l1,l2,u:integer;
    t:extended;
begin
 l1:=0; result:=true;
 while (result)and(l1 < equs) do begin
   u:=greatestpivot(l1);
   t:=M[u,l1];
   if (abs(t) < 1e-60) then result:=false
   else begin
     if (u < > l1)then swaplines(u,l1);
     TEVect(vec.items[l1]).mult(1/t);   // [l1,l1]:=1
     for l2:=0 to equs-1 do begin
       if (l1 < > l2) then begin
         t:=M[l2,l1];
         if (abs(t) > 1e-60) then
           TEVect(vec.Items[l2]).submul(t,TEVect(vec.items[l1]));
        end;
      end;
     inc(l1);
   end;
 end;
end;

findings

a standard approach is implemented. 2x2 and 3x3 systems can be
solved faster with a simpler approach without loops.
While rather big systems can be solved, it should be noted
that special approaches exist to solve sparse systems and
for those with bands.