fn_AR_eigen.hpp

Go to the documentation of this file.

/*
 * BELFEM -- The Berkeley Lab Finite Element Framework
 * Copyright (c) 2026, The Regents of the University of California,
 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
 * approvals from the U.S. Dept. of Energy).  All rights reserved.
 *
 * Developers: Christian Messe, Gregory Giard
 *
 * See the top-level LICENSE file for the complete license and disclaimer.
 */
 
#ifndef BELFEM_FN_AR_EIGEN_HPP
#define BELFEM_FN_AR_EIGEN_HPP
#include "typedefs.hpp"
#include "assert.hpp"
#include "cl_AR_Vector.hpp"
#include "cl_AR_Matrix.hpp"
 
namespace belfem
{
//------------------------------------------------------------------------------
 
    inline int_t
    eigen( const Matrix< real > &  aMatrix,
           Vector< real > & aValues,
           const bool aAbortOnComplex = true )
    {
        BELFEM_ASSERT( aMatrix.n_cols() == aMatrix.n_rows(),
                      "Matrix must be quadratic");
 
        arma::cx_vec tValues ;
 
        arma::eig_gen( tValues, aMatrix.matrix_data()  );
 
        size_t tN = aMatrix.n_cols() ;
        aValues.set_size( aMatrix.n_cols() );
 
        int_t tNumComplex = 0 ;
 
        for( size_t k=0; k<tN; ++k )
        {
            if( std::abs( std::imag( tValues( k ) ) ) > BELFEM_EPSILON )
            {
                BELFEM_ERROR( ! aAbortOnComplex,
                    "eigen(): eigenvalue %u of this matrix is complex ( %g %+g i ), and a real "
                    "Vector cannot hold it. Use eigen_sym() if the matrix is symmetric, or pass "
                    "aAbortOnComplex = false to receive a count and NaN entries instead.",
                    ( unsigned int ) k,
                    ( double ) std::real( tValues( k ) ),
                    ( double ) std::imag( tValues( k ) ) );
 
                aValues( k ) = BELFEM_QUIET_NAN ;
                ++tNumComplex ;
            }
            else
            {
                aValues( k ) = std::real( tValues( k ) );
            }
        }
 
        return tNumComplex ;
    }
 
//------------------------------------------------------------------------------
 
    inline void
    eigen_sym( const Matrix< real > & aMatrix,
               Vector< real > & aValues )
    {
        BELFEM_ASSERT( aMatrix.n_cols() == aMatrix.n_rows(),
                      "Matrix must be quadratic");
 
        // arma::eig_sym writes into an arma::Col, but Vector< T > is backed by an
        // arma::Mat, so the eigenvalues are collected in a local and copied out
        arma::Col< real > tValues ;
 
        BELFEM_ERROR( arma::eig_sym( tValues, aMatrix.matrix_data() ),
                      "eigen_sym(): the symmetric eigenvalue decomposition did not converge" );
 
        size_t tN = aMatrix.n_cols() ;
        aValues.set_size( tN );
 
        for( size_t k=0; k<tN; ++k )
        {
            aValues( k ) = tValues( k );
        }
    }
 
//------------------------------------------------------------------------------
}
 
#endif //BELFEM_FN_AR_EIGEN_HPP