testblasunit.pas

maths lib with source

    CTErrors := Was1;
    
    //
    // Testing MatrixMatrixMultiply
    //
    N := 10;
    SetLength(A, 2*N+1, 2*N+1);
    SetLength(B, 2*N+1, 2*N+1);
    SetLength(C1, 2*N+1, 2*N+1);
    SetLength(C2, 2*N+1, 2*N+1);
    SetLength(X1, N+1);
    SetLength(X2, N+1);
    I:=1;
    while I<=2*N do
    begin
        J:=1;
        while J<=2*N do
        begin
            A[I,J] := RandomReal;
            B[I,J] := RandomReal;
            Inc(J);
        end;
        Inc(I);
    end;
    PassCount := 1000;
    Was1 := False;
    Pass:=1;
    while Pass<=PassCount do
    begin
        I:=1;
        while I<=2*N do
        begin
            J:=1;
            while J<=2*N do
            begin
                C1[I,J] := 2.1*I+3.1*J;
                C2[I,J] := C1[I,J];
                Inc(J);
            end;
            Inc(I);
        end;
        L := 1+RandomInteger(N);
        K := 1+RandomInteger(N);
        R := 1+RandomInteger(N);
        I1 := 1+RandomInteger(N);
        J1 := 1+RandomInteger(N);
        I2 := 1+RandomInteger(N);
        J2 := 1+RandomInteger(N);
        I3 := 1+RandomInteger(N);
        J3 := 1+RandomInteger(N);
        Trans1 := RandomReal>0.5;
        Trans2 := RandomReal>0.5;
        if Trans1 then
        begin
            Col1 := L;
            Row1 := K;
        end
        else
        begin
            Col1 := K;
            Row1 := L;
        end;
        if Trans2 then
        begin
            Col2 := K;
            Row2 := R;
        end
        else
        begin
            Col2 := R;
            Row2 := K;
        end;
        Scl1 := RandomReal;
        Scl2 := RandomReal;
        MatrixMatrixMultiply(A, I1, I1+Row1-1, J1, J1+Col1-1, Trans1, B, I2, I2+Row2-1, J2, J2+Col2-1, Trans2, Scl1, C1, I3, I3+L-1, J3, J3+R-1, Scl2, X1);
        NaiveMatrixMatrixMultiply(A, I1, I1+Row1-1, J1, J1+Col1-1, Trans1, B, I2, I2+Row2-1, J2, J2+Col2-1, Trans2, Scl1, C2, I3, I3+L-1, J3, J3+R-1, Scl2);
        Err := 0;
        I:=1;
        while I<=L do
        begin
            J:=1;
            while J<=R do
            begin
                Err := Max(Err, AbsReal(C1[I3+I-1,J3+J-1]-C2[I3+I-1,J3+J-1]));
                Inc(J);
            end;
            Inc(I);
        end;
        if Err>100*MachineEpsilon then
        begin
            Was1 := True;
            Break;
        end;
        Inc(Pass);
    end;
    MMErrors := Was1;
    
    //
    // report
    //
    WasErrors := N2Errors or AMaxErrors or HSNErrors or MVErrors or ITErrors or CTErrors or MMErrors;
    if  not Silent then
    begin
        Write(Format('TESTING BLAS'#13#10'',[]));
        Write(Format('VectorNorm2:                             ',[]));
        if N2Errors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('AbsMax (vector/row/column):              ',[]));
        if AMaxErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('UpperHessenberg1Norm:                    ',[]));
        if HSNErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('MatrixVectorMultiply:                    ',[]));
        if MVErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('InplaceTranspose:                        ',[]));
        if ITErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('CopyAndTranspose:                        ',[]));
        if CTErrors then
        begin
            Write(Format('FAILED'#13#10'',[]));
        end
        else
        begin
            Write(Format('OK'#13#10'',[]));
        end;
        Write(Format('MatrixMatrixMultiply:                    ',[]));
        if MMErrors 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;


procedure NaiveMatrixMatrixMultiply(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;
     Alpha : Double;
     var C : TReal2DArray;
     CI1 : Integer;
     CI2 : Integer;
     CJ1 : Integer;
     CJ2 : Integer;
     Beta : Double);
var
    ARows : Integer;
    ACols : Integer;
    BRows : Integer;
    BCols : Integer;
    I : Integer;
    J : Integer;
    K : Integer;
    L : Integer;
    R : Integer;
    V : Double;
    X1 : TReal1DArray;
    X2 : TReal1DArray;
    i_ : Integer;
    i1_ : Integer;
begin
    
    //
    // 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, 'NaiveMatrixMatrixMultiply: incorrect matrix sizes!');
    if (ARows<=0) or (ACols<=0) or (BRows<=0) or (BCols<=0) then
    begin
        Exit;
    end;
    L := ARows;
    R := BCols;
    K := ACols;
    SetLength(X1, K+1);
    SetLength(X2, K+1);
    I:=1;
    while I<=L do
    begin
        J:=1;
        while J<=R do
        begin
            if  not TransA then
            begin
                if  not TransB then
                begin
                    i1_ := (AJ1)-(BI1);
                    V := 0.0;
                    for i_ := BI1 to BI2 do
                    begin
                        V := V + B[i_,BJ1+J-1]*A[AI1+I-1,i_+i1_];
                    end;
                end
                else
                begin
                    V := APVDotProduct(@B[BI1+J-1][0], BJ1, BJ2, @A[AI1+I-1][0], AJ1, AJ2);
                end;
            end
            else
            begin
                if  not TransB then
                begin
                    i1_ := (AI1)-(BI1);
                    V := 0.0;
                    for i_ := BI1 to BI2 do
                    begin
                        V := V + B[i_,BJ1+J-1]*A[i_+i1_,AJ1+I-1];
                    end;
                end
                else
                begin
                    i1_ := (AI1)-(BJ1);
                    V := 0.0;
                    for i_ := BJ1 to BJ2 do
                    begin
                        V := V + B[BI1+J-1,i_]*A[i_+i1_,AJ1+I-1];
                    end;
                end;
            end;
            if Beta=0 then
            begin
                C[CI1+I-1,CJ1+J-1] := Alpha*V;
            end
            else
            begin
                C[CI1+I-1,CJ1+J-1] := Beta*C[CI1+I-1,CJ1+J-1]+Alpha*V;
            end;
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testblasunit_test_silent():Boolean;
begin
    Result := TestBLAS(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testblasunit_test():Boolean;
begin
    Result := TestBLAS(False);
end;


end.