rsolve.pas

maths lib with source

unit rsolve;
interface
uses Math, Ap, Sysutils, lu;

function RMatrixLUSolve(const A : TReal2DArray;
     const Pivots : TInteger1DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
function RMatrixSolve(A : TReal2DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
function SolveSystemLU(const A : TReal2DArray;
     const Pivots : TInteger1DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
function SolveSystem(A : TReal2DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;

implementation

(*************************************************************************
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 RMatrixLUSolve(const A : TReal2DArray;
     const Pivots : TInteger1DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
var
    Y : TReal1DArray;
    I : Integer;
    J : Integer;
    V : Double;
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 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 := APVDotProduct(@A[I][0], 0, I-1, @Y[0], 0, I-1);
        Y[I] := B[I]-V;
        Inc(I);
    end;
    
    //
    // Ux = y
    //
    X[N-1] := Y[N-1]/A[N-1,N-1];
    I:=N-2;
    while I>=0 do
    begin
        V := APVDotProduct(@A[I][0], I+1, N-1, @X[0], I+1, N-1);
        X[I] := (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 RMatrixSolve(A : TReal2DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
var
    Pivots : TInteger1DArray;
    I : Integer;
begin
    A := DynamicArrayCopy(A);
    B := DynamicArrayCopy(B);
    RMatrixLU(A, N, N, Pivots);
    Result := RMatrixLUSolve(A, Pivots, B, N, X);
end;


(*************************************************************************
Obsolete 1-based subroutine
*************************************************************************)
function SolveSystemLU(const A : TReal2DArray;
     const Pivots : TInteger1DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
var
    Y : TReal1DArray;
    I : Integer;
    J : Integer;
    V : Double;
    IP1 : Integer;
    IM1 : Integer;
begin
    B := DynamicArrayCopy(B);
    SetLength(Y, N+1);
    SetLength(X, N+1);
    Result := True;
    I:=1;
    while I<=N do
    begin
        if 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 := APVDotProduct(@A[I][0], 1, IM1, @Y[0], 1, IM1);
        Y[I] := B[I]-V;
        Inc(I);
    end;
    
    //
    // Ux = y
    //
    X[N] := Y[N]/A[N,N];
    I:=N-1;
    while I>=1 do
    begin
        IP1 := I+1;
        V := APVDotProduct(@A[I][0], IP1, N, @X[0], IP1, N);
        X[I] := (Y[I]-V)/A[I,I];
        Dec(I);
    end;
end;


(*************************************************************************
Obsolete 1-based subroutine
*************************************************************************)
function SolveSystem(A : TReal2DArray;
     B : TReal1DArray;
     N : Integer;
     var X : TReal1DArray):Boolean;
var
    Pivots : TInteger1DArray;
    I : Integer;
begin
    A := DynamicArrayCopy(A);
    B := DynamicArrayCopy(B);
    LUDecomposition(A, N, N, Pivots);
    Result := SolveSystemLU(A, Pivots, B, N, X);
end;


end.