testcinvunit.pas

maths lib with source

unit testcinvunit;
interface
uses Math, Ap, Sysutils, clu, ctrinverse, cinverse;

function TestCInv(Silent : Boolean):Boolean;
function testcinvunit_test_silent():Boolean;
function testcinvunit_test():Boolean;

implementation

var
    InvErrors : Boolean;
    Threshold : Double;

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


(*************************************************************************
Main unittest subroutine
*************************************************************************)
function TestCInv(Silent : Boolean):Boolean;
var
    MaxN : Integer;
    GPassCount : Integer;
    A : TComplex2DArray;
    N : Integer;
    GPass : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    WasErrors : Boolean;
begin
    InvErrors := False;
    Threshold := 5*1000*MachineEpsilon;
    WasErrors := False;
    MaxN := 10;
    GPassCount := 30;
    
    //
    // Different problems
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(A, N-1+1, N-1+1);
        GPass:=1;
        while GPass<=GPassCount do
        begin
            
            //
            // diagonal matrix, several cases
            //
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A[I,J] := C_Complex(0);
                    Inc(J);
                end;
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                A[I,I].X := 2*RandomReal-1;
                A[I,I].Y := 2*RandomReal-1;
                Inc(I);
            end;
            TestProblem(A, N);
            
            //
            // shifted diagonal matrix, several cases
            //
            K := RandomInteger(N);
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A[I,J] := C_Complex(0);
                    Inc(J);
                end;
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                A[I,(I+K) mod N].X := 2*RandomReal-1;
                A[I,(I+K) mod N].Y := 2*RandomReal-1;
                Inc(I);
            end;
            TestProblem(A, N);
            
            //
            // Dense matrices
            //
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A[I,J].X := 2*RandomReal-1;
                    A[I,J].Y := 2*RandomReal-1;
                    Inc(J);
                end;
                Inc(I);
            end;
            TestProblem(A, N);
            Inc(GPass);
        end;
        Inc(N);
    end;
    
    //
    // report
    //
    WasErrors := InvErrors;
    if  not Silent then
    begin
        Write(Format('TESTING COMPLEX INVERSE'#13#10'',[]));
        Write(Format('INVERSE ERRORS                           ',[]));
        if InvErrors 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;


(*************************************************************************
Problem testing
*************************************************************************)
procedure TestProblem(const A : TComplex2DArray; N : Integer);
var
    B : TComplex2DArray;
    BLU : TComplex2DArray;
    T1 : TComplex2DArray;
    P : TInteger1DArray;
    I : Integer;
    J : Integer;
    V : Complex;
begin
    
    //
    // Decomposition
    //
    MakeACopy(A, N, N, B);
    CMatrixInverse(B, N);
    MakeACopy(A, N, N, BLU);
    CMatrixLU(BLU, N, N, P);
    CMatrixLUInverse(BLU, P, N);
    
    //
    // Test
    //
    SetLength(T1, N-1+1, N-1+1);
    InternalMatrixMatrixMultiply(A, 0, N-1, 0, N-1, False, B, 0, N-1, 0, N-1, False, T1, 0, N-1, 0, N-1);
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := T1[I,J];
            if I=J then
            begin
                V := C_SubR(V,1);
            end;
            InvErrors := InvErrors or (AbsComplex(V)>Threshold);
            Inc(J);
        end;
        Inc(I);
    end;
    InternalMatrixMatrixMultiply(A, 0, N-1, 0, N-1, False, BLU, 0, N-1, 0, N-1, False, T1, 0, N-1, 0, N-1);
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := T1[I,J];
            if I=J then
            begin
                V := C_SubR(V,1);
            end;
            InvErrors := InvErrors or (AbsComplex(V)>Threshold);
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Copy
*************************************************************************)
procedure MakeACopy(const A : TComplex2DArray;
     M : Integer;
     N : Integer;
     var B : TComplex2DArray);
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;


(*************************************************************************
Copy
*************************************************************************)
procedure MakeACopyOldMem(const A : TComplex2DArray;
     M : Integer;
     N : Integer;
     var B : TComplex2DArray);
var
    I : Integer;
    J : Integer;
begin
    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;


(*************************************************************************
Diff
*************************************************************************)
function MatrixDiff(const A : TComplex2DArray;
     const B : TComplex2DArray;
     M : Integer;
     N : Integer):Double;
var
    I : Integer;
    J : Integer;
begin
    Result := 0;
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            Result := Max(Result, AbsComplex(C_Sub(B[I,J],A[I,J])));
            Inc(J);
        end;