testdetunit.pas

maths lib with source

unit testdetunit;
interface
uses Math, Ap, Sysutils, lu, det;

function TestDet(Silent : Boolean):Boolean;
function testdetunit_test_silent():Boolean;
function testdetunit_test():Boolean;

implementation

var
    DetErrors : Boolean;

procedure FillSparseA(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     Sparcity : Double);forward;
procedure TestProblem(const A : TReal2DArray; N : Integer);forward;
procedure MakeACopy(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);forward;
function DetTriangle(A : TReal2DArray; const N : Integer):Double;forward;


(*************************************************************************
Main unittest subroutine
*************************************************************************)
function TestDet(Silent : Boolean):Boolean;
var
    MaxN : Integer;
    GPassCount : Integer;
    Threshold : Double;
    A : TReal2DArray;
    N : Integer;
    GPass : Integer;
    I : Integer;
    J : Integer;
    WasErrors : Boolean;
begin
    DetErrors := False;
    WasErrors := False;
    MaxN := 8;
    GPassCount := 5;
    Threshold := 5*100*MachineEpsilon;
    SetLength(A, MaxN-1+1, MaxN-1+1);
    
    //
    // Different problems
    //
    GPass:=1;
    while GPass<=GPassCount do
    begin
        
        //
        // zero matrix, several cases
        //
        I:=0;
        while I<=MaxN-1 do
        begin
            J:=0;
            while J<=MaxN-1 do
            begin
                A[I,J] := 0;
                Inc(J);
            end;
            Inc(I);
        end;
        I:=1;
        while I<=MaxN do
        begin
            TestProblem(A, I);
            Inc(I);
        end;
        
        //
        // Dense matrices
        //
        N:=1;
        while N<=MaxN do
        begin
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A[I,J] := 2*RandomReal-1;
                    Inc(J);
                end;
                Inc(I);
            end;
            TestProblem(A, N);
            Inc(N);
        end;
        
        //
        // Sparse matrices, very sparse matrices, incredible sparse matrices
        //
        N:=1;
        while N<=MaxN do
        begin
            FillSparseA(A, N, N, 0.8);
            TestProblem(A, N);
            FillSparseA(A, N, N, 0.9);
            TestProblem(A, N);
            FillSparseA(A, N, N, 0.95);
            TestProblem(A, N);
            Inc(N);
        end;
        Inc(GPass);
    end;
    
    //
    // report
    //
    WasErrors := DetErrors;
    if  not Silent then
    begin
        Write(Format('TESTING DET'#13#10'',[]));
        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;


(*************************************************************************
Sparse matrix
*************************************************************************)
procedure FillSparseA(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     Sparcity : Double);
var
    I : Integer;
    J : Integer;
begin
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            if RandomReal>=Sparcity then
            begin
                A[I,J] := 2*RandomReal-1;
            end
            else
            begin
                A[I,J] := 0;
            end;
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Problem testing
*************************************************************************)
procedure TestProblem(const A : TReal2DArray; N : Integer);
var
    I : Integer;
    J : Integer;
    B : TReal2DArray;
    C : TReal2DArray;
    Pivots : TInteger1DArray;
    V1 : Double;
    V2 : Double;
    VE : Double;
begin
    SetLength(B, N+1, N+1);
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            B[I,J] := A[I-1,J-1];
            Inc(J);
        end;
        Inc(I);
    end;
    VE := DetTriangle(B, N);
    V1 := RMatrixDet(A, N);
    MakeACopy(A, N, N, C);
    RMatrixLU(C, N, N, Pivots);
    V2 := RMatrixLUDet(C, Pivots, N);
    if VE<>0 then
    begin
        DetErrors := DetErrors or (AbsReal((V1-VE)/VE)>1.0E-9);
        DetErrors := DetErrors or (AbsReal((V2-VE)/VE)>1.0E-9);
    end
    else
    begin
        DetErrors := DetErrors or (V1<>VE);
        DetErrors := DetErrors or (V2<>VE);
    end;
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;


(*************************************************************************
Basic det subroutine
*************************************************************************)
function DetTriangle(A : TReal2DArray; const N : Integer):Double;
var
    i : Integer;
    j : Integer;
    k : Integer;
    l : Integer;
    f : Integer;
    z : Integer;
    t : Double;
    m1 : Double;
begin
    A := DynamicArrayCopy(A);
    Result := 1;
    k := 1;
    repeat
        m1 := 0;
        i := k;
        while i<=n do
        begin
            t := a[i,k];
            if AbsReal(t)>AbsReal(m1) then
            begin
                m1 := t;
                j := i;
            end;
            i := i+1;
        end;
        if AbsReal(m1)=0 then
        begin
            Result := 0;
            k := n+1;
        end
        else
        begin
            if j<>k then
            begin
                Result := -Result;
                l := k;
                while l<=n do
                begin
                    t := a[j,l];
                    a[j,l] := a[k,l];
                    a[k,l] := t;
                    l := l+1;
                end;
            end;
            f := k+1;
            while f<=n do
            begin
                t := a[f,k]/m1;
                z := k+1;
                while z<=n do
                begin
                    a[f,z] := a[f,z]-t*a[k,z];
                    z := z+1;
                end;
                f := f+1;
            end;
            Result := Result*a[k,k];
        end;
        k := k+1;
    until  not (k<=n);
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testdetunit_test_silent():Boolean;
begin
    Result := TestDet(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testdetunit_test():Boolean;
begin
    Result := TestDet(False);
end;


end.