I am working with Eigen 3 to do linear algebra for a numerical algebraic geometry tool, and am using my own complex class wrapped around Boost.Multiprecision's mpfr_float class. I can use Eigen to do LU solves on linear algebra problems using my own type, and it passes my tests. What I am stuck on now is computing the norm of a vector.

Here's my code calling norm():

Code: Select all

Eigen::Matrix<bertini::complex, Eigen::Dynamic, Eigen::Dynamic> A(1,3);
A << bertini::complex("1.0","1.0"), bertini::complex("1.0","1.0"), bertini::complex("1.0","1.0");
boost::multiprecision::mpfr_float n = A.norm();

I have defined the NumTraits for mpfr_float, and for my own type, the bertini::complex. I further implemented what I believed were the necessary free functions on my complex type, including abs and abs2. Yet, the compiler error seems to refer to the abs2 function, and complains that it wants a conversion to mpfr_float from complex. This is despite my abs2(bertini::complex) function returning an mpfr_float.

Here is the compiler error I get:

Code: Select all

In file included from /usr/local/include/eigen3/Eigen/Core:259:
/usr/local/include/eigen3/Eigen/src/Core/MathFunctions.h:227:12: error: no viable conversion from 'bertini::complex' to 'RealScalar' (aka 'number<mpfr_float_backend<0> >')
    return x*x;
           ^~~
/usr/local/include/eigen3/Eigen/src/Core/MathFunctions.h:601:45: note: in instantiation of member function 'Eigen::internal::abs2_impl<bertini::complex>::run' requested here
  return EIGEN_MATHFUNC_IMPL(abs2, Scalar)::run(x);
                                            ^
/usr/local/include/eigen3/Eigen/src/Core/Functors.h:313:93: note: in instantiation of function template specialization 'Eigen::numext::abs2<bertini::complex>' requested here
  EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return numext::abs2(a); }
                                                                                            ^
/usr/local/include/eigen3/Eigen/src/Core/CwiseUnaryOp.h:103:14: note: in instantiation of member function 'Eigen::internal::scalar_abs2_op<bertini::complex>::operator()' requested here
      return derived().functor()(derived().nestedExpression().coeff(rowId, colId));
             ^
/usr/local/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:98:24: note: in instantiation of member function 'Eigen::CwiseUnaryOpImpl<Eigen::internal::scalar_abs2_op<bertini::complex>, const
      Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1>, Eigen::Dense>::coeff' requested here
      return derived().coeff(row, col);
                       ^
/usr/local/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:103:14: note: in instantiation of member function
      'Eigen::DenseCoeffsBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<bertini::complex>, const Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >, 0>::coeff' requested here
      return coeff(rowIndexByOuterInner(outer, inner),
             ^
/usr/local/include/eigen3/Eigen/src/Core/Redux.h:177:15: note: (skipping 1 context in backtrace; use -ftemplate-backtrace-limit=0 to see all)
    res = mat.coeffByOuterInner(0, 0);
              ^
/usr/local/include/eigen3/Eigen/src/Core/Redux.h:330:15: note: in instantiation of member function
      'Eigen::internal::redux_impl<Eigen::internal::scalar_sum_op<boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>, 1> >,
      Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<bertini::complex>, const Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >, 0, 0>::run' requested here
            ::run(derived(), func);
              ^
/usr/local/include/eigen3/Eigen/src/Core/Redux.h:363:16: note: in instantiation of function template specialization
      'Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<bertini::complex>, const Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >
      >::redux<Eigen::internal::scalar_sum_op<boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>, 1> > >' requested here
  return this->redux(Eigen::internal::scalar_sum_op<Scalar>());
               ^
/usr/local/include/eigen3/Eigen/src/Core/Dot.h:115:43: note: in instantiation of member function 'Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_abs2_op<bertini::complex>, const
      Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> > >::sum' requested here
  return numext::real((*this).cwiseAbs2().sum());
                                          ^
/usr/local/include/eigen3/Eigen/src/Core/Dot.h:128:15: note: in instantiation of member function 'Eigen::MatrixBase<Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >::squaredNorm'
      requested here
  return sqrt(squaredNorm());
              ^
src/test/classes/complex_test.cpp:583:43: note: in instantiation of member function 'Eigen::MatrixBase<Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >::norm' requested here
                boost::multiprecision::mpfr_float n = A.norm();
                                                        ^
/usr/local/include/boost/multiprecision/number.hpp:44:41: note: candidate constructor not viable: no known conversion from 'bertini::complex' to 'const
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>, 1> &' for 1st argument
   BOOST_MP_FORCEINLINE BOOST_CONSTEXPR number(const number& e) BOOST_NOEXCEPT_IF(noexcept(Backend(static_cast<const Backend&>(std::declval<Backend>())))) : m_backend(e.m_backend){}
                                        ^
/usr/local/include/boost/multiprecision/number.hpp:192:41: note: candidate constructor not viable: no known conversion from 'bertini::complex' to
      'boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>, 1> &&' for 1st argument
   BOOST_MP_FORCEINLINE BOOST_CONSTEXPR number(number&& r)
                                        ^
/usr/local/include/boost/multiprecision/number.hpp:46:25: note: candidate template ignored: substitution failure [with V = bertini::complex]: no type named 'type' in
      'boost::enable_if_c<false, void>'
   BOOST_MP_FORCEINLINE number(const V& v, typename boost::enable_if_c<
                        ^
/usr/local/include/boost/multiprecision/number.hpp:55:41: note: candidate template ignored: substitution failure [with V = bertini::complex]: no type named 'type' in
      'boost::enable_if_c<false, void>'
   BOOST_MP_FORCEINLINE BOOST_CONSTEXPR number(const V& v, typename boost::enable_if_c<
                                        ^
/usr/local/include/boost/multiprecision/number.hpp:91:41: note: candidate template ignored: could not match 'number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      ExpressionTemplates>' against 'bertini::complex'
   BOOST_MP_FORCEINLINE BOOST_CONSTEXPR number(const number<Backend, ET>& val)
                                        ^
/usr/local/include/boost/multiprecision/number.hpp:95:25: note: candidate template ignored: could not match 'number<type-parameter-0-0, ExpressionTemplates>' against 'bertini::complex'
   BOOST_MP_FORCEINLINE number(const number<Other, ET>& val,
                        ^
/usr/local/include/boost/multiprecision/number.hpp:179:4: note: candidate template ignored: could not match 'expression<type-parameter-0-0, type-parameter-0-1, type-parameter-0-2,
      type-parameter-0-3, type-parameter-0-4>' against 'bertini::complex'
   number(const detail::expression<tag, Arg1, Arg2, Arg3, Arg4>& e, typename boost::enable_if_c<is_convertible<typename detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type, s...
   ^
1 error generated.

I feel like I am missing something. In my searching for a solution to this problem, I saw other similar posts (but none exactly dealing with this particular compiler error) where the user was instructed to consult the MathFunctions.h file. In other posts, it seems like users were instructed to provide ei_abs and other related functions.

I add that my abs() and other functions are given in my bertini:: namespace, the same in which the class is defined. I further provide them as both free and member functions for my complex class.

What am I missing? Or, more precisely, what is the set of (names of) functions I need to provide for my number class in order to get 100% compatibility with Eigen? Looking at the online documentation for NumTraits<> and MathFunctions.h has not been enough for me. Do I have to provide ei_*() functions?

For bonus points, help me provide the correct function for using the Eigen::...::Random(bla) function to generate random matrices.

You probably have to add something like:

Code: Select all

namespace Eigen {
namespace numext {
template<>
struct abs2_impl<bertini::complex>
{
  EIGEN_DEVICE_FUNC
  static inline mpfr_float run(const bertini::complex& x)
  {
    return real(x)*real(x) + imag(x)*imag(x);
  }
};
}
}


In the future, we should replace calls like "numext::abs2(x);" to "using numext::abs2; abs2(x);" so that your overload of abs2 can be accessible.

And/or maybe we could specialize for NumTraits<Scalar>::IsComplex instead of std::complex so that our own implementation can work on arbitrary complex types providing real()/imag() accessors.

btw, why cannot you directly use std::complex<mpfr_float > ?

btw, why cannot you directly use std::complex<mpfr_float > ?

because std::complex only supports float, double, and long double. while you can stuff a boost::multiprecision::mpfr_float in there, the internal fields will convert into low-precision numbers (double, etc) when computing things like abs(). not only that, but it will do so silently, with no warning that the resulting computation is essentially garbage. that is, the std::implementation does not use overloading based on templated type in all functions except for with the three standard floating types. i currently view this as the major failure of the standard double type.

at least it is clearly documented, e.g. at http://en.cppreference.com/w/cpp/numeric/complex, where it says

The effect of instantiating the template complex for any other type [than float, double, long double] is unspecified.

(types from previous paragraph added by me)

As far as the actual fix in Eigen you self-suggest,

In the future, we should replace calls like "numext::abs2(x);" to "using numext::abs2; abs2(x);" so that your overload of abs2 can be accessible.

And/or maybe we could specialize for NumTraits<Scalar>::IsComplex instead of std::complex so that our own implementation can work on arbitrary complex types providing real()/imag() accessors.

I had believed (assumed) before reading MathFunctions that this was the case -- it is not documented that implementation is contrary. As a consequence, I have struggled to get all the bits and pieces into place for my own complex number type. Now I am having problems with negation of a vector... But I will ask about it separately if I cannot solve it myself by reading MathFunctions.hpp.

Perhaps the thing I struggle with the most in extending Eigen to work with bertini::complex is that there is no authoritative and comprehensive list of dependencies of Eigen methods on implementations of specializations. That is, I still yearn to be able to look at a single webpage, and read off exactly what I have to add to my code in the eigen:: namespace in order to achieve 100% compatibility with Eigen, including vectorization and other features. Even the NumTraits documentation is pretty incomplete.

In absence of my desired list, I thank you for your quick reply. For now, I was able to get my code to compile by adding the following in the Eigen:: namespace:

Code: Select all

namespace internal {
      template<>
      struct abs2_impl<bertini::complex>
      {
         static inline boost::multiprecision::mpfr_float run(const bertini::complex& x)
         {
            return real(x)*real(x) + imag(x)*imag(x);
         }
      };
   } // re: namespace internal

If the above code is incorrect or violates principles of Eigen, I would appreciate a response so indicating, so that it doesn't make it into production code. I conjecture that I am not done with my class and Eigen yet, as there are other similar structs I have not implemented for my type... Perhaps you could speak to that? Thanks again for your time!

ggael wrote: we could specialize for NumTraits<Scalar>::IsComplex instead of std::complex so that our own implementation can work on arbitrary complex types providing real()/imag() accessors.

for the record, this was done in changeset 1ff681424639