testblasunit.pas

maths lib with source

unit testblasunit;
interface
uses Math, Ap, Sysutils, blas;

function TestBLAS(Silent : Boolean):Boolean;
function testblasunit_test_silent():Boolean;
function testblasunit_test():Boolean;

implementation

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);forward;


function TestBLAS(Silent : Boolean):Boolean;
var
    Pass : Integer;
    PassCount : Integer;
    N : Integer;
    I : Integer;
    I1 : Integer;
    I2 : Integer;
    J : Integer;
    J1 : Integer;
    J2 : Integer;
    L : Integer;
    K : Integer;
    R : Integer;
    I3 : Integer;
    J3 : Integer;
    Col1 : Integer;
    Col2 : Integer;
    Row1 : Integer;
    Row2 : Integer;
    X1 : TReal1DArray;
    X2 : TReal1DArray;
    A : TReal2DArray;
    B : TReal2DArray;
    C1 : TReal2DArray;
    C2 : TReal2DArray;
    Err : Double;
    E1 : Double;
    E2 : Double;
    E3 : Double;
    V : Double;
    Scl1 : Double;
    Scl2 : Double;
    Scl3 : Double;
    Was1 : Boolean;
    Was2 : Boolean;
    Trans1 : Boolean;
    Trans2 : Boolean;
    N2Errors : Boolean;
    HSNErrors : Boolean;
    AMaxErrors : Boolean;
    MVErrors : Boolean;
    ITErrors : Boolean;
    CTErrors : Boolean;
    MMErrors : Boolean;
    WasErrors : Boolean;
    i_ : Integer;
begin
    N2Errors := False;
    AMaxErrors := False;
    HSNErrors := False;
    MVErrors := False;
    ITErrors := False;
    CTErrors := False;
    MMErrors := False;
    WasErrors := False;
    
    //
    // Test Norm2
    //
    PassCount := 1000;
    E1 := 0;
    E2 := 0;
    E3 := 0;
    Scl2 := 0.5*MaxRealNumber;
    Scl3 := 2*MinRealNumber;
    Pass:=1;
    while Pass<=PassCount do
    begin
        N := 1+RandomInteger(1000);
        I1 := RandomInteger(10);
        I2 := N+I1-1;
        SetLength(X1, I2+1);
        SetLength(X2, I2+1);
        I:=I1;
        while I<=I2 do
        begin
            X1[I] := 2*RandomReal-1;
            Inc(I);
        end;
        V := 0;
        I:=I1;
        while I<=I2 do
        begin
            V := V+Sqr(X1[I]);
            Inc(I);
        end;
        V := Sqrt(V);
        E1 := Max(E1, AbsReal(V-VectorNorm2(X1, I1, I2)));
        I:=I1;
        while I<=I2 do
        begin
            X2[I] := Scl2*X1[I];
            Inc(I);
        end;
        E2 := Max(E2, AbsReal(V*Scl2-VectorNorm2(X2, I1, I2)));
        I:=I1;
        while I<=I2 do
        begin
            X2[I] := Scl3*X1[I];
            Inc(I);
        end;
        E3 := Max(E3, AbsReal(V*Scl3-VectorNorm2(X2, I1, I2)));
        Inc(Pass);
    end;
    E2 := E2/Scl2;
    E3 := E3/Scl3;
    N2Errors := (E1>=10000*MachineEpsilon) or (E2>=10000*MachineEpsilon) or (E3>=10000*MachineEpsilon);
    
    //
    // Testing VectorAbsMax, Column/Row AbsMax
    //
    SetLength(X1, 5+1);
    X1[1] := 2.0;
    X1[2] := 0.2;
    X1[3] := -1.3;
    X1[4] := 0.7;
    X1[5] := -3.0;
    AMaxErrors := (VectorIdxAbsMax(X1, 1, 5)<>5) or (VectorIdxAbsMax(X1, 1, 4)<>1) or (VectorIdxAbsMax(X1, 2, 4)<>3);
    N := 30;
    SetLength(X1, N+1);
    SetLength(A, N+1, N+1);
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            A[I,J] := 2*RandomReal-1;
            Inc(J);
        end;
        Inc(I);
    end;
    Was1 := False;
    Was2 := False;
    Pass:=1;
    while Pass<=1000 do
    begin
        J := 1+RandomInteger(N);
        I1 := 1+RandomInteger(N);
        I2 := I1+RandomInteger(N+1-I1);
        for i_ := I1 to I2 do
        begin
            X1[i_] := A[i_,J];
        end;
        if VectorIdxAbsMax(X1, I1, I2)<>ColumnIdxAbsMax(A, I1, I2, J) then
        begin
            Was1 := True;
        end;
        I := 1+RandomInteger(N);
        J1 := 1+RandomInteger(N);
        J2 := J1+RandomInteger(N+1-J1);
        APVMove(@X1[0], J1, J2, @A[I][0], J1, J2);
        if VectorIdxAbsMax(X1, J1, J2)<>RowIdxAbsMax(A, J1, J2, I) then
        begin
            Was2 := True;
        end;
        Inc(Pass);
    end;
    AMaxErrors := AMaxErrors or Was1 or Was2;
    
    //
    // Testing upper Hessenberg 1-norm
    //
    SetLength(A, 3+1, 3+1);
    SetLength(X1, 3+1);
    A[1,1] := 2;
    A[1,2] := 3;
    A[1,3] := 1;
    A[2,1] := 4;
    A[2,2] := -5;
    A[2,3] := 8;
    A[3,1] := 99;
    A[3,2] := 3;
    A[3,3] := 1;
    HSNErrors := AbsReal(UpperHessenberg1Norm(A, 1, 3, 1, 3, X1)-11)>10*MachineEpsilon;
    
    //
    // Testing MatrixVectorMultiply
    //
    SetLength(A, 3+1, 5+1);
    SetLength(X1, 3+1);
    SetLength(X2, 2+1);
    A[2,3] := 2;
    A[2,4] := -1;
    A[2,5] := -1;
    A[3,3] := 1;
    A[3,4] := -2;
    A[3,5] := 2;
    X1[1] := 1;
    X1[2] := 2;
    X1[3] := 1;
    X2[1] := -1;
    X2[2] := -1;
    MatrixVectorMultiply(A, 2, 3, 3, 5, False, X1, 1, 3, 1.0, X2, 1, 2, 1.0);
    MatrixVectorMultiply(A, 2, 3, 3, 5, True, X2, 1, 2, 1.0, X1, 1, 3, 1.0);
    E1 := AbsReal(X1[1]+5)+AbsReal(X1[2]-8)+AbsReal(X1[3]+1)+AbsReal(X2[1]+2)+AbsReal(X2[2]+2);
    X1[1] := 1;
    X1[2] := 2;
    X1[3] := 1;
    X2[1] := -1;
    X2[2] := -1;
    MatrixVectorMultiply(A, 2, 3, 3, 5, False, X1, 1, 3, 1.0, X2, 1, 2, 0.0);
    MatrixVectorMultiply(A, 2, 3, 3, 5, True, X2, 1, 2, 1.0, X1, 1, 3, 0.0);
    E2 := AbsReal(X1[1]+3)+AbsReal(X1[2]-3)+AbsReal(X1[3]+1)+AbsReal(X2[1]+1)+AbsReal(X2[2]+1);
    MVErrors := E1+E2>=50*MachineEpsilon;
    
    //
    // testing inplace transpose
    //
    N := 10;
    SetLength(A, N+1, N+1);
    SetLength(B, N+1, N+1);
    SetLength(X1, N-1+1);
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            A[I,J] := RandomReal;
            Inc(J);
        end;
        Inc(I);
    end;
    PassCount := 10000;
    Was1 := False;
    Pass:=1;
    while Pass<=PassCount do
    begin
        I1 := 1+RandomInteger(N);
        I2 := I1+RandomInteger(N-I1+1);
        J1 := 1+RandomInteger(N-(I2-I1));
        J2 := J1+(I2-I1);
        CopyMatrix(A, I1, I2, J1, J2, B, I1, I2, J1, J2);
        InplaceTranspose(B, I1, I2, J1, J2, X1);
        I:=I1;
        while I<=I2 do
        begin
            J:=J1;
            while J<=J2 do
            begin
                if A[I,J]<>B[I1+(J-J1),J1+(I-I1)] then
                begin
                    Was1 := True;
                end;
                Inc(J);
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;
    ITErrors := Was1;
    
    //
    // testing copy and transpose
    //
    N := 10;
    SetLength(A, N+1, N+1);
    SetLength(B, N+1, N+1);
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            A[I,J] := RandomReal;
            Inc(J);
        end;
        Inc(I);
    end;
    PassCount := 10000;
    Was1 := False;
    Pass:=1;
    while Pass<=PassCount do
    begin
        I1 := 1+RandomInteger(N);
        I2 := I1+RandomInteger(N-I1+1);
        J1 := 1+RandomInteger(N);
        J2 := J1+RandomInteger(N-J1+1);
        CopyAndTranspose(A, I1, I2, J1, J2, B, J1, J2, I1, I2);
        I:=I1;
        while I<=I2 do
        begin
            J:=J1;
            while J<=J2 do
            begin
                if A[I,J]<>B[J,I] then
                begin
                    Was1 := True;
                end;
                Inc(J);
            end;
            Inc(I);
        end;
        Inc(Pass);
    end;