trlinsolve.pas

maths lib with source

unit trlinsolve;
interface
uses Math, Ap, Sysutils;

procedure RMatrixTRSafeSolve(const A : TReal2DArray;
     N : Integer;
     var X : TReal1DArray;
     var S : Double;
     IsUpper : Boolean;
     IsTrans : Boolean;
     IsUnit : Boolean);
procedure SafeSolveTriangular(const A : TReal2DArray;
     N : Integer;
     var X : TReal1DArray;
     var S : Double;
     IsUpper : Boolean;
     IsTrans : Boolean;
     IsUnit : Boolean;
     NORMIN : Boolean;
     var CNORM : TReal1DArray);

implementation

(*************************************************************************
Utility subroutine performing the "safe" solution of system of linear
equations with triangular coefficient matrices.

The subroutine uses scaling and solves the scaled system A*x=s*b (where  s
is  a  scalar  value)  instead  of  A*x=b,  choosing  s  so  that x can be
represented by a floating-point number. The closer the system  gets  to  a
singular, the less s is. If the system is singular, s=0 and x contains the
non-trivial solution of equation A*x=0.

The feature of an algorithm is that it could not cause an  overflow  or  a
division by zero regardless of the matrix used as the input.

The algorithm can solve systems of equations with  upper/lower  triangular
matrices,  with/without unit diagonal, and systems of type A*x=b or A'*x=b
(where A' is a transposed matrix A).

Input parameters:
    A       -   system matrix. Array whose indexes range within [0..N-1, 0..N-1].
    N       -   size of matrix A.
    X       -   right-hand member of a system.
                Array whose index ranges within [0..N-1].
    IsUpper -   matrix type. If it is True, the system matrix is the upper
                triangular and is located in  the  corresponding  part  of
                matrix A.
    Trans   -   problem type. If it is True, the problem to be  solved  is
                A'*x=b, otherwise it is A*x=b.
    IsUnit  -   matrix type. If it is True, the system matrix has  a  unit
                diagonal (the elements on the main diagonal are  not  used
                in the calculation process), otherwise the matrix is considered
                to be a general triangular matrix.

Output parameters:
    X       -   solution. Array whose index ranges within [0..N-1].
    S       -   scaling factor.

  -- LAPACK auxiliary routine (version 3.0) --
     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
     Courant Institute, Argonne National Lab, and Rice University
     June 30, 1992
*************************************************************************)
procedure RMatrixTRSafeSolve(const A : TReal2DArray;
     N : Integer;
     var X : TReal1DArray;
     var S : Double;
     IsUpper : Boolean;
     IsTrans : Boolean;
     IsUnit : Boolean);
var
    NORMIN : Boolean;
    CNORM : TReal1DArray;
    A1 : TReal2DArray;
    X1 : TReal1DArray;
    I : Integer;
begin
    
    //
    // From 0-based to 1-based
    //
    NORMIN := False;
    SetLength(A1, N+1, N+1);
    SetLength(X1, N+1);
    I:=1;
    while I<=N do
    begin
        APVMove(@A1[I][0], 1, N, @A[I-1][0], 0, N-1);
        Inc(I);
    end;
    APVMove(@X1[0], 1, N, @X[0], 0, N-1);
    
    //
    // Solve 1-based
    //
    SafeSolveTriangular(A1, N, X1, S, IsUpper, IsTrans, IsUnit, NORMIN, CNORM);
    
    //
    // From 1-based to 0-based
    //
    APVMove(@X[0], 0, N-1, @X1[0], 1, N);
end;


(*************************************************************************
Obsolete 1-based subroutine.
See RMatrixTRSafeSolve for 0-based replacement.
*************************************************************************)
procedure SafeSolveTriangular(const A : TReal2DArray;
     N : Integer;
     var X : TReal1DArray;
     var S : Double;
     IsUpper : Boolean;
     IsTrans : Boolean;
     IsUnit : Boolean;
     NORMIN : Boolean;
     var CNORM : TReal1DArray);
var
    I : Integer;
    IMAX : Integer;
    J : Integer;
    JFIRST : Integer;
    JINC : Integer;
    JLAST : Integer;
    JM1 : Integer;
    JP1 : Integer;
    IP1 : Integer;
    IM1 : Integer;
    K : Integer;
    Flg : Integer;
    V : Double;
    VD : Double;
    BIGNUM : Double;
    GROW : Double;
    REC : Double;
    SMLNUM : Double;
    SUMJ : Double;
    TJJ : Double;
    TJJS : Double;
    TMAX : Double;
    TSCAL : Double;
    USCAL : Double;
    XBND : Double;
    XJ : Double;
    XMAX : Double;
    NOTRAN : Boolean;
    UPPER : Boolean;
    NOUNIT : Boolean;
    i_ : Integer;
begin
    UPPER := IsUpper;
    NOTRAN :=  not IsTrans;
    NOUNIT :=  not IsUnit;
    
    //
    // Quick return if possible
    //
    if N=0 then
    begin
        Exit;
    end;
    
    //
    // Determine machine dependent parameters to control overflow.
    //
    SMLNUM := MinRealNumber/(MachineEpsilon*2);
    BIGNUM := 1/SMLNUM;
    S := 1;
    if  not NORMIN then
    begin
        SetLength(CNORM, N+1);
        
        //
        // Compute the 1-norm of each column, not including the diagonal.
        //
        if UPPER then
        begin
            
            //
            // A is upper triangular.
            //
            J:=1;
            while J<=N do
            begin
                V := 0;
                K:=1;
                while K<=J-1 do
                begin
                    V := V+AbsReal(A[K,J]);
                    Inc(K);
                end;
                CNORM[J] := V;
                Inc(J);
            end;
        end
        else
        begin
            
            //
            // A is lower triangular.
            //
            J:=1;
            while J<=N-1 do
            begin
                V := 0;
                K:=J+1;
                while K<=N do
                begin
                    V := V+AbsReal(A[K,J]);
                    Inc(K);
                end;
                CNORM[J] := V;
                Inc(J);
            end;
            CNORM[N] := 0;
        end;
    end;
    
    //
    // Scale the column norms by TSCAL if the maximum element in CNORM is
    // greater than BIGNUM.
    //
    IMAX := 1;
    K:=2;
    while K<=N do
    begin
        if CNORM[K]>CNORM[IMAX] then
        begin
            IMAX := K;
        end;
        Inc(K);
    end;
    TMAX := CNORM[IMAX];
    if TMAX<=BIGNUM then
    begin
        TSCAL := 1;
    end
    else
    begin
        TSCAL := 1/(SMLNUM*TMAX);
        APVMul(@CNORM[0], 1, N, TSCAL);
    end;
    
    //
    // Compute a bound on the computed solution vector to see if the
    // Level 2 BLAS routine DTRSV can be used.
    //
    J := 1;
    K:=2;
    while K<=N do
    begin
        if AbsReal(X[K])>AbsReal(X[J]) then
        begin
            J := K;
        end;
        Inc(K);
    end;
    XMAX := ABSReal(X[J]);
    XBND := XMAX;
    if NOTRAN then
    begin
        
        //
        // Compute the growth in A * x = b.
        //
        if UPPER then
        begin
            JFIRST := N;
            JLAST := 1;
            JINC := -1;
        end
        else
        begin
            JFIRST := 1;
            JLAST := N;
            JINC := 1;
        end;
        if TSCAL<>1 then
        begin
            GROW := 0;
        end
        else
        begin
            if NOUNIT then
            begin
                
                //
                // A is non-unit triangular.
                //
                // Compute GROW = 1/G(j) and XBND = 1/M(j).
                // Initially, G(0) = max{x(i), i=1,...,n}.
                //
                GROW := 1/Max(XBND, SMLNUM);
                XBND := GROW;
                J := JFIRST;
                while (JINC>0) and (J<=JLAST) or (JINC<0) and (J>=JLAST) do
                begin
                    
                    //
                    // Exit the loop if the growth factor is too small.
                    //
                    if GROW<=SMLNUM then
                    begin
                        Break;
                    end;
                    
                    //
                    // M(j) = G(j-1) / abs(A(j,j))
                    //
                    TJJ := ABSReal(A[J,J]);
                    XBND := Min(XBND, Min(1, TJJ)*GROW);
                    if TJJ+CNORM[J]>=SMLNUM then
                    begin
                        
                        //
                        // G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )
                        //
                        GROW := GROW*(TJJ/(TJJ+CNORM[J]));
                    end
                    else
                    begin
                        
                        //
                        // G(j) could overflow, set GROW to 0.
                        //
                        GROW := 0;
                    end;
                    if J=JLAST then
                    begin
                        GROW := XBND;
                    end;
                    J := J+JINC;