testinverseupdateunit.pas

maths lib with source

            Lambda := APVDotProduct(@V[0], 1, S, @V[0], 1, S);
        until Lambda<>0;
        Lambda := 2/Lambda;
        
        //
        // A * (I - 2 vv'/v'v ) =
        //   = A - (2/v'v) * A * v * v' =
        //   = A - (2/v'v) * w * v'
        //  where w = Av
        //
        I:=1;
        while I<=S do
        begin
            T := APVDotProduct(@A[I][0], 1, S, @V[0], 1, S);
            W[I] := T;
            Inc(I);
        end;
        I:=1;
        while I<=S do
        begin
            T := W[I]*Lambda;
            APVSub(@A[I][0], 1, S, @V[0], 1, S, T);
            Inc(I);
        end;
        Inc(S);
    end;
    
    //
    //
    //
    I:=1;
    while I<=N do
    begin
        J:=1;
        while J<=N do
        begin
            A0[I-1,J-1] := A[I,J];
            Inc(J);
        end;
        Inc(I);
    end;
end;


procedure GenerateRandomMatrixCond(var A0 : TReal2DArray;
     N : Integer;
     C : Double);
var
    L1 : Double;
    L2 : Double;
    Q1 : TReal2DArray;
    Q2 : TReal2DArray;
    CC : TReal1DArray;
    I : Integer;
    J : Integer;
    K : Integer;
begin
    GenerateRandomOrthogonalMatrix(Q1, N);
    GenerateRandomOrthogonalMatrix(Q2, N);
    SetLength(CC, N-1+1);
    L1 := 0;
    L2 := Ln(1/C);
    CC[0] := Exp(L1);
    I:=1;
    while I<=N-2 do
    begin
        CC[I] := Exp(RandomReal*(L2-L1)+L1);
        Inc(I);
    end;
    CC[N-1] := Exp(L2);
    SetLength(A0, N-1+1, N-1+1);
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            A0[I,J] := 0;
            K:=0;
            while K<=N-1 do
            begin
                A0[I,J] := A0[I,J]+Q1[I,K]*CC[K]*Q2[J,K];
                Inc(K);
            end;
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
triangular inverse
*************************************************************************)
function InvMatTR(var A : TReal2DArray;
     N : Integer;
     IsUpper : Boolean;
     IsUnitTriangular : Boolean):Boolean;
var
    NOUNIT : Boolean;
    I : Integer;
    J : Integer;
    V : Double;
    AJJ : Double;
    T : TReal1DArray;
    i_ : Integer;
begin
    Result := True;
    SetLength(T, N-1+1);
    
    //
    // Test the input parameters.
    //
    NOUNIT :=  not IsUnitTriangular;
    if IsUpper then
    begin
        
        //
        // Compute inverse of upper triangular matrix.
        //
        J:=0;
        while J<=N-1 do
        begin
            if NOUNIT then
            begin
                if A[J,J]=0 then
                begin
                    Result := False;
                    Exit;
                end;
                A[J,J] := 1/A[J,J];
                AJJ := -A[J,J];
            end
            else
            begin
                AJJ := -1;
            end;
            
            //
            // Compute elements 1:j-1 of j-th column.
            //
            if J>0 then
            begin
                for i_ := 0 to J-1 do
                begin
                    T[i_] := A[i_,J];
                end;
                I:=0;
                while I<=J-1 do
                begin
                    if I<J-1 then
                    begin
                        V := APVDotProduct(@A[I][0], I+1, J-1, @T[0], I+1, J-1);
                    end
                    else
                    begin
                        V := 0;
                    end;
                    if NOUNIT then
                    begin
                        A[I,J] := V+A[I,I]*T[I];
                    end
                    else
                    begin
                        A[I,J] := V+T[I];
                    end;
                    Inc(I);
                end;
                for i_ := 0 to J-1 do
                begin
                    A[i_,J] := AJJ*A[i_,J];
                end;
            end;
            Inc(J);
        end;
    end
    else
    begin
        
        //
        // Compute inverse of lower triangular matrix.
        //
        J:=N-1;
        while J>=0 do
        begin
            if NOUNIT then
            begin
                if A[J,J]=0 then
                begin
                    Result := False;
                    Exit;
                end;
                A[J,J] := 1/A[J,J];
                AJJ := -A[J,J];
            end
            else
            begin
                AJJ := -1;
            end;
            if J<N-1 then
            begin
                
                //
                // Compute elements j+1:n of j-th column.
                //
                for i_ := J+1 to N-1 do
                begin
                    T[i_] := A[i_,J];
                end;
                I:=J+1;
                while I<=N-1 do
                begin
                    if I>J+1 then
                    begin
                        V := APVDotProduct(@A[I][0], J+1, I-1, @T[0], J+1, I-1);
                    end
                    else
                    begin
                        V := 0;
                    end;
                    if NOUNIT then
                    begin
                        A[I,J] := V+A[I,I]*T[I];
                    end
                    else
                    begin
                        A[I,J] := V+T[I];
                    end;
                    Inc(I);
                end;
                for i_ := J+1 to N-1 do
                begin
                    A[i_,J] := AJJ*A[i_,J];
                end;
            end;
            Dec(J);
        end;
    end;
end;


(*************************************************************************
LU inverse
*************************************************************************)
function InvMatLU(var A : TReal2DArray;
     const Pivots : TInteger1DArray;
     N : Integer):Boolean;
var
    WORK : TReal1DArray;
    I : Integer;
    IWS : Integer;
    J : Integer;
    JB : Integer;
    JJ : Integer;
    JP : Integer;
    V : Double;
    i_ : Integer;
begin
    Result := True;
    
    //
    // Quick return if possible
    //
    if N=0 then
    begin
        Exit;
    end;
    SetLength(WORK, N-1+1);
    
    //
    // Form inv(U)
    //
    if  not InvMatTR(A, N, True, False) then
    begin
        Result := False;
        Exit;
    end;
    
    //
    // Solve the equation inv(A)*L = inv(U) for inv(A).
    //
    J:=N-1;
    while J>=0 do
    begin
        
        //
        // Copy current column of L to WORK and replace with zeros.
        //
        I:=J+1;
        while I<=N-1 do
        begin
            WORK[I] := A[I,J];
            A[I,J] := 0;
            Inc(I);
        end;
        
        //
        // Compute current column of inv(A).
        //
        if J<N-1 then
        begin
            I:=0;
            while I<=N-1 do
            begin
                V := APVDotProduct(@A[I][0], J+1, N-1, @WORK[0], J+1, N-1);
                A[I,J] := A[I,J]-V;
                Inc(I);
            end;
        end;
        Dec(J);
    end;
    
    //
    // Apply column interchanges.
    //
    J:=N-2;
    while J>=0 do
    begin
        JP := Pivots[J];
        if JP<>J then
        begin
            for i_ := 0 to N-1 do
            begin
                WORK[i_] := A[i_,J];
            end;
            for i_ := 0 to N-1 do
            begin
                A[i_,J] := A[i_,JP];
            end;
            for i_ := 0 to N-1 do
            begin
                A[i_,JP] := WORK[i_];
            end;
        end;
        Dec(J);
    end;
end;


(*************************************************************************
Matrix inverse
*************************************************************************)
function InvMat(var A : TReal2DArray; N : Integer):Boolean;
var
    Pivots : TInteger1DArray;
begin
    MatLU(A, N, N, Pivots);
    Result := InvMatLU(A, Pivots, N);
end;


(*************************************************************************
Diff
*************************************************************************)
function MatrixDiff(const A : TReal2DArray;
     const B : TReal2DArray;
     M : Integer;
     N : Integer):Double;
var
    I : Integer;
    J : Integer;
begin
    Result := 0;
    I:=0;
    while I<=M-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            Result := Max(Result, AbsReal(B[I,J]-A[I,J]));
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Update and inverse
*************************************************************************)
function UpdAndInv(var A : TReal2DArray;
     const U : TReal1DArray;
     const V : TReal1DArray;
     N : Integer):Boolean;
var
    Pivots : TInteger1DArray;
    I : Integer;
    R : Double;
begin
    I:=0;
    while I<=N-1 do
    begin
        R := U[I];
        APVAdd(@A[I][0], 0, N-1, @V[0], 0, N-1, R);
        Inc(I);
    end;
    MatLU(A, N, N, Pivots);
    Result := InvMatLU(A, Pivots, N);
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testinverseupdateunit_test_silent():Boolean;
begin
    Result := TestInverseUpdate(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testinverseupdateunit_test():Boolean;
begin
    Result := TestInverseUpdate(False);
end;


end.