testlinsolve.pas

maths lib with source

unit testlinsolve;
interface
uses Math, Ap, Sysutils, lu, rsolve;

function TestSolveSystem(Silent : Boolean):Boolean;
function testlinsolve_test_silent():Boolean;
function testlinsolve_test():Boolean;

implementation

function TestSolveSystem(Silent : Boolean):Boolean;
var
    N : Integer;
    MaxN : Integer;
    I : Integer;
    J : Integer;
    PassCount : Integer;
    Pass : Integer;
    WasErrors : Boolean;
    WFailed : Boolean;
    Err : Double;
    V : Double;
    Threshold : Double;
    A : TReal2DArray;
    TX : TReal1DArray;
    X : TReal1DArray;
    B : TReal1DArray;
    P : TInteger1DArray;
begin
    Err := 0;
    WFailed := False;
    WasErrors := False;
    MaxN := 30;
    PassCount := 10;
    Threshold := 1.0E-8;
    
    //
    // N = 1 .. 30
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(A, N-1+1, N-1+1);
        SetLength(TX, N-1+1);
        SetLength(B, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            
            //
            // init A, TX
            //
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A[I,J] := 2*RandomReal-1;
                    Inc(J);
                end;
                Inc(I);
            end;
            A[RandomInteger(N),RandomInteger(N)] := 10;
            A[RandomInteger(N),RandomInteger(N)] := 7;
            I:=0;
            while I<=N-1 do
            begin
                TX[I] := 2*RandomReal-1;
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                V := APVDotProduct(@A[I][0], 0, N-1, @TX[0], 0, N-1);
                B[I] := V;
                Inc(I);
            end;
            
            //
            // solve and test normal
            //
            if  not RMatrixSolve(A, B, N, X) then
            begin
                WFailed := True;
                Inc(Pass);
                Continue;
            end;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsReal(TX[I]-X[I]));
                Inc(I);
            end;
            
            //
            // solve and test LU
            //
            RMatrixLU(A, N, N, P);
            if  not RMatrixLUSolve(A, P, B, N, X) then
            begin
                WFailed := True;
                Inc(Pass);
                Continue;
            end;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsReal(TX[I]-X[I]));
                Inc(I);
            end;
            Inc(Pass);
        end;
        Inc(N);
    end;
    
    //
    // report
    //
    WasErrors := (Err>Threshold) or WFailed;
    if  not Silent then
    begin
        Write(Format('TESTING LINSOLVE'#13#10'',[]));
        Write(Format('Error:                                   %5.4e'#13#10'',[
            Err]));
        Write(Format('Always succeeded:                        ',[]));
        if  not WFailed then
        begin
            Write(Format('YES'#13#10'',[]));
        end
        else
        begin
            Write(Format('NO'#13#10'',[]));
        end;
        Write(Format('Threshold:                               %5.4e'#13#10'',[
            Threshold]));
        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;


(*************************************************************************
Silent unit test
*************************************************************************)
function testlinsolve_test_silent():Boolean;
begin
    Result := TestSolveSystem(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testlinsolve_test():Boolean;
begin
    Result := TestSolveSystem(False);
end;


end.