testhevdbiunit.pas

maths lib with source

    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 AFull : TComplex2DArray;
     const AL : TComplex2DArray;
     const AU : TComplex2DArray;
     N : Integer;
     var ValErr : Double;
     var VecErr : Double;
     var WNSorted : Boolean;
     var FailC : Integer);
var
    Lambda : TReal1DArray;
    LambdaRef : TReal1DArray;
    Z : TComplex2DArray;
    ZRef : TComplex2DArray;
    WSucc : Boolean;
    I : Integer;
    J : Integer;
    K : Integer;
    M : Integer;
    I1 : Integer;
    I2 : Integer;
    V : Complex;
    A : Double;
    B : Double;
    i_ : Integer;
begin
    SetLength(LambdaRef, N-1+1);
    SetLength(ZRef, N-1+1, N-1+1);
    
    //
    // Reference EVD
    //
    if  not RefEVD(AFull, 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, lower A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDR(AL, N, 0, False, A, B, M, Lambda, 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 interval, no vectors, upper A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDR(AU, N, 0, True, A, B, M, Lambda, 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, lower A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDI(AL, N, 0, False, I1, I2, Lambda, 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 indexes, no vectors, upper A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDI(AU, N, 0, True, I1, I2, Lambda, 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, do not transform vectors, lower A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDR(AL, N, 1, False, A, B, M, Lambda, 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 := C_Complex(0.0);
                for i_ := 0 to N-1 do
                begin
                    V := C_Add(V,C_Mul(Z[i_,J],Conj(ZRef[i_,I1+J])));
                end;
                V := C_RDiv(1,V);
                for i_ := 0 to N-1 do
                begin
                    Z[i_,J] := C_Mul(V, Z[i_,J]);
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsComplex(C_Sub(Z[I,J],ZRef[I,I1+J])));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test interval, do not transform vectors, upper A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDR(AU, N, 1, True, A, B, M, Lambda, 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 := C_Complex(0.0);
                for i_ := 0 to N-1 do
                begin
                    V := C_Add(V,C_Mul(Z[i_,J],Conj(ZRef[i_,I1+J])));
                end;
                V := C_RDiv(1,V);
                for i_ := 0 to N-1 do
                begin
                    Z[i_,J] := C_Mul(V, Z[i_,J]);
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsComplex(C_Sub(Z[I,J],ZRef[I,I1+J])));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test indexes, do not transform vectors, lower A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDI(AL, N, 1, False, I1, I2, Lambda, 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 := C_Complex(0.0);
                for i_ := 0 to N-1 do
                begin
                    V := C_Add(V,C_Mul(Z[i_,J],Conj(ZRef[i_,I1+J])));
                end;
                V := C_RDiv(1,V);
                for i_ := 0 to N-1 do
                begin
                    Z[i_,J] := C_Mul(V, Z[i_,J]);
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsComplex(C_Sub(Z[I,J],ZRef[I,I1+J])));
                    Inc(J);
                end;
                Inc(I);
            end;
            
            //
            // Test indexes, do not transform vectors, upper A
            //
            Unset1D(Lambda);
            Unset2D(Z);
            if  not HMatrixEVDI(AU, N, 1, True, I1, I2, Lambda, 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 := C_Complex(0.0);
                for i_ := 0 to N-1 do
                begin
                    V := C_Add(V,C_Mul(Z[i_,J],Conj(ZRef[i_,I1+J])));
                end;
                V := C_RDiv(1,V);
                for i_ := 0 to N-1 do
                begin
                    Z[i_,J] := C_Mul(V, Z[i_,J]);
                end;
                Inc(J);
            end;
            I:=0;
            while I<=N-1 do
            begin
                J:=0;
                while J<=M-1 do
                begin
                    VecErr := Max(VecErr, AbsComplex(C_Sub(Z[I,J],ZRef[I,I1+J])));
                    Inc(J);
                end;
                Inc(I);
            end;
            Inc(I2);
        end;
        Inc(I1);
    end;
end;


(*************************************************************************
Silent unit test
*************************************************************************)
function testhevdbiunit_test_silent():Boolean;
begin
    Result := TestHEVDBI(True);
end;


(*************************************************************************
Unit test
*************************************************************************)
function testhevdbiunit_test():Boolean;
begin
    Result := TestHEVDBI(False);
end;


end.