testtdevdbiunit.pas

maths lib with source

    begin
        E[I] := 2*RandomReal-1;
        Inc(I);
    end;
end;


(*************************************************************************
Tests Z*Lambda*Z' against tridiag(D,E).
Returns relative error.
*************************************************************************)
function TestProduct(const D : TReal1DArray;
     const E : TReal1DArray;
     N : Integer;
     const Z : TReal2DArray;
     const Lambda : TReal1DArray):Double;
var
    I : Integer;
    J : Integer;
    K : Integer;
    V : Double;
    MX : Double;
begin
    Result := 0;
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            
            //
            // Calculate V = A[i,j], A = Z*Lambda*Z'
            //
            V := 0;
            K:=0;
            while K<=N-1 do
            begin
                V := V+Z[I,K]*Lambda[K]*Z[J,K];
                Inc(K);
            end;
            
            //
            // Compare
            //
            if AbsInt(I-J)=0 then
            begin
                Result := Max(Result, AbsReal(V-D[I]));
            end;
            if AbsInt(I-J)=1 then
            begin
                Result := Max(Result, AbsReal(V-E[Min(I, J)]));
            end;
            if AbsInt(I-J)>1 then
            begin
                Result := Max(Result, AbsReal(V));
            end;
            Inc(J);
        end;
        Inc(I);
    end;
    MX := 0;
    I:=0;
    while I<=N-1 do
    begin
        MX := Max(MX, AbsReal(D[I]));
        Inc(I);
    end;
    I:=0;
    while I<=N-2 do
    begin
        MX := Max(MX, AbsReal(E[I]));
        Inc(I);
    end;
    if MX=0 then
    begin
        MX := 1;
    end;
    Result := Result/MX;
end;


(*************************************************************************
Tests Z*Z' against diag(1...1)
Returns absolute error.
*************************************************************************)
function TestOrt(const Z : TReal2DArray; N : Integer):Double;
var
    I : Integer;
    J : Integer;
    V : Double;
    i_ : Integer;
begin
    Result := 0;
    I:=0;
    while I<=N-1 do
    begin
        J:=0;
        while J<=N-1 do
        begin
            V := 0.0;
            for i_ := 0 to N-1 do
            begin
                V := V + Z[i_,I]*Z[i_,J];
            end;
            if I=J then
            begin
                V := V-1;
            end;
            Result := Max(Result, AbsReal(V));
            Inc(J);
        end;
        Inc(I);
    end;
end;


(*************************************************************************
Tests EVD problem
*************************************************************************)
procedure TestEVDProblem(const D : TReal1DArray;
     const E : TReal1DArray;
     N : Integer;
     var ValErr : Double;
     var VecErr : Double;
     var WNSorted : Boolean;
     var FailC : Integer);
var
    Lambda : TReal1DArray;
    LambdaRef : TReal1DArray;
    Z : TReal2DArray;
    ZRef : TReal2DArray;
    A1 : TReal2DArray;
    A2 : TReal2DArray;
    AR : TReal2DArray;
    WSucc : Boolean;
    I : Integer;
    J : Integer;
    K : Integer;
    M : Integer;
    I1 : Integer;
    I2 : Integer;
    V : Double;
    A : Double;
    B : Double;
    i_ : Integer;
begin
    SetLength(LambdaRef, N-1+1);
    SetLength(ZRef, N-1+1, N-1+1);
    SetLength(A1, N-1+1, N-1+1);
    SetLength(A2, N-1+1, N-1+1);
    
    //
    // Reference EVD
    //
    if  not RefEVD(D, E, N, LambdaRef, ZRef) then
    begin
        FailC := FailC+1;
        Exit;
    end;
    
    //
    // Test different combinations
    //
    I1:=0;
    while I1<=N-1 do
    begin
        I2:=I1;
        while I2<=N-1 do
        begin
            
            //
            // Select A, B
            //
            if I1>0 then
            begin
                A := 0.5*(LambdaRef[I1]+LambdaRef[I1-1]);
            end
            else
            begin
                A := LambdaRef[0]-1;
            end;
            if I2<N-1 then
            begin
                B := 0.5*(LambdaRef[I2]+LambdaRef[I2+1]);
            end
            else
            begin
                B := LambdaRef[N-1]+1;
            end;
            
            //
            // Test interval, no vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            if  not SMatrixTDEVDR(Lambda, E, N, 0, A, B, M, Z) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            if M<>I2-I1+1 then
            begin
                FailC := FailC+1;
                Exit;
            end;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            
            //
            // Test indexes, no vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            if  not SMatrixTDEVDI(Lambda, E, N, 0, I1, I2, Z) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            M := I2-I1+1;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            
            //
            // Test interval, transform vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            SetLength(A1, N-1+1, N-1+1);
            SetLength(A2, N-1+1, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A1[I,J] := 2*RandomReal-1;
                    A2[I,J] := A1[I,J];
                    Inc(J);
                end;
                Inc(I);
            end;
            if  not SMatrixTDEVDR(Lambda, E, N, 1, A, B, M, A1) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            if M<>I2-I1+1 then
            begin
                FailC := FailC+1;
                Exit;
            end;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            SetLength(AR, N-1+1, M-1+1);
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    V := 0.0;
                    for i_ := 0 to N-1 do
                    begin
                        V := V + A2[I,i_]*ZRef[i_,I1+J];
                    end;
                    AR[I,J] := V;
                    Inc(J);
                end;
                Inc(I);
            end;
            J:=0;
            while J<=M-1 do
            begin
                V := 0.0;
                for i_ := 0 to N-1 do
                begin
                    V := V + A1[i_,J]*AR[i_,J];
                end;
                if V<0 then
                begin
                    for i_ := 0 to N-1 do
                    begin
                        AR[i_,J] := -1*AR[i_,J];
                    end;
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsReal(A1[I,J]-AR[I,J]));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test indexes, transform vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            SetLength(A1, N-1+1, N-1+1);
            SetLength(A2, N-1+1, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=N-1 do
                begin
                    A1[I,J] := 2*RandomReal-1;
                    A2[I,J] := A1[I,J];
                    Inc(J);
                end;
                Inc(I);
            end;
            if  not SMatrixTDEVDI(Lambda, E, N, 1, I1, I2, A1) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            M := I2-I1+1;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            SetLength(AR, N-1+1, M-1+1);
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    V := 0.0;
                    for i_ := 0 to N-1 do
                    begin
                        V := V + A2[I,i_]*ZRef[i_,I1+J];
                    end;
                    AR[I,J] := V;
                    Inc(J);
                end;
                Inc(I);
            end;
            J:=0;
            while J<=M-1 do
            begin
                V := 0.0;
                for i_ := 0 to N-1 do
                begin
                    V := V + A1[i_,J]*AR[i_,J];
                end;
                if V<0 then
                begin
                    for i_ := 0 to N-1 do
                    begin
                        AR[i_,J] := -1*AR[i_,J];
                    end;
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsReal(A1[I,J]-AR[I,J]));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test interval, do not transform vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            SetLength(Z, 0+1, 0+1);
            if  not SMatrixTDEVDR(Lambda, E, N, 2, A, B, M, Z) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            if M<>I2-I1+1 then
            begin
                FailC := FailC+1;
                Exit;
            end;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            J:=0;
            while J<=M-1 do
            begin
                V := 0.0;
                for i_ := 0 to N-1 do
                begin
                    V := V + Z[i_,J]*ZRef[i_,I1+J];
                end;
                if V<0 then
                begin
                    for i_ := 0 to N-1 do
                    begin
                        Z[i_,J] := -1*Z[i_,J];
                    end;
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsReal(Z[I,J]-ZRef[I,I1+J]));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test interval, do not transform vectors
            //
            SetLength(Lambda, N-1+1);
            I:=0;
            while I<=N-1 do
            begin
                Lambda[I] := D[I];
                Inc(I);
            end;
            SetLength(Z, 0+1, 0+1);
            if  not SMatrixTDEVDI(Lambda, E, N, 2, I1, I2, Z) then
            begin
                FailC := FailC+1;
                Exit;
            end;
            M := I2-I1+1;
            K:=0;
            while K<=M-1 do
            begin
                ValErr := Max(ValErr, AbsReal(Lambda[K]-LambdaRef[I1+K]));
                Inc(K);
            end;
            J:=0;
            while J<=M-1 do
            begin
                V := 0.0;
                for i_ := 0 to N-1 do
                begin
                    V := V + Z[i_,J]*ZRef[i_,I1+J];
                end;
                if V<0 then
                begin
                    for i_ := 0 to N-1 do
                    begin
                        Z[i_,J] := -1*Z[i_,J];
                    end;
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsReal(Z[I,J]-ZRef[I,I1+J]));
                    Inc(J);
                end;
                Inc(I);
            end;
            Inc(I2);
        end;
        Inc(I1);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testtdevdbiunit_test_silent():Boolean;
begin
    Result := TestTDEVDBI(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testtdevdbiunit_test():Boolean;
begin
    Result := TestTDEVDBI(False);
end;


end.