testspdrcondunit.pas

maths lib with source

unit testspdrcondunit;
interface
uses Math, Ap, Sysutils, trlinsolve, cholesky, estnorm, spdrcond;

function TestRCondCholesky(Silent : Boolean):Boolean;
function testspdrcondunit_test_silent():Boolean;
function testspdrcondunit_test():Boolean;

implementation

procedure MakeACopy(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);forward;
function InvMatTR(var A : TReal2DArray;
     N : Integer;
     IsUpper : Boolean;
     IsUnitTriangular : Boolean):Boolean;forward;
function InvMatLU(var A : TReal2DArray;
     const Pivots : TInteger1DArray;
     N : Integer):Boolean;forward;
procedure MatLU(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     var Pivots : TInteger1DArray);forward;
function InvMat(var A : TReal2DArray; N : Integer):Boolean;forward;
procedure RefRCond(const A : TReal2DArray;
     N : Integer;
     var RC1 : Double;
     var RCInf : Double);forward;


function TestRCondCholesky(Silent : Boolean):Boolean;
var
    L : TReal2DArray;
    A : TReal2DArray;
    A2 : TReal2DArray;
    MinIJ : Integer;
    P : TInteger1DArray;
    N : Integer;
    MaxN : Integer;
    I : Integer;
    J : Integer;
    Pass : Integer;
    PassCount : Integer;
    WasErrors : Boolean;
    Err50 : Boolean;
    Err90 : Boolean;
    ErrLess : Boolean;
    ERC1 : Double;
    ERCInf : Double;
    Q50 : TReal1DArray;
    Q90 : TReal1DArray;
    V : Double;
    Threshold50 : Double;
    Threshold90 : Double;
    HTask : Integer;
    UpperIn : Boolean;
    i_ : Integer;
begin
    WasErrors := False;
    Err50 := False;
    Err90 := False;
    ErrLess := False;
    MaxN := 10;
    PassCount := 100;
    Threshold50 := 0.5;
    Threshold90 := 0.1;
    SetLength(Q50, 1+1);
    SetLength(Q90, 1+1);
    
    //
    // process
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(L, N-1+1, N-1+1);
        SetLength(A, N-1+1, N-1+1);
        SetLength(A2, N-1+1, N-1+1);
        I:=0;
        while I<=1 do
        begin
            Q50[I] := 0;
            Q90[I] := 0;
            Inc(I);
        end;
        HTask:=0;
        while HTask<=1 do
        begin
            Pass:=1;
            while Pass<=PassCount do
            begin
                UpperIn := HTask=0;
                
                //
                // Prepare task:
                // * A contains full matrix A
                // * L contains its Cholesky factor (upper or lower)
                // * A2 contains upper/lower triangle of A
                //
                I:=0;
                while I<=N-1 do
                begin
                    J:=I+1;
                    while J<=N-1 do
                    begin
                        L[I,J] := 2*RandomReal-1;
                        L[J,I] := L[I,J];
                        Inc(J);
                    end;
                    L[I,I] := 1.1+RandomReal;
                    Inc(I);
                end;
                I:=0;
                while I<=N-1 do
                begin
                    J:=0;
                    while J<=N-1 do
                    begin
                        MinIJ := Min(I, J);
                        V := 0.0;
                        for i_ := 0 to MinIJ do
                        begin
                            V := V + L[I,i_]*L[i_,J];
                        end;
                        A[I,J] := V;
                        A2[I,J] := V;
                        Inc(J);
                    end;
                    Inc(I);
                end;
                I:=0;
                while I<=N-1 do
                begin
                    J:=0;
                    while J<=N-1 do
                    begin
                        if UpperIn then
                        begin
                            if J<I then
                            begin
                                L[I,J] := 0;
                                A2[I,J] := 0;
                            end;
                        end
                        else
                        begin
                            if I<J then
                            begin
                                L[I,J] := 0;
                                A2[I,J] := 0;
                            end;
                        end;
                        Inc(J);
                    end;
                    Inc(I);
                end;
                RefRCond(A, N, ERC1, ERCInf);
                
                //
                // normal
                //
                V := 1/SPDMatrixRCond(A, N, UpperIn);
                if V>=Threshold50*ERC1 then
                begin
                    Q50[0] := Q50[0]+1/PassCount;
                end;
                if V>=Threshold90*ERC1 then
                begin
                    Q90[0] := Q90[0]+1/PassCount;
                end;
                ErrLess := ErrLess or (V>ERC1*1.001);
                
                //
                // Cholesky
                //
                V := 1/SPDMatrixCholeskyRCond(L, N, UpperIn);
                if V>=Threshold50*ERC1 then
                begin
                    Q50[1] := Q50[1]+1/PassCount;
                end;
                if V>=Threshold90*ERC1 then
                begin
                    Q90[1] := Q90[1]+1/PassCount;
                end;
                ErrLess := ErrLess or (V>ERC1*1.001);
                Inc(Pass);
            end;
            Inc(HTask);
        end;
        I:=0;
        while I<=1 do
        begin
            Err50 := Err50 or (Q50[I]<0.50);
            Err90 := Err90 or (Q90[I]<0.90);
            Inc(I);
        end;
        Inc(N);
    end;
    
    //
    // report
    //
    WasErrors := Err50 or Err90 or ErrLess;
    if  not Silent then
    begin
        Write(Format('TESTING RCOND (SYMMETRIC)'#13#10'',[]));
        Write(Format('50%% within [0.5*cond,cond]:              ',[]));
        if Err50 or ErrLess then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('90%% within [0.1*cond,cond]               ',[]));
        if Err90 or ErrLess then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        if WasErrors then
        begin
            Write(Format('TEST FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('TEST PASSED'#13#10'',[]));
        end;
        Write(Format(''#13#10''#13#10'',[]));
    end;
    Result :=  not WasErrors;
end;


(*************************************************************************
Copy
*************************************************************************)
procedure MakeACopy(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);
var
    I : Integer;
    J : Integer;
begin
    SetLength(B, M-1+1, N-1+1);
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            B[I,J] := A[I,J];
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
triangular inverse
*************************************************************************)
function InvMatTR(var A : TReal2DArray;
     N : Integer;
     IsUpper : Boolean;
     IsUnitTriangular : Boolean):Boolean;
var
    NOUNIT : Boolean;
    I : Integer;
    J : Integer;
    V : Double;
    AJJ : Double;
    T : TReal1DArray;
    i_ : Integer;
begin
    Result := True;
    SetLength(T, N-1+1);
    
    //
    // Test the input parameters.
    //
    NOUNIT :=  not IsUnitTriangular;
    if IsUpper then
    begin
        
        //
        // Compute inverse of upper triangular matrix.
        //
        J:=0;
        while J<=N-1 do
        begin
            if NOUNIT then
            begin
                if A[J,J]=0 then
                begin
                    Result := False;
                    Exit;
                end;
                A[J,J] := 1/A[J,J];
                AJJ := -A[J,J];
            end
            else
            begin
                AJJ := -1;
            end;
            
            //
            // Compute elements 1:j-1 of j-th column.
            //
            if J>0 then
            begin
                for i_ := 0 to J-1 do
                begin
                    T[i_] := A[i_,J];
                end;
                I:=0;
                while I<=J-1 do
                begin
                    if I<J-1 then
                    begin
                        V := APVDotProduct(@A[I][0], I+1, J-1, @T[0], I+1, J-1);
                    end
                    else
                    begin
                        V := 0;
                    end;
                    if NOUNIT then
                    begin
                        A[I,J] := V+A[I,I]*T[I];
                    end
                    else
                    begin
                        A[I,J] := V+T[I];
                    end;
                    Inc(I);
                end;
                for i_ := 0 to J-1 do
                begin
                    A[i_,J] := AJJ*A[i_,J];
                end;
            end;
            Inc(J);
        end;
    end
    else
    begin
        
        //
        // Compute inverse of lower triangular matrix.
        //
        J:=N-1;
        while J>=0 do
        begin
            if NOUNIT then
            begin
                if A[J,J]=0 then