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;