csolve.pas

maths lib with source

unit csolve;
interface
uses Math, Ap, Sysutils, clu;

function CMatrixLUSolve(const A : TComplex2DArray;
     const Pivots : TInteger1DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
function CMatrixSolve(A : TComplex2DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
function ComplexSolveSystemLU(const A : TComplex2DArray;
     const Pivots : TInteger1DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
function ComplexSolveSystem(A : TComplex2DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;

implementation

procedure TestComplexSolveSystemLU();forward;


(*************************************************************************
Solving a system of linear equations with a system matrix given by its
LU decomposition.

The algorithm solves a system of linear equations whose matrix is given by
its LU decomposition. In case of a singular matrix, the algorithm  returns
False.

The algorithm solves systems with a square matrix only.

Input parameters:
    A       -   LU decomposition of a system matrix in compact  form  (the
                result of the RMatrixLU subroutine).
    Pivots  -   row permutation table (the result of a
                RMatrixLU subroutine).
    B       -   right side of a system.
                Array whose index ranges within [0..N-1].
    N       -   size of matrix A.

Output parameters:
    X       -   solution of a system.
                Array whose index ranges within [0..N-1].

Result:
    True, if the matrix is not singular.
    False, if the matrux is singular. In this case, X doesn't contain a
solution.

  -- ALGLIB --
     Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************)
function CMatrixLUSolve(const A : TComplex2DArray;
     const Pivots : TInteger1DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
var
    Y : TComplex1DArray;
    I : Integer;
    J : Integer;
    V : Complex;
    i_ : Integer;
begin
    B := DynamicArrayCopy(B);
    SetLength(Y, N-1+1);
    SetLength(X, N-1+1);
    Result := True;
    I:=0;
    while I<=N-1 do
    begin
        if C_EqualR(A[I,I],0) then
        begin
            Result := False;
            Exit;
        end;
        Inc(I);
    end;
    
    //
    // pivots
    //
    I:=0;
    while I<=N-1 do
    begin
        if Pivots[I]<>I then
        begin
            V := B[I];
            B[I] := B[Pivots[I]];
            B[Pivots[I]] := V;
        end;
        Inc(I);
    end;
    
    //
    // Ly = b
    //
    Y[0] := B[0];
    I:=1;
    while I<=N-1 do
    begin
        V := C_Complex(0.0);
        for i_ := 0 to I-1 do
        begin
            V := C_Add(V,C_Mul(A[I,i_],Y[i_]));
        end;
        Y[I] := C_Sub(B[I],V);
        Inc(I);
    end;
    
    //
    // Ux = y
    //
    X[N-1] := C_Div(Y[N-1],A[N-1,N-1]);
    I:=N-2;
    while I>=0 do
    begin
        V := C_Complex(0.0);
        for i_ := I+1 to N-1 do
        begin
            V := C_Add(V,C_Mul(A[I,i_],X[i_]));
        end;
        X[I] := C_Div(C_Sub(Y[I],V),A[I,I]);
        Dec(I);
    end;
end;


(*************************************************************************
Solving a system of linear equations.

The algorithm solves a system of linear equations by using the
LU decomposition. The algorithm solves systems with a square matrix only.

Input parameters:
    A   -   system matrix.
            Array whose indexes range within [0..N-1, 0..N-1].
    B   -   right side of a system.
            Array whose indexes range within [0..N-1].
    N   -   size of matrix A.

Output parameters:
    X   -   solution of a system.
            Array whose index ranges within [0..N-1].

Result:
    True, if the matrix is not singular.
    False, if the matrix is singular. In this case, X doesn't contain a
solution.

  -- ALGLIB --
     Copyright 2005-2008 by Bochkanov Sergey
*************************************************************************)
function CMatrixSolve(A : TComplex2DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
var
    Pivots : TInteger1DArray;
    I : Integer;
begin
    A := DynamicArrayCopy(A);
    B := DynamicArrayCopy(B);
    CMatrixLU(A, N, N, Pivots);
    Result := CMatrixLUSolve(A, Pivots, B, N, X);
end;


(*************************************************************************
Obsolete 1-based subroutine
*************************************************************************)
function ComplexSolveSystemLU(const A : TComplex2DArray;
     const Pivots : TInteger1DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
var
    Y : TComplex1DArray;
    I : Integer;
    V : Complex;
    IP1 : Integer;
    IM1 : Integer;
    i_ : Integer;
begin
    B := DynamicArrayCopy(B);
    SetLength(Y, N+1);
    SetLength(X, N+1);
    Result := True;
    I:=1;
    while I<=N do
    begin
        if C_EqualR(A[I,I],0) then
        begin
            Result := False;
            Exit;
        end;
        Inc(I);
    end;
    
    //
    // pivots
    //
    I:=1;
    while I<=N do
    begin
        if Pivots[I]<>I then
        begin
            V := B[I];
            B[I] := B[Pivots[I]];
            B[Pivots[I]] := V;
        end;
        Inc(I);
    end;
    
    //
    // Ly = b
    //
    Y[1] := B[1];
    I:=2;
    while I<=N do
    begin
        IM1 := I-1;
        V := C_Complex(0.0);
        for i_ := 1 to IM1 do
        begin
            V := C_Add(V,C_Mul(A[I,i_],Y[i_]));
        end;
        Y[I] := C_Sub(B[I],V);
        Inc(I);
    end;
    
    //
    // Ux = y
    //
    X[N] := C_Div(Y[N],A[N,N]);
    I:=N-1;
    while I>=1 do
    begin
        IP1 := I+1;
        V := C_Complex(0.0);
        for i_ := IP1 to N do
        begin
            V := C_Add(V,C_Mul(A[I,i_],X[i_]));
        end;
        X[I] := C_Div(C_Sub(Y[I],V),A[I,I]);
        Dec(I);
    end;
end;


(*************************************************************************
Obsolete 1-based subroutine
*************************************************************************)
function ComplexSolveSystem(A : TComplex2DArray;
     B : TComplex1DArray;
     N : Integer;
     var X : TComplex1DArray):Boolean;
var
    Pivots : TInteger1DArray;
begin
    A := DynamicArrayCopy(A);
    B := DynamicArrayCopy(B);
    ComplexLUDecomposition(A, N, N, Pivots);
    Result := ComplexSolveSystemLU(A, Pivots, B, N, X);
end;


procedure TestComplexSolveSystemLU();
var
    I : Integer;
    J : Integer;
    Err : Double;
    V : Complex;
    A : TComplex2DArray;
    TX : TComplex1DArray;
    X : TComplex1DArray;
    B : TComplex1DArray;
    N : Integer;
    Pass : Integer;
    PassCount : Integer;
    i_ : Integer;
begin
    Err := 0;
    PassCount := 1000;
    Pass:=1;
    while Pass<=PassCount do
    begin
        N := 1+RandomInteger(20);
        SetLength(A, N+1, N+1);
        SetLength(TX, N+1);
        SetLength(B, N+1);
        
        //
        // init A, TX
        //
        I:=1;
        while I<=N do
        begin
            J:=1;
            while J<=N do
            begin
                A[I,J].X := 2*RandomReal-1;
                A[I,J].Y := 2*RandomReal-1;
                Inc(J);
            end;
            Inc(I);
        end;
        A[1+RandomInteger(N),1+RandomInteger(N)] := C_Complex(10);
        A[1+RandomInteger(N),1+RandomInteger(N)] := C_Complex(7);
        I:=1;
        while I<=N do
        begin
            TX[I].X := 2*RandomReal-1;
            TX[I].Y := 2*RandomReal-1;
            Inc(I);
        end;
        I:=1;
        while I<=N do
        begin
            V := C_Complex(0.0);
            for i_ := 1 to N do
            begin
                V := C_Add(V,C_Mul(A[I,i_],TX[i_]));
            end;
            B[I] := V;
            Inc(I);
        end;
        
        //
        // solve
        //
        ComplexSolveSystem(A, B, N, X);
        
        //
        // test
        //
        I:=1;
        while I<=N do
        begin
            Err := Max(Err, AbsComplex(C_Sub(TX[I],X[I])));
            Inc(I);
        end;
        Inc(Pass);
    end;
    Write(Format('TESTING COMPLEX SOLVE SYSTEM'#13#10'',[]));
    Write(Format('Pass count is %0d'#13#10'SolveSystem absolute error is %5.4e'#13#10'',[
        PassCount,
        Err]));
end;


end.