testdetldltunit.pas

maths lib with source

unit testdetldltunit;
interface
uses Math, Ap, Sysutils, ldlt, sdet;

function TestDetLDLT(Silent : Boolean):Boolean;
function testdetldltunit_test_silent():Boolean;
function testdetldltunit_test():Boolean;

implementation

function RefDet(const A0 : TReal2DArray; N : Integer):Double;forward;
procedure GenerateMatrix(var A : TReal2DArray;
     N : Integer;
     Task : Integer);forward;
procedure RestoreMatrix(var A : TReal2DArray;
     N : Integer;
     UpperIn : Boolean);forward;


function TestDetLDLT(Silent : Boolean):Boolean;
var
    A : TReal2DArray;
    A2 : TReal2DArray;
    A3 : TReal2DArray;
    P : TInteger1DArray;
    N : Integer;
    Pass : Integer;
    MTask : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    MinIJ : Integer;
    UpperIn : Boolean;
    CR : Boolean;
    V1 : Double;
    V2 : Double;
    Err : Double;
    WasErrors : Boolean;
    PassCount : Integer;
    MaxN : Integer;
    HTask : Integer;
    Threshold : Double;
begin
    Err := 0;
    PassCount := 10;
    MaxN := 20;
    Threshold := 100000*MachineEpsilon;
    WasErrors := False;
    
    //
    // Test
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(A, N-1+1, N-1+1);
        SetLength(A2, N-1+1, N-1+1);
        SetLength(A3, N-1+1, N-1+1);
        MTask:=0;
        while MTask<=2 do
        begin
            HTask:=0;
            while HTask<=1 do
            begin
                Pass:=1;
                while Pass<=PassCount do
                begin
                    UpperIn := HTask=0;
                    
                    //
                    // Prepare task:
                    // * A contains symmetric matrix
                    // * A2, A3 contains its upper (or lower) half
                    //
                    GenerateMatrix(A, N, MTask);
                    I:=0;
                    while I<=N-1 do
                    begin
                        J:=0;
                        while J<=N-1 do
                        begin
                            A2[I,J] := A[I,J];
                            A3[I,J] := A[I,J];
                            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
                                    A2[I,J] := 0;
                                    A3[I,J] := 0;
                                end;
                            end
                            else
                            begin
                                if I<J then
                                begin
                                    A2[I,J] := 0;
                                    A3[I,J] := 0;
                                end;
                            end;
                            Inc(J);
                        end;
                        Inc(I);
                    end;
                    
                    //
                    // Test 1: det(A2)
                    //
                    V1 := SMatrixDet(A2, N, UpperIn);
                    V2 := RefDet(A, N);
                    if V2<>0 then
                    begin
                        Err := Max(Err, AbsReal((V1-V2)/V2));
                    end
                    else
                    begin
                        Err := Max(Err, AbsReal(V1));
                    end;
                    
                    //
                    // Test 2: inv(LDLt(A3))
                    //
                    SMatrixLDLT(A3, N, UpperIn, P);
                    V1 := SMatrixLDLTDet(A3, P, N, UpperIn);
                    V2 := RefDet(A, N);
                    if V2<>0 then
                    begin
                        Err := Max(Err, AbsReal((V1-V2)/V2));
                    end
                    else
                    begin
                        Err := Max(Err, AbsReal(V1));
                    end;
                    Inc(Pass);
                end;
                Inc(HTask);
            end;
            Inc(MTask);
        end;
        Inc(N);
    end;
    
    //
    // report
    //
    WasErrors := Err>Threshold;
    if  not Silent then
    begin
        Write(Format('TESTING LDLT DET'#13#10'',[]));
        Write(Format('ERROR:                                   %5.4e'#13#10'',[
            Err]));
        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;


(*************************************************************************
Reference det subroutine
*************************************************************************)
function RefDet(const A0 : TReal2DArray; N : Integer):Double;
var
    i : Integer;
    j : Integer;
    k : Integer;
    l : Integer;
    f : Integer;
    z : Integer;
    t : Double;
    m1 : Double;
    A : TReal2DArray;
begin
    SetLength(A, N+1, N+1);
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            A[I,J] := A0[I-1,J-1];
            Inc(J);
        end;
        Inc(I);
    end;
    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;


procedure GenerateMatrix(var A : TReal2DArray; N : Integer; Task : Integer);
var
    I : Integer;
    J : Integer;
begin
    if Task=0 then
    begin
        
        //
        // Zero matrix
        //
        I:=0;
        while I<=N-1 do
        begin
            J:=0;
            while J<=N-1 do
            begin
                A[I,J] := 0;
                Inc(J);
            end;
            Inc(I);
        end;
    end;
    if Task=1 then
    begin
        
        //
        // Sparse matrix
        //
        I:=0;
        while I<=N-1 do
        begin
            J:=I+1;
            while J<=N-1 do
            begin
                if RandomReal>0.95 then
                begin
                    A[I,J] := 2*RandomReal-1;
                end
                else
                begin
                    A[I,J] := 0;
                end;
                A[J,I] := A[I,J];
                Inc(J);
            end;
            if RandomReal>0.95 then
            begin
                A[I,I] := (2*RandomInteger(2)-1)*(0.8+RandomReal);
            end
            else
            begin
                A[I,I] := 0;
            end;
            Inc(I);
        end;
    end;
    if Task=2 then
    begin
        
        //
        // Dense matrix
        //
        I:=0;
        while I<=N-1 do
        begin
            J:=I+1;
            while J<=N-1 do
            begin
                A[I,J] := 2*RandomReal-1;
                A[J,I] := A[I,J];
                Inc(J);
            end;
            A[I,I] := (2*RandomInteger(2)-1)*(0.7+RandomReal);
            Inc(I);
        end;
    end;
end;


procedure RestoreMatrix(var A : TReal2DArray; N : Integer; UpperIn : Boolean);
var
    I : Integer;
    J : Integer;
begin
    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
                    A[I,J] := A[J,I];
                end;
            end
            else
            begin
                if I<J then
                begin
                    A[I,J] := A[J,I];
                end;
            end;
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testdetldltunit_test_silent():Boolean;
begin
    Result := TestDetLDLT(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testdetldltunit_test():Boolean;
begin
    Result := TestDetLDLT(False);
end;


end.