Search code examples
c++linear-algebralapackarmadilloblas

Can Armadillo do eig_gen on Mat<float>?


I am using Armadillo. I have these variables:

arma::Mat<float> m_matrix;
arma::cx_vec     m_eigenvalues;
arma::cx_mat     m_eigenvectors;

I want to do this:

void calculate_eigens ()
{
    arma :: eig_gen (m_eigenvalues, m_eigenvectors, m_matrix);
}

but the function prototype for eig_gen does not allow these parameter types, it wants the third argument to be a double-precision matrix.

Here are the gcc errors:

error: no matching function for call to ‘eig_gen(arma::cx_vec&, arma::cx_mat&, arma::Mat&)’

note: candidate: template typename arma::enable_if2<arma::is_supported_blas_type::value, arma::Col<std::complex > >::result arma::eig_gen(const arma::Base<typename T1::elem_type, T1>&)

armadillo_bits/fn_eig_gen.hpp:25:1: note: template argument deduction/substitution failed:

note: candidate expects 1 argument, 3 provided

This fixes it.

void calculate_eigens ()
{
    auto tmp = arma::conv_to<arma::Mat<double>>::from(m_matrix);

    arma :: eig_gen (m_eigenvalues, m_eigenvectors, tmp);
}

I don't want to do this conversion. I looked into Armadillo's source and it seems that eig_gen defers its main work to this function:

  template<typename eT>
  inline
  void
  geev(char* jobvl, char* jobvr, blas_int* N, eT* a, blas_int* lda, eT* wr, eT* wi, eT* vl, blas_int* ldvl, eT* vr, blas_int* ldvr, eT* work, blas_int* lwork, blas_int* info)
    {
    arma_type_check(( is_supported_blas_type<eT>::value == false ));

    if(is_float<eT>::value)
      {
      typedef float T;
      arma_fortran(arma_sgeev)(jobvl, jobvr, N, (T*)a, lda, (T*)wr, (T*)wi, (T*)vl, ldvl, (T*)vr, ldvr, (T*)work, lwork, info);
      }
    else
    if(is_double<eT>::value)
      {
      typedef double T;
      arma_fortran(arma_dgeev)(jobvl, jobvr, N, (T*)a, lda, (T*)wr, (T*)wi, (T*)vl, ldvl, (T*)vr, ldvr, (T*)work, lwork, info);
      }
    }

The eT*a parameter is passed in from Mat<T>::get_ref on m_matrix, so it seems to me like the library author intended for eig_gen to work happily on Mat<float> as well as Mat<double> but it doesn't work in practice.

How do I calculate the eigenvalues/eigenvectors on a Mat<float> without first converting to Mat<double>?


Solution

  • You need to ensure that all matrices and vectors have the same precision, not just the input matrix. The correct usage for single precision would be:

    arma::fmat     X;
    arma::cx_fvec  eigvals;
    arma::cx_fmat  eigvecs;
    
    // ... fill X with values ...
    
    arma::eig_gen(eigvals, eigvecs, X);