[eigen] Re: Eigen interop with boost::multiprecision broken

[Ian Bell <ian.h.bell@xxxxxxxxx>]

It seems that at version 3.3.3 things were all good, and since then there have a been quite a few regressions for extended precision things. I'm trying to track down a couple of bugs on my side. My hypothesis is that there is somewhere a precision-destroying cast, but I'm not quite sure.

On Thu, Dec 31, 2020 at 3:47 PM Ian Bell <ian.h.bell@xxxxxxxxx> wrote:

I have no idea whether this is a problem with boost or Eigen, but putting boost::multiprecision into Eigen Array/Matrix worked in older versions of Eigen from circa 2017, but no longer does.  For instance this example (https://www.boost.org/doc/libs/1_75_0/libs/multiprecision/doc/html/boost_multiprecision/tut/eigen.html ; copied below) used to work, but now when I compile with master, I get complaints about the NumTraits. I'm at commit fdf2ee62c5174441076fb64c9737d89bbe102759.  I tried to look into this a bit, but template-fu needs work.

On Visual Studio 2019:

Build started...
1>------ Build started: Project: boosteigen, Configuration: Release x64 ------
1>boosteigentest.cpp
1>C:\Users\ihb\PapersInProgress\Working\03a. UlrichCubicSuperAncillary\code\bld\_deps\boost-src\boost/multiprecision/eigen.hpp(147,23): error C2977: 'Eigen::internal::conj_impl': too many template arguments
1>C:\Users\ihb\PapersInProgress\Working\03a. UlrichCubicSuperAncillary\code\bld\_deps\boost-src\boost/multiprecision/eigen.hpp(147): message : see declaration of 'Eigen::internal::conj_impl'
1>C:\Users\ihb\PapersInProgress\Working\03a. UlrichCubicSuperAncillary\code\bld\_deps\boost-src\boost/multiprecision/eigen.hpp(164,8): error C2977: 'Eigen::internal::conj_impl': too many template arguments
1>C:\Users\ihb\PapersInProgress\Working\03a. UlrichCubicSuperAncillary\code\bld\_deps\boost-src\boost/multiprecision/eigen.hpp(157): message : see declaration of 'Eigen::internal::conj_impl'

And on gcc 9.3.0:

[100%] Building CXX object CMakeFiles/boosteigen.dir/boosteigentest.cpp.o
In file included from /input/boosteigentest.cpp:3:
/bld/_deps/boost-src/boost/multiprecision/eigen.hpp:147:14: error: redeclared with 2 template parameters
  147 |       struct conj_impl;
      |              ^~~~~~~~~
In file included from /input/externals/ChebTools/externals/Eigen/Eigen/Core:164,
                 from /bld/_deps/boost-src/boost/multiprecision/eigen.hpp:10,
                 from /input/boosteigentest.cpp:3:
/input/externals/ChebTools/externals/Eigen/Eigen/src/Core/MathFunctions.h:263:34: note: previous declaration 'template<class Scalar> struct Eigen::internal::conj_impl' used 1 template parameter
  263 | template<typename Scalar> struct conj_impl : conj_default_impl<Scalar> {};
      |                                  ^~~~~~~~~
In file included from /input/boosteigentest.cpp:3:
/bld/_deps/boost-src/boost/multiprecision/eigen.hpp:157:100: error: wrong number of template arguments (2, should be 1)
  157 |       struct conj_impl<boost::multiprecision::detail::_expression_<tag, Arg1, Arg2, Arg3, Arg4>, true>
      |                                                                                                    ^
In file included from /input/externals/ChebTools/externals/Eigen/Eigen/Core:164,
                 from /bld/_deps/boost-src/boost/multiprecision/eigen.hpp:10,
                 from /input/boosteigentest.cpp:3:
/input/externals/ChebTools/externals/Eigen/Eigen/src/Core/MathFunctions.h:263:34: note: provided for 'template<class Scalar> struct Eigen::internal::conj_impl'
  263 | template<typename Scalar> struct conj_impl : conj_default_impl<Scalar> {};
      |                                  ^~~~~~~~~

Example:

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

int main()
{
   using namespace Eigen;
   typedef boost::multiprecision::cpp_complex_quad complex_type;
   //
   // We want to solve Ax = b for x,
   // define A and b first:
   //
   Matrix<complex_type, 2, 2> A, b;
   A << complex_type(2, 3), complex_type(-1, -2), complex_type(-1, -4), complex_type(3, 6);
   b << 1, 2, 3, 1;
   std::cout << "Here is the matrix A:\n" << A << std::endl;
   std::cout << "Here is the right hand side b:\n" << b << std::endl;
   //
   // Solve for x:
   //
   Matrix<complex_type, 2, 2> x = A.fullPivHouseholderQr().solve(b);
   std::cout << "The solution is:\n" << x << std::endl;
   //
   // Compute the error in the solution by using the norms of Ax - b and b:
   //
   complex_type::value_type relative_error = (A*x - b).norm() / b.norm();
   std::cout << "The relative error is: " << relative_error << std::endl;
   return 0;
}