testbdunit.pas

maths lib with source

unit testbdunit;
interface
uses Math, Ap, Sysutils, reflections, bidiagonal;

function TestBD(Silent : Boolean):Boolean;
function testbdunit_test_silent():Boolean;
function testbdunit_test():Boolean;

implementation

var
    DecompErrors : Boolean;
    PropErrors : Boolean;
    PartErrors : Boolean;
    MulErrors : Boolean;
    Threshold : Double;

procedure FillSparseA(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     Sparcity : Double);forward;
procedure TestProblem(const A : TReal2DArray;
     M : Integer;
     N : Integer);forward;
procedure MakeACopy(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);forward;
procedure MakeACopyOldMem(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);forward;
function MatrixDiff(const A : TReal2DArray;
     const B : TReal2DArray;
     M : Integer;
     N : Integer):Double;forward;
procedure InternalMatrixMatrixMultiply(const A : TReal2DArray;
     AI1 : Integer;
     AI2 : Integer;
     AJ1 : Integer;
     AJ2 : Integer;
     TransA : Boolean;
     const B : TReal2DArray;
     BI1 : Integer;
     BI2 : Integer;
     BJ1 : Integer;
     BJ2 : Integer;
     TransB : Boolean;
     var C : TReal2DArray;
     CI1 : Integer;
     CI2 : Integer;
     CJ1 : Integer;
     CJ2 : Integer);forward;


(*************************************************************************
Main unittest subroutine
*************************************************************************)
function TestBD(Silent : Boolean):Boolean;
var
    ShortMN : Integer;
    MaxMN : Integer;
    GPassCount : Integer;
    A : TReal2DArray;
    M : Integer;
    N : Integer;
    GPass : Integer;
    I : Integer;
    J : Integer;
    WasErrors : Boolean;
begin
    DecompErrors := False;
    PropErrors := False;
    MulErrors := False;
    PartErrors := False;
    Threshold := 5*100*MachineEpsilon;
    WasErrors := False;
    ShortMN := 5;
    MaxMN := 10;
    GPassCount := 5;
    SetLength(A, MaxMN-1+1, MaxMN-1+1);
    
    //
    // Different problems
    //
    GPass:=1;
    while GPass<=GPassCount do
    begin
        
        //
        // zero matrix, several cases
        //
        I:=0;
        while I<=MaxMN-1 do
        begin
            J:=0;
            while J<=MaxMN-1 do
            begin
                A[I,J] := 0;
                Inc(J);
            end;
            Inc(I);
        end;
        I:=1;
        while I<=MaxMN do
        begin
            J:=1;
            while J<=MaxMN do
            begin
                TestProblem(A, I, J);
                Inc(J);
            end;
            Inc(I);
        end;
        
        //
        // Long dense matrix
        //
        I:=0;
        while I<=MaxMN-1 do
        begin
            J:=0;
            while J<=ShortMN-1 do
            begin
                A[I,J] := 2*RandomReal-1;
                Inc(J);
            end;
            Inc(I);
        end;
        I:=ShortMN+1;
        while I<=MaxMN do
        begin
            TestProblem(A, I, ShortMN);
            Inc(I);
        end;
        I:=0;
        while I<=ShortMN-1 do
        begin
            J:=0;
            while J<=MaxMN-1 do
            begin
                A[I,J] := 2*RandomReal-1;
                Inc(J);
            end;
            Inc(I);
        end;
        J:=ShortMN+1;
        while J<=MaxMN do
        begin
            TestProblem(A, ShortMN, J);
            Inc(J);
        end;
        
        //
        // Dense matrices
        //
        M:=1;
        while M<=MaxMN do
        begin
            N:=1;
            while N<=MaxMN do
            begin
                I:=0;
                while I<=M-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, M, N);
                Inc(N);
            end;
            Inc(M);
        end;
        
        //
        // Sparse matrices, very sparse matrices, incredible sparse matrices
        //
        M:=1;
        while M<=MaxMN do
        begin
            N:=1;
            while N<=MaxMN do
            begin
                FillSparseA(A, M, N, 0.8);
                TestProblem(A, M, N);
                FillSparseA(A, M, N, 0.9);
                TestProblem(A, M, N);
                FillSparseA(A, M, N, 0.95);
                TestProblem(A, M, N);
                Inc(N);
            end;
            Inc(M);
        end;
        Inc(GPass);
    end;
    
    //
    // report
    //
    WasErrors := DecompErrors or PropErrors or PartErrors or MulErrors;
    if  not Silent then
    begin
        Write(Format('TESTING 2BIDIAGONAL'#13#10'',[]));
        Write(Format('DECOMPOSITION ERRORS                     ',[]));
        if DecompErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('MATRIX PROPERTIES                        ',[]));
        if PropErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('PARTIAL UNPACKING                        ',[]));
        if PartErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('MULTIPLICATION TEST                      ',[]));
        if MulErrors 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;


(*************************************************************************
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; M : Integer; N : Integer);
var
    I : Integer;
    J : Integer;
    K : Integer;
    MX : Double;
    T : TReal2DArray;
    PT : TReal2DArray;
    Q : TReal2DArray;
    R : TReal2DArray;
    BD : TReal2DArray;
    X : TReal2DArray;
    R1 : TReal2DArray;
    R2 : TReal2DArray;
    TauP : TReal1DArray;
    TauQ : TReal1DArray;
    D : TReal1DArray;
    E : TReal1DArray;
    Up : Boolean;
    V : Double;
    MTSize : Integer;
    i_ : Integer;
begin
    
    //
    // MX - estimate of the matrix norm
    //
    MX := 0;
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            if AbsReal(A[I,J])>MX then
            begin
                MX := AbsReal(A[I,J]);
            end;
            Inc(J);
        end;
        Inc(I);
    end;
    if MX=0 then
    begin
        MX := 1;
    end;
    
    //
    // Bidiagonal decomposition error
    //
    MakeACopy(A, M, N, T);
    RMatrixBD(T, M, N, TauQ, TauP);
    RMatrixBDUnpackQ(T, M, N, TauQ, M, Q);
    RMatrixBDUnpackPT(T, M, N, TauP, N, PT);
    RMatrixBDUnpackDiagonals(T, M, N, Up, D, E);
    SetLength(BD, M-1+1, N-1+1);
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            BD[I,J] := 0;
            Inc(J);
        end;
        Inc(I);
    end;
    I:=0;
    while I<=Min(M, N)-1 do
    begin
        BD[I,I] := D[I];
        Inc(I);
    end;
    if Up then
    begin
        I:=0;
        while I<=Min(M, N)-2 do
        begin
            BD[I,I+1] := E[I];
            Inc(I);
        end;
    end
    else
    begin
        I:=0;
        while I<=Min(M, N)-2 do
        begin
            BD[I+1,I] := E[I];
            Inc(I);
        end;
    end;
    SetLength(R, M-1+1, N-1+1);
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := 0.0;
            for i_ := 0 to M-1 do
            begin
                V := V + Q[I,i_]*BD[i_,J];
            end;
            R[I,J] := V;
            Inc(J);
        end;
        Inc(I);
    end;
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := 0.0;
            for i_ := 0 to N-1 do
            begin
                V := V + R[I,i_]*PT[i_,J];
            end;
            DecompErrors := DecompErrors or (AbsReal(V-A[I,J])>Threshold);
            Inc(J);
        end;
        Inc(I);
    end;
    
    //
    // Orthogonality test for Q/PT
    //
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=M-1 do
        begin
            V := 0.0;
            for i_ := 0 to M-1 do
            begin
                V := V + Q[i_,I]*Q[i_,J];
            end;
            if I=J then
            begin
                PropErrors := PropErrors or (AbsReal(V-1)>Threshold);
            end
            else
            begin
                PropErrors := PropErrors or (AbsReal(V)>Threshold);
            end;
            Inc(J);
        end;
        Inc(I);
    end;
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := APVDotProduct(@PT[I][0], 0, N-1, @PT[J][0], 0, N-1);
            if I=J then
            begin
                PropErrors := PropErrors or (AbsReal(V-1)>Threshold);
            end
            else
            begin
                PropErrors := PropErrors or (AbsReal(V)>Threshold);
            end;