test_Matrix3.cpp

Go to the documentation of this file.

#include <vector>
#include <cassert>
#include "moab/Matrix3.hpp"
#include "moab/Core.hpp"
#include "moab/CartVect.hpp"
#include "TestUtil.hpp"

using namespace moab;

#define ACOS( x ) acos( std::min( std::max( x, -1.0 ), 1.0 ) )

double find_angle( const moab::CartVect& A, const moab::CartVect& B )
{
    const double lA = A.length(), lB = B.length();
    assert( lA > 0.0 );
    assert( lB > 0.0 );
    const double dPI = 3.14159265;
    const double dot = ( A[0] * B[0] + A[1] * B[1] + A[2] * B[2] );
    return ACOS( dot / ( lA * lB ) ) * 180.0 / dPI;
}

#define CHECK_EIGVECREAL_EQUAL( EXP, ACT, EPS ) \
    check_equal_eigvect( ( EXP ), ( ACT ), ( EPS ), #EXP, #ACT, __LINE__, __FILE__ )
void check_equal_eigvect( const moab::CartVect& A,
                          const moab::CartVect& B,
                          double eps,
                          const char* sA,
                          const char* sB,
                          int line,
                          const char* file )
{
    check_equal( A.length(), B.length(), eps, sA, sB, line, file );

    double angle = find_angle( A, B );

    if( ( fabs( A[0] - B[0] ) <= eps || fabs( A[0] + B[0] ) <= eps ) &&
        ( fabs( A[1] - B[1] ) <= eps || fabs( A[1] + B[1] ) <= eps ) &&
        ( fabs( A[2] - B[2] ) <= eps || fabs( A[2] + B[2] ) <= eps ) &&
        ( angle <= eps || fabs( angle - 180.0 ) <= eps ) )
        return;

    std::cout << "Equality Test Failed: " << sA << " == " << sB << std::endl;
    std::cout << "  at line " << line << " of '" << file << "'" << std::endl;

    std::cout << "  Expected: ";
    std::cout << A << std::endl;

    std::cout << "  Actual:   ";
    std::cout << B << std::endl;

    std::cout << "Angle between vectors := " << angle << std::endl;

    flag_error();
}

void test_EigenDecomp();

int main()
{

    int result = 0;

    result += RUN_TEST( test_EigenDecomp );

    return result;
}

// test to ensure the Eigenvalues/vectors are calculated correctly and returned properly
// from the Matrix3 class for a simple case
void test_EigenDecomp()
{
    // Create a matrix
    moab::Matrix3 mat;

    mat( 0 ) = 2;
    mat( 1 ) = -1;
    mat( 2 ) = 0;
    mat( 3 ) = -1;
    mat( 4 ) = 2;
    mat( 5 ) = -1;
    mat( 6 ) = 0;
    mat( 7 ) = -1;
    mat( 8 ) = 2;

    // now do the Eigen Decomposition of this Matrix

    CartVect lamda;
    Matrix3 vectors;
    moab::ErrorCode rval = mat.eigen_decomposition( lamda, vectors );CHECK_ERR( rval );
    for( int i = 0; i < 3; ++i )
        vectors.col( i ).normalize();

    // Hardcoded check values for the results
    double lamda_check[3];
    lamda_check[0] = 0.585786;
    lamda_check[1] = 2.0;
    lamda_check[2] = 3.41421;

    moab::CartVect vec0_check( 0.5, 0.707107, 0.5 );
    moab::CartVect vec1_check( 0.707107, 3.37748e-17, -0.707107 );
    moab::CartVect vec2_check( 0.5, -0.707107, 0.5 );

    // now verfy that the returns Eigenvalues and Eigenvectors are correct (within some tolerance)
    double tol = 1e-04;

    // use a deprecated constructor
    Matrix3 mat3( CartVect( 2, -1, 0 ), CartVect( -1, 2, -1 ), CartVect( 0, -1, 2 ), true );
    CHECK_REAL_EQUAL( mat( 1 ), mat3( 1 ), tol );

    vec0_check.normalize();
    vec1_check.normalize();
    vec2_check.normalize();

    // check that the correct Eigenvalues are returned correctly (in order)
    CHECK_REAL_EQUAL( lamda[0], lamda_check[0], tol );
    CHECK_REAL_EQUAL( lamda[1], lamda_check[1], tol );
    CHECK_REAL_EQUAL( lamda[2], lamda_check[2], tol );

    // check the Eigenvector values (order should correspond to the Eigenvalues)
    // first vector
    CHECK_EIGVECREAL_EQUAL( vectors.col( 0 ), vec0_check, tol );

    // sceond vector
    CHECK_EIGVECREAL_EQUAL( vectors.col( 1 ), vec1_check, tol );

    // third vector
    CHECK_EIGVECREAL_EQUAL( vectors.col( 2 ), vec2_check, tol );

    // another check to ensure the result is valid (AM-kM = 0)
    for( unsigned i = 0; i < 3; ++i )
    {
        moab::CartVect v = moab::Matrix::matrix_vector( mat, vectors.col( i ) ) - lamda[i] * vectors.col( i );
        CHECK_REAL_EQUAL( v.length(), 0, tol );
    }

    // for a real, symmetric matrix the Eigenvectors should be orthogonal
    CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 1 ), 0, tol );
    CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 2 ), 0, tol );

    return;
}