testinverseupdateunit.pas

maths lib with source

unit testinverseupdateunit;
interface
uses Math, Ap, Sysutils, inverseupdate;

function TestInverseUpdate(Silent : Boolean):Boolean;
function testinverseupdateunit_test_silent():Boolean;
function testinverseupdateunit_test():Boolean;

implementation

procedure MakeACopy(const A : TReal2DArray;
     M : Integer;
     N : Integer;
     var B : TReal2DArray);forward;
procedure MatLU(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     var Pivots : TInteger1DArray);forward;
procedure GenerateRandomOrthogonalMatrix(var A0 : TReal2DArray;
     N : Integer);forward;
procedure GenerateRandomMatrixCond(var A0 : TReal2DArray;
     N : Integer;
     C : Double);forward;
function InvMatTR(var A : TReal2DArray;
     N : Integer;
     IsUpper : Boolean;
     IsUnitTriangular : Boolean):Boolean;forward;
function InvMatLU(var A : TReal2DArray;
     const Pivots : TInteger1DArray;
     N : Integer):Boolean;forward;
function InvMat(var A : TReal2DArray; N : Integer):Boolean;forward;
function MatrixDiff(const A : TReal2DArray;
     const B : TReal2DArray;
     M : Integer;
     N : Integer):Double;forward;
function UpdAndInv(var A : TReal2DArray;
     const U : TReal1DArray;
     const V : TReal1DArray;
     N : Integer):Boolean;forward;


function TestInverseUpdate(Silent : Boolean):Boolean;
var
    A : TReal2DArray;
    InvA : TReal2DArray;
    B1 : TReal2DArray;
    B2 : TReal2DArray;
    U : TReal1DArray;
    V : TReal1DArray;
    N : Integer;
    MaxN : Integer;
    I : Integer;
    J : Integer;
    UpdRow : Integer;
    UpdCol : Integer;
    Val : Double;
    Pass : Integer;
    PassCount : Integer;
    WasErrors : Boolean;
    Threshold : Double;
    C : Double;
begin
    WasErrors := False;
    MaxN := 10;
    PassCount := 100;
    Threshold := 1.0E+7*MachineEpsilon;
    
    //
    // process
    //
    N:=1;
    while N<=MaxN do
    begin
        SetLength(A, N-1+1, N-1+1);
        SetLength(B1, N-1+1, N-1+1);
        SetLength(B2, N-1+1, N-1+1);
        SetLength(U, N-1+1);
        SetLength(V, N-1+1);
        Pass:=1;
        while Pass<=PassCount do
        begin
            C := Exp(RandomReal*Ln(10));
            GenerateRandomMatrixCond(A, N, C);
            MakeACopy(A, N, N, InvA);
            if  not InvMat(InvA, N) then
            begin
                WasErrors := True;
                Break;
            end;
            
            //
            // Test simple update
            //
            UpdRow := RandomInteger(N);
            UpdCol := RandomInteger(N);
            Val := 0.1*(2*RandomReal-1);
            I:=0;
            while I<=N-1 do
            begin
                if I=UpdRow then
                begin
                    U[I] := Val;
                end
                else
                begin
                    U[I] := 0;
                end;
                if I=UpdCol then
                begin
                    V[I] := 1;
                end
                else
                begin
                    V[I] := 0;
                end;
                Inc(I);
            end;
            MakeACopy(A, N, N, B1);
            if  not UpdAndInv(B1, U, V, N) then
            begin
                WasErrors := True;
                Break;
            end;
            MakeACopy(InvA, N, N, B2);
            RMatrixInvUpdateSimple(B2, N, UpdRow, UpdCol, Val);
            WasErrors := WasErrors or (MatrixDiff(B1, B2, N, N)>Threshold);
            
            //
            // Test row update
            //
            UpdRow := RandomInteger(N);
            I:=0;
            while I<=N-1 do
            begin
                if I=UpdRow then
                begin
                    U[I] := 1;
                end
                else
                begin
                    U[I] := 0;
                end;
                V[I] := 0.1*(2*RandomReal-1);
                Inc(I);
            end;
            MakeACopy(A, N, N, B1);
            if  not UpdAndInv(B1, U, V, N) then
            begin
                WasErrors := True;
                Break;
            end;
            MakeACopy(InvA, N, N, B2);
            RMatrixInvUpdateRow(B2, N, UpdRow, V);
            WasErrors := WasErrors or (MatrixDiff(B1, B2, N, N)>Threshold);
            
            //
            // Test column update
            //
            UpdCol := RandomInteger(N);
            I:=0;
            while I<=N-1 do
            begin
                if I=UpdCol then
                begin
                    V[I] := 1;
                end
                else
                begin
                    V[I] := 0;
                end;
                U[I] := 0.1*(2*RandomReal-1);
                Inc(I);
            end;
            MakeACopy(A, N, N, B1);
            if  not UpdAndInv(B1, U, V, N) then
            begin
                WasErrors := True;
                Break;
            end;
            MakeACopy(InvA, N, N, B2);
            RMatrixInvUpdateColumn(B2, N, UpdCol, U);
            WasErrors := WasErrors or (MatrixDiff(B1, B2, N, N)>Threshold);
            
            //
            // Test full update
            //
            I:=0;
            while I<=N-1 do
            begin
                V[I] := 0.1*(2*RandomReal-1);
                U[I] := 0.1*(2*RandomReal-1);
                Inc(I);
            end;
            MakeACopy(A, N, N, B1);
            if  not UpdAndInv(B1, U, V, N) then
            begin
                WasErrors := True;
                Break;
            end;
            MakeACopy(InvA, N, N, B2);
            RMatrixInvUpdateUV(B2, N, U, V);
            WasErrors := WasErrors or (MatrixDiff(B1, B2, N, N)>Threshold);
            Inc(Pass);
        end;
        Inc(N);
    end;
    
    //
    // report
    //
    if  not Silent then
    begin
        Write(Format('TESTING INVERSE UPDATE (REAL)'#13#10'',[]));
        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;


(*************************************************************************
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;


(*************************************************************************
LU decomposition
*************************************************************************)
procedure MatLU(var A : TReal2DArray;
     M : Integer;
     N : Integer;
     var Pivots : TInteger1DArray);
var
    I : Integer;
    J : Integer;
    JP : Integer;
    T1 : TReal1DArray;
    s : Double;
    i_ : Integer;
begin
    SetLength(Pivots, Min(M-1, N-1)+1);
    SetLength(T1, Max(M-1, N-1)+1);
    Assert((M>=0) and (N>=0), 'Error in LUDecomposition: incorrect function arguments');
    
    //
    // Quick return if possible
    //
    if (M=0) or (N=0) then
    begin
        Exit;
    end;
    J:=0;
    while J<=Min(M-1, N-1) do
    begin
        
        //
        // Find pivot and test for singularity.
        //
        JP := J;
        I:=J+1;
        while I<=M-1 do
        begin
            if AbsReal(A[I,J])>AbsReal(A[JP,J]) then
            begin
                JP := I;
            end;
            Inc(I);
        end;
        Pivots[J] := JP;
        if A[JP,J]<>0 then
        begin
            
            //
            //Apply the interchange to rows
            //
            if JP<>J then
            begin
                APVMove(@T1[0], 0, N-1, @A[J][0], 0, N-1);
                APVMove(@A[J][0], 0, N-1, @A[JP][0], 0, N-1);
                APVMove(@A[JP][0], 0, N-1, @T1[0], 0, N-1);
            end;
            
            //
            //Compute elements J+1:M of J-th column.
            //
            if J<M then
            begin
                JP := J+1;
                S := 1/A[J,J];
                for i_ := JP to M-1 do
                begin
                    A[i_,J] := S*A[i_,J];
                end;
            end;
        end;
        if J<Min(M, N)-1 then
        begin
            
            //
            //Update trailing submatrix.
            //
            JP := J+1;
            I:=J+1;
            while I<=M-1 do
            begin
                S := A[I,J];
                APVSub(@A[I][0], JP, N-1, @A[J][0], JP, N-1, S);
                Inc(I);
            end;
        end;
        Inc(J);
    end;
end;


(*************************************************************************
Generate matrix with given condition number C (2-norm)
*************************************************************************)
procedure GenerateRandomOrthogonalMatrix(var A0 : TReal2DArray; N : Integer);
var
    T : Double;
    Lambda : Double;
    S : Integer;
    I : Integer;
    J : Integer;
    U1 : Double;
    U2 : Double;
    W : TReal1DArray;
    V : TReal1DArray;
    A : TReal2DArray;
    SM : Double;
begin
    if N<=0 then
    begin
        Exit;
    end;
    SetLength(W, N+1);
    SetLength(V, N+1);
    SetLength(A, N+1, N+1);
    SetLength(A0, N-1+1, N-1+1);
    
    //
    // Prepare A
    //
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            if I=J then
            begin
                A[I,J] := 1;
            end
            else
            begin
                A[I,J] := 0;
            end;
            Inc(J);
        end;
        Inc(I);
    end;
    
    //
    // Calculate A using Stewart algorithm
    //
    S:=2;
    while S<=N do
    begin
        
        //
        // Prepare v and Lambda = v'*v
        //
        repeat
            I := 1;
            while i<=S do
            begin
                U1 := 2*RandomReal-1;
                U2 := 2*RandomReal-1;
                SM := U1*U1+U2*U2;
                if (SM=0) or (SM>1) then
                begin
                    Continue;
                end;
                SM := Sqrt(-2*Ln(SM)/SM);
                V[I] := U1*SM;
                if I+1<=S then
                begin
                    V[I+1] := U2*SM;
                end;
                I := I+2;
            end;