CG solver doesn't compile for __float128 element type because of ambiguous overload for sqrt(RealScalar)

The following program fails to compile with the current git

#include <Eigen/Sparse>

int main()
{
  using DType = __float128;
  long L = 1000;
  using SparseMatrixType = Eigen::SparseMatrix<DType, Eigen::RowMajor, int64_t>;
  SparseMatrixType I_Q(L, L);
  I_Q.reserve(Eigen::Matrix<long, Eigen::Dynamic, 1>::Constant(L, 386));
  Eigen::ConjugateGradient<SparseMatrixType, Eigen::Lower|Eigen::Upper> cg;
  cg.compute(I_Q);
  Eigen::Matrix<DType, Eigen::Dynamic, 1> c(L);
  Eigen::Matrix<DType, Eigen::Dynamic, 1> x(L);
  std::fill(c.begin(), c.end(), 1);
  x = cg.solve(c);
}

using g++ -std=gnu++17 -Ieigen -c test2l.cc with g++ 8.5.0

The error is

(...)
eigen2/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h:91:19: error: call of overloaded ‘sqrt(RealScalar)’ is ambiguous
   tol_error = sqrt(residualNorm2 / rhsNorm2);
               ~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /usr/include/features.h:438,
                 from /usr/include/c++/8/x86_64-redhat-linux/bits/os_defines.h:39,
                 from /usr/include/c++/8/x86_64-redhat-linux/bits/c++config.h:2470,
                 from /usr/include/c++/8/cmath:41,
                 from eigen2/Eigen/src/Core/util/Macros.h:679,
                 from eigen2/Eigen/Core:19,
                 from eigen2/Eigen/SparseCore:11,
                 from eigen2/Eigen/Sparse:26,
                 from test2l.cc:1:
/usr/include/bits/mathcalls.h:143:1: note: candidate: ‘double sqrt(double)’
 __MATHCALL (sqrt,, (_Mdouble_ __x));
 ^~~~~~~~~~
In file included from eigen2/Eigen/src/Core/util/Macros.h:679,
                 from eigen2/Eigen/Core:19,
                 from eigen2/Eigen/SparseCore:11,
                 from eigen2/Eigen/Sparse:26,
                 from test2l.cc:1:
/usr/include/c++/8/cmath:467:3: note: candidate: ‘constexpr long double std::sqrt(long double)’
   sqrt(long double __x)
   ^~~~
/usr/include/c++/8/cmath:463:3: note: candidate: ‘constexpr float std::sqrt(float)’
   sqrt(float __x)
   ^~~~

If I replace the call to sqrt(x) in the CG with the C version (::sqrt) it works. According to Antonio, it should be using Eigen::numext::sqrt, but replacing using std::sqrt with using Eigen::numext::sqrt doesn't compile, either.

Possibly relevant here is Boost's float128 type. I note that the following program, which uses Boost's float128 instead of GNUs __float128, compiles:

#include <boost/cstdfloat.hpp> // For float_64_t, float128_t. Must be first include!
#include <boost/multiprecision/float128.hpp>

#include <Eigen/Sparse>
 
int main()
{
  using DType = boost::multiprecision::float128;
  long L = 1000;
  using SparseMatrixType = Eigen::SparseMatrix<DType, Eigen::RowMajor, int64_t>;
  SparseMatrixType I_Q(L, L);
  I_Q.reserve(Eigen::Matrix<long, Eigen::Dynamic, 1>::Constant(L, 386));
  Eigen::ConjugateGradient<SparseMatrixType, Eigen::Lower|Eigen::Upper> cg;
  cg.compute(I_Q);
  Eigen::Matrix<DType, Eigen::Dynamic, 1> c(L);
  Eigen::Matrix<DType, Eigen::Dynamic, 1> x(L);
  std::fill(c.begin(), c.end(), 1);
  x = cg.solve(c);
}

Also, the CG code already uses the using std:sqrt pattern to try to use ADL. However, __float128 is of course a GNU extension that apparently doesn't have an overload in the std namespace (not sure why?). I note that Boost's documentation explicitly mentions compatibility with std::sqrt as one of the selling points of their float128 type.

I suppose Eigen could either support __float128 or rely on users to provide a suitable wrapper type, such as Boost's, but it would be kind of nice if this worked out of the box without requiring another library in Eigen.

Edited by Godmar Back