Gaussian elimination

From: Michael Trofimov (spamtrap_at_crayne.org)
Date: 01/10/05 

Date: Mon, 10 Jan 2005 21:33:39 +0000 (UTC)

Hi,

I am looking for a fast code for the Gaussian elimination (classical method
of linear algebra for solving systems of linear equations). Target platform
is Pentium4 with HT. I think to replace the following code (or its part) in
Delphi 7 for optimization and to use build-in assembler for it:

procedure Gauss (var A : TLinSMat; var B,X : TLinSVect);

// solve system AX=B, sz is the system size =2..41

// TLinSMat = array [1..41,1..41] of extended; -- may be real (8 bytes, if
it would be faster)?

// TLinSVect = array [1..41] of extended; -- may be real?

var

i,j,k : integer;

q : extended; //-- may be real?

begin

for k :=1 to sz do

begin

// if A[k,k]<>0 then -- redundant for my task

for i:=k+1 to sz do

begin

q := - A[i,k]/A[k,k];

for j:=1 to k-1 do

A[i,j] := A[i,j]+q*A[k,j];

A[i,k] := 0.0;

for j:=k+1 to sz do

A[i,j] := A[i,j]+q*A[k,j];

B[i] := B[i]+q*B[k];

end;

end;

X[sz] := B[sz]/A[sz,sz];

for i:=sz-1 downto 1 do

begin

q := 0.0;

for j:=1 to sz-i do

q := q + A[i,i+j]*X[i+j];

X[i] := (B[i]-q)/A[i,i];

end;

end;

Any ideas, references etc:

Many thanks,

Michael Trofimov.