testcrcondunit.pas

maths lib with source

                    if NOUNIT then
                    begin
                        A[I,J] := C_Add(V,C_Mul(A[I,I],T[I]));
                    end
                    else
                    begin
                        A[I,J] := C_Add(V,T[I]);
                    end;
                    Inc(I);
                end;
                for i_ := 0 to J-1 do
                begin
                    A[i_,J] := C_Mul(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 C_EqualR(A[J,J],0) then
                begin
                    Result := False;
                    Exit;
                end;
                A[J,J] := C_RDiv(1,A[J,J]);
                AJJ := C_Opposite(A[J,J]);
            end
            else
            begin
                AJJ := C_Complex(-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 := C_Complex(0.0);
                        for i_ := J+1 to I-1 do
                        begin
                            V := C_Add(V,C_Mul(A[I,i_],T[i_]));
                        end;
                    end
                    else
                    begin
                        V := C_Complex(0);
                    end;
                    if NOUNIT then
                    begin
                        A[I,J] := C_Add(V,C_Mul(A[I,I],T[I]));
                    end
                    else
                    begin
                        A[I,J] := C_Add(V,T[I]);
                    end;
                    Inc(I);
                end;
                for i_ := J+1 to N-1 do
                begin
                    A[i_,J] := C_Mul(AJJ, A[i_,J]);
                end;
            end;
            Dec(J);
        end;
    end;
end;


(*************************************************************************
LU inverse
*************************************************************************)
function InvMatLU(var A : TComplex2DArray;
     const Pivots : TInteger1DArray;
     N : Integer):Boolean;
var
    WORK : TComplex1DArray;
    I : Integer;
    IWS : Integer;
    J : Integer;
    JB : Integer;
    JJ : Integer;
    JP : Integer;
    V : Complex;
    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] := C_Complex(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 := C_Complex(0.0);
                for i_ := J+1 to N-1 do
                begin
                    V := C_Add(V,C_Mul(A[I,i_],WORK[i_]));
                end;
                A[I,J] := C_Sub(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 : TComplex2DArray; N : Integer):Boolean;
var
    Pivots : TInteger1DArray;
begin
    CMatrixLU(A, N, N, Pivots);
    Result := InvMatLU(A, Pivots, N);
end;


(*************************************************************************
reference RCond
*************************************************************************)
procedure RefRCond(const A : TComplex2DArray;
     N : Integer;
     var RC1 : Double;
     var RCInf : Double);
var
    InvA : TComplex2DArray;
    Nrm1A : Double;
    NrmInfA : Double;
    Nrm1InvA : Double;
    NrmInfInvA : Double;
    V : Double;
    K : Integer;
    I : Integer;
begin
    
    //
    // inv A
    //
    MakeACopy(A, N, N, InvA);
    if  not InvMat(InvA, N) then
    begin
        RC1 := 0;
        RCInf := 0;
        Exit;
    end;
    
    //
    // norm A
    //
    Nrm1A := 0;
    NrmInfA := 0;
    K:=0;
    while K<=N-1 do
    begin
        V := 0;
        I:=0;
        while I<=N-1 do
        begin
            V := V+AbsComplex(A[I,K]);
            Inc(I);
        end;
        Nrm1A := Max(Nrm1A, V);
        V := 0;
        I:=0;
        while I<=N-1 do
        begin
            V := V+AbsComplex(A[K,I]);
            Inc(I);
        end;
        NrmInfA := Max(NrmInfA, V);
        Inc(K);
    end;
    
    //
    // norm inv A
    //
    Nrm1InvA := 0;
    NrmInfInvA := 0;
    K:=0;
    while K<=N-1 do
    begin
        V := 0;
        I:=0;
        while I<=N-1 do
        begin
            V := V+AbsComplex(InvA[I,K]);
            Inc(I);
        end;
        Nrm1InvA := Max(Nrm1InvA, V);
        V := 0;
        I:=0;
        while I<=N-1 do
        begin
            V := V+AbsComplex(InvA[K,I]);
            Inc(I);
        end;
        NrmInfInvA := Max(NrmInfInvA, V);
        Inc(K);
    end;
    
    //
    // result
    //
    RC1 := Nrm1InvA*Nrm1A;
    RCInf := NrmInfInvA*NrmInfA;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testcrcondunit_test_silent():Boolean;
begin
    Result := TestCRCond(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testcrcondunit_test():Boolean;
begin
    Result := TestCRCond(False);
end;


end.