testbdunit.pas

maths lib with source

            Inc(J);
        end;
        Inc(I);
    end;
    
    //
    // Partial unpacking test
    //
    K:=1;
    while K<=M-1 do
    begin
        RMatrixBDUnpackQ(T, M, N, TauQ, K, R);
        I:=0;
        while I<=M-1 do
        begin
            J:=0;
            while J<=K-1 do
            begin
                PartErrors := PartErrors or (AbsReal(R[I,J]-Q[I,J])>10*MachineEpsilon);
                Inc(J);
            end;
            Inc(I);
        end;
        Inc(K);
    end;
    K:=1;
    while K<=N-1 do
    begin
        RMatrixBDUnpackPT(T, M, N, TauP, K, R);
        I:=0;
        while I<=K-1 do
        begin
            J:=0;
            while J<=N-1 do
            begin
                PartErrors := PartErrors or (R[I,J]-PT[I,J]<>0);
                Inc(J);
            end;
            Inc(I);
        end;
        Inc(K);
    end;
    
    //
    // Multiplication test
    //
    SetLength(X, Max(M, N)-1+1, Max(M, N)-1+1);
    SetLength(R, Max(M, N)-1+1, Max(M, N)-1+1);
    SetLength(R1, Max(M, N)-1+1, Max(M, N)-1+1);
    SetLength(R2, Max(M, N)-1+1, Max(M, N)-1+1);
    I:=0;
    while I<=Max(M, N)-1 do
    begin
        J:=0;
        while J<=Max(M, N)-1 do
        begin
            X[I,J] := 2*RandomReal-1;
            Inc(J);
        end;
        Inc(I);
    end;
    MTSize := 1+RandomInteger(Max(M, N));
    MakeACopyOldMem(X, MTSize, M, R);
    InternalMatrixMatrixMultiply(R, 0, MTSize-1, 0, M-1, False, Q, 0, M-1, 0, M-1, False, R1, 0, MTSize-1, 0, M-1);
    MakeACopyOldMem(X, MTSize, M, R2);
    RMatrixBDMultiplyByQ(T, M, N, TauQ, R2, MTSize, M, True, False);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, MTSize, M)>Threshold);
    MakeACopyOldMem(X, MTSize, M, R);
    InternalMatrixMatrixMultiply(R, 0, MTSize-1, 0, M-1, False, Q, 0, M-1, 0, M-1, True, R1, 0, MTSize-1, 0, M-1);
    MakeACopyOldMem(X, MTSize, M, R2);
    RMatrixBDMultiplyByQ(T, M, N, TauQ, R2, MTSize, M, True, True);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, MTSize, M)>Threshold);
    MakeACopyOldMem(X, M, MTSize, R);
    InternalMatrixMatrixMultiply(Q, 0, M-1, 0, M-1, False, R, 0, M-1, 0, MTSize-1, False, R1, 0, M-1, 0, MTSize-1);
    MakeACopyOldMem(X, M, MTSize, R2);
    RMatrixBDMultiplyByQ(T, M, N, TauQ, R2, M, MTSize, False, False);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, M, MTSize)>Threshold);
    MakeACopyOldMem(X, M, MTSize, R);
    InternalMatrixMatrixMultiply(Q, 0, M-1, 0, M-1, True, R, 0, M-1, 0, MTSize-1, False, R1, 0, M-1, 0, MTSize-1);
    MakeACopyOldMem(X, M, MTSize, R2);
    RMatrixBDMultiplyByQ(T, M, N, TauQ, R2, M, MTSize, False, True);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, M, MTSize)>Threshold);
    MakeACopyOldMem(X, MTSize, N, R);
    InternalMatrixMatrixMultiply(R, 0, MTSize-1, 0, N-1, False, PT, 0, N-1, 0, N-1, True, R1, 0, MTSize-1, 0, N-1);
    MakeACopyOldMem(X, MTSize, N, R2);
    RMatrixBDMultiplyByP(T, M, N, TauP, R2, MTSize, N, True, False);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, MTSize, N)>Threshold);
    MakeACopyOldMem(X, MTSize, N, R);
    InternalMatrixMatrixMultiply(R, 0, MTSize-1, 0, N-1, False, PT, 0, N-1, 0, N-1, False, R1, 0, MTSize-1, 0, N-1);
    MakeACopyOldMem(X, MTSize, N, R2);
    RMatrixBDMultiplyByP(T, M, N, TauP, R2, MTSize, N, True, True);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, MTSize, N)>Threshold);
    MakeACopyOldMem(X, N, MTSize, R);
    InternalMatrixMatrixMultiply(PT, 0, N-1, 0, N-1, True, R, 0, N-1, 0, MTSize-1, False, R1, 0, N-1, 0, MTSize-1);
    MakeACopyOldMem(X, N, MTSize, R2);
    RMatrixBDMultiplyByP(T, M, N, TauP, R2, N, MTSize, False, False);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, N, MTSize)>Threshold);
    MakeACopyOldMem(X, N, MTSize, R);
    InternalMatrixMatrixMultiply(PT, 0, N-1, 0, N-1, False, R, 0, N-1, 0, MTSize-1, False, R1, 0, N-1, 0, MTSize-1);
    MakeACopyOldMem(X, N, MTSize, R2);
    RMatrixBDMultiplyByP(T, M, N, TauP, R2, N, MTSize, False, True);
    MulErrors := MulErrors or (MatrixDiff(R1, R2, N, MTSize)>Threshold);
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;


(*************************************************************************
Copy
*************************************************************************)
procedure MakeACopyOldMem(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);
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 : TReal2DArray;
     const B : TReal2DArray;
     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, AbsReal(B[I,J]-A[I,J]));
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Matrix multiplication
*************************************************************************)
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);
var
    ARows : Integer;
    ACols : Integer;
    BRows : Integer;
    BCols : Integer;
    CRows : Integer;
    CCols : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    L : Integer;
    R : Integer;
    V : Double;
    WORK : TReal1DArray;
    Beta : Double;
    Alpha : Double;
    i_ : Integer;
    i1_ : Integer;
begin
    
    //
    // Pre-setup
    //
    K := Max(AI2-AI1+1, AJ2-AJ1+1);
    K := Max(K, BI2-BI1+1);
    K := Max(K, BJ2-BJ1+1);
    SetLength(WORK, K+1);
    Beta := 0;
    Alpha := 1;
    
    //
    // Setup
    //
    if  not TransA then
    begin
        ARows := AI2-AI1+1;
        ACols := AJ2-AJ1+1;
    end
    else
    begin
        ARows := AJ2-AJ1+1;
        ACols := AI2-AI1+1;
    end;
    if  not TransB then
    begin
        BRows := BI2-BI1+1;
        BCols := BJ2-BJ1+1;
    end
    else
    begin
        BRows := BJ2-BJ1+1;
        BCols := BI2-BI1+1;
    end;
    Assert(ACols=BRows, 'MatrixMatrixMultiply: incorrect matrix sizes!');
    if (ARows<=0) or (ACols<=0) or (BRows<=0) or (BCols<=0) then
    begin
        Exit;
    end;
    CRows := ARows;
    CCols := BCols;
    
    //
    // Test WORK
    //
    I := Max(ARows, ACols);
    I := Max(BRows, I);
    I := Max(I, BCols);
    Work[1] := 0;
    Work[I] := 0;
    
    //
    // Prepare C
    //
    if Beta=0 then
    begin
        I:=CI1;
        while I<=CI2 do
        begin
            J:=CJ1;
            while J<=CJ2 do
            begin
                C[I,J] := 0;
                Inc(J);
            end;
            Inc(I);
        end;
    end
    else
    begin
        I:=CI1;
        while I<=CI2 do
        begin
            APVMul(@C[I][0], CJ1, CJ2, Beta);
            Inc(I);
        end;
    end;
    
    //
    // A*B
    //
    if  not TransA and  not TransB then
    begin
        L:=AI1;
        while L<=AI2 do
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                V := Alpha*A[L,AJ1+R-BI1];
                K := CI1+L-AI1;
                APVAdd(@C[K][0], CJ1, CJ2, @B[R][0], BJ1, BJ2, V);
                Inc(R);
            end;
            Inc(L);
        end;
        Exit;
    end;
    
    //
    // A*B'
    //
    if  not TransA and TransB then
    begin
        if ARows*ACols<BRows*BCols then
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                L:=AI1;
                while L<=AI2 do
                begin
                    V := APVDotProduct(@A[L][0], AJ1, AJ2, @B[R][0], BJ1, BJ2);
                    C[CI1+L-AI1,CJ1+R-BI1] := C[CI1+L-AI1,CJ1+R-BI1]+Alpha*V;
                    Inc(L);
                end;
                Inc(R);
            end;
            Exit;
        end
        else
        begin
            L:=AI1;
            while L<=AI2 do
            begin
                R:=BI1;
                while R<=BI2 do
                begin
                    V := APVDotProduct(@A[L][0], AJ1, AJ2, @B[R][0], BJ1, BJ2);
                    C[CI1+L-AI1,CJ1+R-BI1] := C[CI1+L-AI1,CJ1+R-BI1]+Alpha*V;
                    Inc(R);
                end;
                Inc(L);
            end;
            Exit;
        end;
    end;
    
    //
    // A'*B
    //
    if TransA and  not TransB then
    begin
        L:=AJ1;
        while L<=AJ2 do
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                V := Alpha*A[AI1+R-BI1,L];
                K := CI1+L-AJ1;
                APVAdd(@C[K][0], CJ1, CJ2, @B[R][0], BJ1, BJ2, V);
                Inc(R);
            end;
            Inc(L);
        end;
        Exit;
    end;
    
    //
    // A'*B'
    //
    if TransA and TransB then
    begin
        if ARows*ACols<BRows*BCols then
        begin
            R:=BI1;
            while R<=BI2 do
            begin
                I:=1;
                while I<=CRows do
                begin
                    WORK[I] := 0.0;
                    Inc(I);
                end;
                L:=AI1;
                while L<=AI2 do
                begin
                    V := Alpha*B[R,BJ1+L-AI1];
                    K := CJ1+R-BI1;
                    APVAdd(@WORK[0], 1, CRows, @A[L][0], AJ1, AJ2, V);
                    Inc(L);
                end;
                i1_ := (1) - (CI1);
                for i_ := CI1 to CI2 do
                begin
                    C[i_,K] := C[i_,K] + WORK[i_+i1_];
                end;
                Inc(R);
            end;
            Exit;
        end
        else
        begin
            L:=AJ1;
            while L<=AJ2 do
            begin
                K := AI2-AI1+1;
                i1_ := (AI1) - (1);
                for i_ := 1 to K do
                begin
                    WORK[i_] := A[i_+i1_,L];
                end;
                R:=BI1;
                while R<=BI2 do
                begin
                    V := APVDotProduct(@WORK[0], 1, K, @B[R][0], BJ1, BJ2);
                    C[CI1+L-AJ1,CJ1+R-BI1] := C[CI1+L-AJ1,CJ1+R-BI1]+Alpha*V;
                    Inc(R);
                end;
                Inc(L);
            end;
            Exit;
        end;
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testbdunit_test_silent():Boolean;
begin
    Result := TestBD(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testbdunit_test():Boolean;
begin
    Result := TestBD(False);
end;


end.