Boost.Multiprecision/Eigen interop challenges

Summary

Currently, Eigen's Eigensolver cannot interoperate with Boost.Multiprecision.

Environment

• Operating System : Mac
• Architecture : Arm64
• Eigen Version : 3.4.0
• Compiler Version : clang
• Compile Flags : -O0
• Vector Extension : NEON

Minimal Example

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <boost/multiprecision/eigen.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>

using boost::multiprecision::cpp_bin_float_100;
int main() {

    Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic> A = Eigen::Matrix<cpp_bin_float_100, Eigen::Dynamic, Eigen::Dynamic>::Identity(3,3);
    Eigen::EigenSolver<decltype(A)> es;
    es.compute(A, /*computeEigenvectors=*/ false);

    auto eigs = es.eigenvalues();
    for (auto eig : eigs) {
      std::cout << eig << "\n";
    }
}

Steps to reproduce

1. Attempt to compile the reproducer

What is the current bug behavior?

The compile fails with:

/opt/homebrew/opt/llvm/bin/clang++ -O0 -g --std=c++17 -Wfatal-errors -fsanitize=address -fsanitize=undefined -I./include -I/opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/ -I/opt/homebrew/opt/llvm/include  -I/opt/homebrew/opt/boost@1.85/include main.cpp
In file included from main.cpp:2:
In file included from /opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/Dense:1:
In file included from /opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/Core:50:
/opt/homebrew/opt/llvm/bin/../include/c++/v1/complex:730:62: fatal error: no matching function for call to '__constexpr_fabs'
  730 |   _Tp __logbw  = std::__constexpr_logb(std::__constexpr_fmax(std::__constexpr_fabs(__c), std::__constexpr_fabs(__d)));
      |                                                              ^~~~~~~~~~~~~~~~~~~~~
/opt/homebrew/Cellar/eigen/3.4.0_1/include/eigen3/Eigen/src/Eigenvalues/EigenSolver.h:542:74: note: in instantiation of function template specialization 'std::operator/<boost::multiprecision::number<boost::multiprecision::backends::cpp_bin_float<100>>>' requested here
  542 |         ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
      |                                                                          ^

What is the expected correct behavior?

The compilation succeeds.

Anything else that might help

We have discussed this issue in boost.math. Basically, std::complex cannot be instantiated with types that are not float, double, or long double, so we need ADL.

marked this issue as related to #1589

The suggestion in the related issue is to extend Eigen::NumTraits<T> to add a Complex type that defaults to std::complex<T>. Then we need to replace all instances of std::complex<Real> with typename Eigen::NumTraits<Real>::Complex.

I'd be willing to look at a merge-request along those lines.

Valid points. However, this compiles and runs (for me) on MSVC. Also on godbolt with the mainstream linux compilers.

The return type is std::complex<cpp_bin_float_100>, which is apparently OK.

OK, replacing instances of typedef std::complex<Scalar> ComplexScalar; with typedef typename NumTraits<Scalar>::Complex ComplexScalar; works, but a little tweaking is necessary for types that are already complex.

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_complex.hpp>
#include <boost/multiprecision/eigen.hpp>

using T = boost::multiprecision::cpp_bin_float_100;
using CT = boost::multiprecision::cpp_complex_100;

template <> struct Eigen::NumTraits<CT> : public Eigen::NumTraitsImp<CT, T> {
  using Complex = CT;
};
using namespace Eigen;
int main() {
  Eigen::Matrix<CT, Eigen::Dynamic, Eigen::Dynamic> A =
      Eigen::Matrix<CT, Eigen::Dynamic, Eigen::Dynamic>::Identity(3, 3);
  Eigen::ComplexEigenSolver<decltype(A)> es;

  es.compute(A, /*computeEigenvectors=*/false);

  auto eigs = es.eigenvalues();
  for (auto eig : eigs) {
    std::cout << eig << "\n";
  }

We could also do something funky like using ComplexScalar = std::conditional_t<NumTraits<Scalar>::IsComplex, Scalar, std::complex<Scalar>> for minimal effort compatibility. That works without having to touch the numtraits of the boost types.

I think that requires 4 line changes in: EigenSolver.h, ComplexEigenSolver.h, RealSchur.h, ComplexSchur.h.

mentioned in issue #2865 (closed)

@cantonios is this fixed with your latest commit?