线性方程组的共轭梯度法
procedure SPARSE(B:array of real; N:integer;var X:array of real; RSQ:real);
label 1;
const
NMAX = 500; EPS = 0.000001;
var
G,H,XI,XJ:array[0..500] of real;
J,ITER,IRST:integer;
EPS2,RP,BSQ,ANUM,ADEN,GG,DGG,GAM:real;
begin
EPS2:=N * Sqr(EPS);
IRST:=0;
1: IRST:=IRST + 1;
ASUB(X, XI);
RP:=0;
BSQ:=0;
For J:=1 To N do
begin
BSQ:=BSQ + Sqr(B[J]);
XI[J]:=XI[J] - B[J];
RP:=RP + Sqr(XI[J]);
end;
ATSUB(XI, G);
For J:=1 To N do
begin
G[J]:=-G[J];
H[J]:=G[J];
end;
For ITER:=1 To 10 * N do
begin
ASUB(H, XI);
ANUM:=0;
ADEN:=0;
For J:=1 To N do
begin
ANUM:=ANUM + G[J] * H[J];
ADEN:=ADEN + Sqr(XI[J]);
end;
If ADEN = 0 Then ShowMessage('very singular matrix');
ANUM:=ANUM / ADEN;
For J:=1 To N do
begin
XI[J]:=X[J];
X[J]:=X[J] + ANUM * H[J];
end;
ASUB(X, XJ);
RSQ:=0;
For J:=1 To N do
begin
XJ[J]:=XJ[J] - B[J];
RSQ:=RSQ + Sqr(XJ[J]);
end;
If (RSQ = RP) Or (RSQ <= BSQ * EPS2) Then Exit;
If RSQ > RP Then
begin
For J:=1 To N do
X[J]:=XI[J];
If IRST >= 3 Then Exit;
GoTo 1
end;
RP:=RSQ;
ATSUB(XJ, XI);
GG:=0;
DGG:=0;
For J:=1 To N do
begin
GG:=GG + Sqr(G[J]);
DGG:=DGG + (XI[J] + G[J]) * XI[J];
end;
If GG = 0 Then Exit;
GAM:=DGG / GG;
For J:=1 To N do
begin
G[J]:=-XI[J];
H[J]:=G[J] + GAM * H[J];
end;
end;
ShowMessage('too many iterations');
end;