testldltunit.pas

maths lib with source

unit testldltunit;
interface
uses Math, Ap, Sysutils, ldlt;

function TestLDLT(Silent : Boolean):Boolean;
function testldltunit_test_silent():Boolean;
function testldltunit_test():Boolean;

implementation

procedure GenerateMatrix(var A : TReal2DArray;
     N : Integer;
     Task : Integer);forward;


function TestLDLT(Silent : Boolean):Boolean;
var
    A : TReal2DArray;
    A2 : TReal2DArray;
    L : TReal2DArray;
    D : TReal2DArray;
    U : TReal2DArray;
    T : TReal2DArray;
    T2 : TReal2DArray;
    P : TInteger1DArray;
    N : Integer;
    Pass : Integer;
    MTask : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    MinIJ : Integer;
    UpperIn : Boolean;
    CR : Boolean;
    V : Double;
    Err : Double;
    WasErrors : Boolean;
    PassCount : Integer;
    MaxN : Integer;
    HTask : Integer;
    Threshold : Double;
    i_ : Integer;
begin
    Err := 0;
    PassCount := 10;
    MaxN := 20;
    Threshold := 1000*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(L, N-1+1, N-1+1);
        SetLength(U, N-1+1, N-1+1);
        SetLength(D, N-1+1, N-1+1);
        SetLength(T, N-1+1, N-1+1);
        SetLength(T2, 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 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];
                            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;
                                end;
                            end
                            else
                            begin
                                if I<J then
                                begin
                                    A2[I,J] := 0;
                                end;
                            end;
                            Inc(J);
                        end;
                        Inc(I);
                    end;
                    
                    //
                    // LDLt
                    //
                    SMatrixLDLT(A2, N, UpperIn, P);
                    
                    //
                    // Test (upper or lower)
                    //
                    if UpperIn then
                    begin
                        
                        //
                        // Unpack D
                        //
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=N-1 do
                            begin
                                D[I,J] := 0;
                                Inc(J);
                            end;
                            Inc(I);
                        end;
                        K := 0;
                        while K<N do
                        begin
                            if P[K]>=0 then
                            begin
                                D[K,K] := A2[K,K];
                                K := K+1;
                            end
                            else
                            begin
                                D[K,K] := A2[K,K];
                                D[K,K+1] := A2[K,K+1];
                                D[K+1,K] := A2[K,K+1];
                                D[K+1,K+1] := A2[K+1,K+1];
                                K := K+2;
                            end;
                        end;
                        
                        //
                        // Unpack U
                        //
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=N-1 do
                            begin
                                U[I,J] := 0;
                                Inc(J);
                            end;
                            U[I,I] := 1;
                            Inc(I);
                        end;
                        K := 0;
                        while K<N do
                        begin
                            
                            //
                            // unpack Uk
                            //
                            I:=0;
                            while I<=N-1 do
                            begin
                                J:=0;
                                while J<=N-1 do
                                begin
                                    T[I,J] := 0;
                                    Inc(J);
                                end;
                                T[I,I] := 1;
                                Inc(I);
                            end;
                            if P[K]>=0 then
                            begin
                                I:=0;
                                while I<=K-1 do
                                begin
                                    T[I,K] := A2[I,K];
                                    Inc(I);
                                end;
                            end
                            else
                            begin
                                I:=0;
                                while I<=K-1 do
                                begin
                                    T[I,K] := A2[I,K];
                                    T[I,K+1] := A2[I,K+1];
                                    Inc(I);
                                end;
                            end;
                            
                            //
                            // multiply U
                            //
                            I:=0;
                            while I<=N-1 do
                            begin
                                J:=0;
                                while J<=N-1 do
                                begin
                                    V := 0.0;
                                    for i_ := 0 to N-1 do
                                    begin
                                        V := V + T[I,i_]*U[i_,J];
                                    end;
                                    T2[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
                                    U[I,J] := T2[I,J];
                                    Inc(J);
                                end;
                                Inc(I);
                            end;
                            
                            //
                            // permutations
                            //
                            if P[K]>=0 then
                            begin
                                J:=0;
                                while J<=N-1 do
                                begin
                                    V := U[K,J];
                                    U[K,J] := U[P[K],J];
                                    U[P[K],J] := V;
                                    Inc(J);
                                end;
                                K := K+1;
                            end
                            else
                            begin
                                J:=0;
                                while J<=N-1 do
                                begin
                                    V := U[K,J];
                                    U[K,J] := U[N+P[K],J];
                                    U[N+P[K],J] := V;
                                    Inc(J);
                                end;
                                K := K+2;
                            end;
                        end;
                        
                        //
                        // Calculate U*D*U', store result in T2
                        //
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=N-1 do
                            begin
                                V := 0.0;
                                for i_ := 0 to N-1 do
                                begin
                                    V := V + U[I,i_]*D[i_,J];
                                end;
                                T[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
                                V := APVDotProduct(@T[I][0], 0, N-1, @U[J][0], 0, N-1);
                                T2[I,J] := V;
                                Inc(J);
                            end;
                            Inc(I);
                        end;
                    end
                    else
                    begin
                        
                        //
                        // Unpack D
                        //
                        I:=0;
                        while I<=N-1 do
                        begin
                            J:=0;
                            while J<=N-1 do
                            begin
                                D[I,J] := 0;
                                Inc(J);
                            end;
                            Inc(I);
                        end;