testclinsolve.pas

maths lib with source

unit testclinsolve;
interface
uses Math, Ap, Sysutils, clu, csolve;

function TestCSolveSystem(Silent : Boolean):Boolean;
function testclinsolve_test_silent():Boolean;
function testclinsolve_test():Boolean;

implementation

function TestCSolveSystem(Silent : Boolean):Boolean;
var
    N : Integer;
    MaxN : Integer;
    I : Integer;
    J : Integer;
    PassCount : Integer;
    Pass : Integer;
    WasErrors : Boolean;
    WFailed : Boolean;
    Err : Double;
    V : Complex;
    Threshold : Double;
    A : TComplex2DArray;
    TX : TComplex1DArray;
    X : TComplex1DArray;
    B : TComplex1DArray;
    P : TInteger1DArray;
    i_ : Integer;
begin
    Err := 0;
    WFailed := False;
    WasErrors := False;
    MaxN := 30;
    PassCount := 10;
    Threshold := 2*MaxN*1000*MachineEpsilon;
    
    //
    // 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].X := 2*RandomReal-1;
                    A[I,J].Y := 2*RandomReal-1;
                    Inc(J);
                end;
                Inc(I);
            end;
            A[RandomInteger(N),RandomInteger(N)] := C_Complex(10);
            A[RandomInteger(N),RandomInteger(N)] := C_Complex(7);
            I:=0;
            while I<=N-1 do
            begin
                TX[I].X := 2*RandomReal-1;
                TX[I].Y := 2*RandomReal-1;
                Inc(I);
            end;
            I:=0;
            while I<=N-1 do
            begin
                V := C_Complex(0.0);
                for i_ := 0 to N-1 do
                begin
                    V := C_Add(V,C_Mul(A[I,i_],TX[i_]));
                end;
                B[I] := V;
                Inc(I);
            end;
            
            //
            // solve and test normal
            //
            if  not CMatrixSolve(A, B, N, X) then
            begin
                WFailed := True;
                Inc(Pass);
                Continue;
            end;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsComplex(C_Sub(TX[I],X[I])));
                Inc(I);
            end;
            
            //
            // solve and test LU
            //
            CMatrixLU(A, N, N, P);
            if  not CMatrixLUSolve(A, P, B, N, X) then
            begin
                WFailed := True;
                Inc(Pass);
                Continue;
            end;
            I:=0;
            while I<=N-1 do
            begin
                Err := Max(Err, AbsComplex(C_Sub(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 COMPLEX 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 testclinsolve_test_silent():Boolean;
begin
    Result := TestCSolveSystem(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testclinsolve_test():Boolean;
begin
    Result := TestCSolveSystem(False);
end;


end.