如何解线性方程组？我愿意出大洋200大伙来看看吧，谢谢！(50分)

以下是求解大型稀疏方程组得全选主元高斯-约当消去法，
（清华出得，c常用算法程序集，徐士良编写得）
的c语言原程序，本人把它翻译为delphi版本的，可是编译通过了
计算有错误，本人愿意出200大洋，各位帮忙看一下了，谢谢。
(其它分本人有其它问题等你去那，到了就给分！！）
http://211.101.4.25/delphibbs/dispq.asp?lid=654773
http://211.101.4.25/delphibbs/dispq.asp?lid=552840
如果那位愿意给出其它方法的原码也好，！
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
int aggje(a,n,b)
int n;
double a[],b[];
{ int *js,i,j,k,is,u,v;
double d,t;
js=malloc(n*sizeof(int));
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ t=fabs(a[i*n+j]);
if (t>d) {d=t; js[k]=j; is=i;}
}
if (d+1.0==1.0)
{ free(js); printf("fail/n"); return(0);}
if (is!=k)
{ for (j=k; j<=n-1; j++)
{ u=k*n+j; v=is*n+j;
t=a; a=a[v]; a[v]=t;
}
t=b[k]; b[k]=b[is]; b[is]=t;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
t=a; a=a[v]; a[v]=t;
}
t=a[k*n+k];
for (j=k+1; j<=n-1; j++)
{ u=k*n+j;
if (a!=0.0) a=a/t;
}
b[k]=b[k]/t;
for (j=k+1; j<=n-1; j++)
{ u=k*n+j;
if (a!=0.0)
{ for (i=0; i<=n-1; i++)
{ v=i*n+k;
if ((i!=k)&&(a[v]!=0.0))
{ is=i*n+j;
a[is]=a[is]-a[v]*a;
}
}
}
}
for (i=0; i<=n-1; i++)
{ u=i*n+k;
if ((i!=k)&&(a!=0.0))
b=b-a*b[k];
}
}
for (k=n-1; k>=0; k--)
if (k!=js[k])
{ t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}
free(js);
return(1);
}
delphi版本的在下面，本人水平差，编译正确，计算总是错误的，各位高手看看，谢谢了
形如AX=b,已知a,b,求x
procedure CalcFC(a:TD2Array; b:TDArray;Var x:TDArray; n:Integer);//高斯消元法解线性方程组
var
A1:array of array of double;//a1放矩阵a,
js:array of integer;
b1:array of double;//b1放矩阵b
i,j,k,isa:integer;
d,t:double;
begin
SetLength(A1,n,n);
SetLength(js,n);
Setlength(b1,n);
for i:=0 to n-1 do
begin
for j:=0 to n-1 do
A1[i,j]:=A.Value[i,j];//取出a的直
b1:=B.Value;//取出b的直，
end;
for k:=0 to n-1 do//开始和c版本一样了，各位看看
begin
d:=0.0;
for i:=k to n-1 do
for j:=k to n-1 do
begin
T:=abs(a1[i,j]);
if t>d then
begin
d:=t;
js[k]:=j;
isa:=i;
end;
end;
if d+1.0=1.0 then showmessage('ok');
if js[k]<>k then
for i:=0 to n-1 do
begin
t:=a1[i,k];
a1[i,k]:=a1[i,js[k]];
a1[i,js[k]]:=t;
end;
if isa<>k then
for j:=k to n do
begin
t:=a1[k,j];
a1[k,j]:=a1[isa,j];
a1[isa,j]:=t;
end;
t:=a1[k,k];
for j:=k+1 to n-1 do
if a1[k,j]<>0 then
A1[k,j]:=A1[k,j]/t;
b1[k]:=b1[k]/t;
for j:=k+1 to n-1 do
if a1[k,j]<>0 then
begin
for i:=0 to n-1 do
if (i<>k) and (a1[i,k]<>0) then
a1[i,j]:=a1[i,j]-a1[i,k]*a1[k,j];
end;
for i:=0 to n-1 do
if (i<>k) and (a1[i,k]<>0) then
b1:=b1-a1[i,k]*b1[k];
end;
for k:=n-1 downto 0 do
begin
if js[k]<>K then
begin
t:=b1[k];
b1[k]:=b1[js[k]];
b1[js[k]]:=t;
end;
end;
for k:=0 to n-1 do
x.value[k]:=b1[k];
end;
各位给个修改意见！谢谢

现在也有《DELPHI常用算法程序集》了，我昨天刚在海淀图书城对面的书店看到，那上面应该
也有这个算法的，你去看看吧。

高斯-约当消去法。倒！高度贫血，一看就头晕！
自己写一个得了！ 算法如下：
1, 1/a[k,k]得到新的a[k,k]
2, a[k,j]a[k,k]得到新的a[k,j] j=1,...,n j<>k
3, a[i,j]-a[i,k]a[k,j]得到a[i.j] i,j=1,...,n i,j<>k
4, -a[i,k]a[k,k]得到a[i,k] i=1,...,n i<>k

你说的方法我知道了，我就是用这个的，不过你想想，大型稀疏矩阵，好多对角线元素为0
如果你直接a[k,k]:=1/a[k,k],要出问题的拉~！！所以，想来想去，还是高斯好啊！
其实你那个也是高斯，不过简化了

我不相信没人能搞定，这可是常见问题啊！！

真的是个难题？？？？？？那有没有比这个更快得办法呢？？！！！！

多人接受答案了。