CODEWITHSUNDEEP
c++boosteigen
carsten
carsten

Reputation: 1365

template resolution fails with boost multiprecision + Eigen

I'm trying to use Eigen to do an eigenvector decomposition on mutliprecision floating types from boost::multiprecision. I'm starting with a very simple example that I've copied together from different sources. Here's the code:

#include <boost/multiprecision/cpp_dec_float.hpp>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/LU>
#include <eigen3/Eigen/Eigenvalues>
#include <iostream>

typedef boost::multiprecision::number<boost::multiprecision::cpp_dec_float<100> > SuperFloat;
typedef std::complex<SuperFloat> SuperComplex;

// this is the first fix I came up with to overcome the problem
// that multiprecision doesn't come with an int() operator
namespace Eigen {
  namespace internal {
    template<typename NewType> 
    struct cast_impl<SuperFloat,NewType> {
      static inline NewType run(const SuperFloat& x)  {
        return x.convert_to<NewType>();
      }
    };
  }
}

typedef Eigen::Matrix<SuperFloat, Eigen::Dynamic, Eigen::Dynamic> EigenMatrixR;
typedef Eigen::Matrix<SuperFloat, Eigen::Dynamic, 1 > EigenVectorR;
typedef Eigen::Matrix<SuperComplex, Eigen::Dynamic, Eigen::Dynamic> EigenMatrixC;
typedef Eigen::Matrix<SuperComplex, Eigen::Dynamic, 1 > EigenVectorC;

int main(){
    int size = 10;   

    EigenMatrixR A = EigenMatrixR::Identity(size, size);

    Eigen::EigenSolver<EigenMatrixR> es(A);
    std::cout << "The eigenvalues of A are:" << std::endl << es.eigenvalues() << std::endl;
    std::cout << "The matrix of eigenvectors, V, is:" << std::endl << es.eigenvectors() << std::endl << std::endl;
    SuperComplex lambda = es.eigenvalues()[0];
    std::cout << "Consider the first eigenvalue, lambda = " << lambda << std::endl;
    EigenVectorC v = es.eigenvectors().col(0);
    std::cout << "If v is the corresponding eigenvector, then lambda * v = " << std::endl << lambda * v << std::endl;
    std::cout << "... and A * v = " << std::endl << A.cast<SuperComplex>() * v << std::endl << std::endl;
    EigenMatrixC D = es.eigenvalues().asDiagonal();
    EigenMatrixC V = es.eigenvectors();
    std::cout << "Finally, V * D * V^(-1) = " << std::endl << V * D * V.inverse() << std::endl;

    return 0;
}

I've overcome the first few pitfalls (like the missing int() operator for the boost::multiprecision types, which use a convert_to method instead), but now I've reached the point where the compiler just spits out error messages about failed template resolutions.

The complete error log is quite long (I put it on pastebin: http://pastebin.com/a2R0NDSA), but the first error is this:

/usr/include/eigen3/Eigen/src/Eigenvalues/EigenSolver.h:549:43: error: no matching function for call to ‘cdiv(boost::multiprecision::detail::expression<boost::multiprecision::detail::negate, boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<100u> >, void, void, void>, boost::multiprecision::detail::expression<boost::multiprecision::detail::negate, boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<100u> >, void, void, void>, Eigen::EigenSolver<Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<100u> >, -1, -1> >::Scalar&, Eigen::EigenSolver<Eigen::Matrix<boost::multiprecision::number<boost::multiprecision::backends::cpp_dec_float<100u> >, -1, -1> >::Scalar&)’
             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
                                           ^
/usr/include/eigen3/Eigen/src/Eigenvalues/EigenSolver.h:422:22: note: candidate: template<class Scalar> std::complex<_Tp> Eigen::cdiv(const Scalar&, const Scalar&, const Scalar&, const Scalar&)
 std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi)

In other words, Eigen tries to use a function taking four scalars, but boost provides expression templates for the first two (and the compiler refuses to implicitly convert those to scalars).

Am I on the right path, or is this endeavour futile? Any suggestions on how to proceed to teach Eigen how to utilize the boost::multiprecision types?

Update

Thanks to the helpful comments below this question, I've been able to resolve the issue by turning off expression templates.

typedef boost::multiprecision::cpp_dec_float<50> mp_backend;
typedef boost::multiprecision::number<mp_backend, boost::multiprecision::et_off> SuperFloat;

The remaining error message about a failing template resolution for check_in_range can be fixed like this:

namespace boost{
  namespace multiprecision {
    namespace default_ops{                                                                          
      template <> inline bool check_in_range<SuperComplex,long double>(const long double& t){
        return false;
      }
    }
  }
}

Upvotes: 3

Views: 614

Answers (1)

Walter
Walter

Reputation: 45424

Your error was caused by boost::multiprecision to return a expression template object from the unary minus operator, rather than another scalar (i.e. the same boost::multiprecision number type).

The obvious solution is to use a multi-precision type that avoids expression templates, either by going for another type altogether (which may be faster anyway) or by switching off the expression templates for boost::multiprecision, see here.

Upvotes: 5

Related Questions