jacobi.txt

用于开发税务票据管理的软件

Procedure JACOBI(A:matrx2; N:integer;var D:array of real;
                                   var V:matrx2;var NROT:integer);
var
    B,Z:array[0..100] of real;  IP,IQ,I,J:integer;
    SM,G,DDD,TRESH,H,T,S,SSS,THETA,C,TAU:real;
begin
    For IP:=1 To N do
    begin
        For IQ:=1 To N do
            V[IP, IQ]:=0; 
        V[IP, IP]:=1; 
    end; 
    For IP:=1 To N do
    begin
        B[IP]:=A[IP, IP];
        D[IP]:=B[IP];
        Z[IP]:=0; 
    end; 
    NROT:=0;
    For I:=1 To 50 do 
    begin
        SM:=0;
        For IP:=1 To N - 1 do
            For IQ:=IP + 1 To N do
                SM:=SM + Abs(A[IP, IQ]);
        If SM = 0  Then Exit; 
        If I < 4 Then
            TRESH:=0.2 * SM / Sqr(N)
        Else
            TRESH:=0; 
        For IP:=1 To N - 1 do
        begin
            For IQ:=IP + 1 To N do
            begin
                G:=100 * Abs(A[IP, IQ]);
                SSS:=Abs(D[IP]) + G;
                DDD:=Abs(D[IQ]) + G;
                If (I > 4) And (SSS=Abs(D[IP])) And (DDD=Abs(D[IQ])) Then
                    A[IP, IQ]:=0 
                Else If Abs(A[IP, IQ]) > TRESH Then
                begin
                    H:=D[IQ] - D[IP];
                    If Abs(H) + G = Abs(H) Then
                        T:=A[IP, IQ] / H
                    Else
                    begin
                        THETA:=0.5 * H / A[IP, IQ];
                        T:=1 / (Abs(THETA) + Sqrt(1 + Sqr(THETA)));
                        If THETA < 0  Then T:=-T;
                    end;
                    C:=1 / Sqrt(1 + Sqr(T));
                    S:=T * C;
                    TAU:=S /(1 + C);
                    H:=T * A[IP, IQ];
                    Z[IP]:=Z[IP] - H;
                    Z[IQ]:=Z[IQ] + H;
                    D[IP]:=D[IP] - H;
                    D[IQ]:=D[IQ] + H;
                    A[IP, IQ]:=0; 
                    For J:=1 To IP - 1 do
                    begin
                        G:=A[J, IP];
                        H:=A[J, IQ];
                        A[J, IP]:=G - S * (H + G * TAU);
                        A[J, IQ]:=H + S * (G - H * TAU);
                    end; 
                    For J:=IP + 1 To IQ - 1 do
                    begin
                        G:=A[IP, J];
                        H:=A[J, IQ];
                        A[IP, J]:=G - S * (H + G * TAU);
                        A[J, IQ]:=H + S * (G - H * TAU);
                    end; 
                    For J:=IQ + 1 To N do
                    begin
                        G:=A[IP, J];
                        H:=A[IQ, J];
                        A[IP, J]:=G - S * (H + G * TAU);
                        A[IQ, J]:=H + S * (G - H * TAU);
                    end; 
                    For J:=1 To N do
                    begin
                        G:=V[J, IP];
                        H:=V[J, IQ];
                        V[J, IP]:=G - S * (H + G * TAU);
                        V[J, IQ]:=H + S * (G - H * TAU);
                    end; 
                    NROT:=NROT + 1;
                end; 
            end; 
        end; 
        For IP:=1 To N do
        begin
            B[IP]:=B[IP] + Z[IP];
            D[IP]:=B[IP];
            Z[IP]:=0;
        end;
    end; 
    ShowMessage(' 50 iterations should never happen');
end;