This forum has been archived. All content is frozen. Please use KDE Discuss instead.

right scalar multiplication by real or integral types

Tags: None
(comma "," separated)
danielbrake
Registered Member
Posts
15
Karma
0
Hi Eigen,

I use Eigen for linear algebra in Bertini, a numerical multivariate polynomial system solver. I have enjoyed great success so far with Eigen 3.2.x, with my custom complex type using Boost.Multiprecision's mpfr_float type.

I have two issues. One is scalar multiplication of matrices of my complex type against integers. The other is a compiler error from Eigen 3.3-beta1.

The first issue: creation of temporaries when converting integral types to my complex type are expensive, and this must be avoided. I am concerned about scalars magically changing type from integers to floating points, incurring the cost of construction of a multiple precision complex number.

Second, when I compile my library against 3.3-beta1, I am getting ambiguous operator * for right scalar multiplication. Not left, just right.

I have not put together a minimal example, as that would take several files. See github.com/bertiniteam/b2 if you want full code. For now, believe that I have a typical complex class defined as one would. I then create the NumTraits, etc, for Eigen. Below is some code that demonstrates my errors:

--

The code enabling scalar multiplication by other types than the main data type, or its Real or Nested type. see viewtopic.php?f=74&t=111176

Code: Select all
      //int
      template<>  // these two (for int) apparently aren't needed, because the lower pair enable int --> mpfr_float conversion
      struct scalar_product_traits<int,bertini::complex>
      {
          enum { Defined = 1 };
          typedef bertini::complex ReturnType;
      };

      template<>  // these two (for int) apparently aren't needed, because the lower pair enable int --> mpfr_float conversion
      struct scalar_product_traits<bertini::complex, int>
      {
          enum { Defined = 1 };
          typedef bertini::complex ReturnType;
      };

      template<> // this pair definitely needed regardless of version
      struct scalar_product_traits<bertini::mpfr_float,bertini::complex>
      {
          enum { Defined = 1 };
          typedef bertini::complex ReturnType;
      };

      template<>  // this pair definitely needed regardless of version
      struct scalar_product_traits<bertini::complex, bertini::mpfr_float>
      {
          enum { Defined = 1 };
          typedef bertini::complex ReturnType;
      };


--

The test code which fails:

Code: Select all
BOOST_AUTO_TEST_CASE(scalar_multiplication_complex_int)
   {
      using data_type = bertini::complex;
      
      data_type q(1);
      int a(1);

      Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> A(2,2);
      A << data_type(2), data_type(1), data_type(1), data_type(2);
      
      Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> B = a*A; // needs scalar_product_traits in 3.2/3.3, but apparently with a conversion to mpfr_float, no!!!!!
      Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> C = A*a; // resolves  correctly in 3.2, ambiguous in 3.3
   }


For the errors, it's the scalar multiplication on the right that I am working on. The error exists independently of the scalar_product_traits templates.

--

The errors I am trying to solve, again, are ambiguous operator errors from scalar multiplication on the right:

Code: Select all
test/classes/eigen_test.cpp:378:65: error: use of overloaded operator '*' is ambiguous (with operand types 'Eigen::Matrix<data_type, Eigen::Dynamic,
      Eigen::Dynamic>' and 'int')
                Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> C = A*a;
                                                                             ~^~
/Users/ofloveandhate/software/eigen-3.3beta1/Eigen/src/Core/util/XprHelper.h:470:3: note: candidate function
  operator*(const OtherScalar& scalar) const
  ^
/Users/ofloveandhate/software/eigen-3.3beta1/Eigen/src/Core/../plugins/CommonCwiseUnaryOps.h:56:1: note: candidate function
operator*(const Scalar& scalar) const


The deeper issue I see here is that I would really like to be able to use scalar multiplication by integers of various types, as construction of the multiprecision complex type from integers is expensive, including a heap allocation from mpfr. So resolving the ambiguous operator is the biggest hurdle I am having right now for users wanting Eigen 3.3, but getting integer-multipleprecision arithmetic enabled without temporaries holding those integers is a high priority too.
--

please note (in an effort to be complete) that the scalar_product_traits are needed for multiplication on the left, as without, i get error (eigen 3.2.7):

Code: Select all
/usr/local/include/eigen3/Eigen/src/Core/Functors.h:520:60: error: no type named 'ReturnType' in 'Eigen::internal::scalar_product_traits<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >'
  typedef typename scalar_product_traits<Scalar1,Scalar2>::ReturnType result_type;
          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~
/usr/local/include/eigen3/Eigen/src/Core/util/Meta.h:133:64: note: in instantiation of template class
      'Eigen::internal::scalar_multiple2_op<bertini::complex, boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0,
      allocate_dynamic>, boost::multiprecision::expression_template_option::et_off> >' requested here
    static has_tr1_result      testFunctor(T const *, typename T::template result<T(ArgType)>::type const * = 0);
                                                               ^
/usr/local/include/eigen3/Eigen/src/Core/util/Meta.h:137:32: note: while substituting deduced template arguments into function template 'testFunctor'
      [with T = Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >]
    enum {FunctorType = sizeof(testFunctor(static_cast<Func*>(0)))};
                               ^
/usr/local/include/eigen3/Eigen/src/Core/CwiseUnaryOp.h:41:20: note: in instantiation of template class
      'Eigen::internal::result_of<Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> > (bertini::complex)>' requested here
  typedef typename result_of<
                   ^
/usr/local/include/eigen3/Eigen/src/Core/util/XprHelper.h:349:56: note: in instantiation of template class
      'Eigen::internal::traits<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >, Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> > >' requested here
template<typename Derived, typename XprKind = typename traits<Derived>::XprKind>
                                                       ^
/usr/local/include/eigen3/Eigen/src/Core/CwiseUnaryOp.h:93:22: note: in instantiation of default argument for
      'dense_xpr_base<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >, Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> > >' required here
  : public internal::dense_xpr_base<CwiseUnaryOp<UnaryOp, XprType> >::type
                     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/usr/local/include/eigen3/Eigen/src/Core/CwiseUnaryOp.h:60:10: note: in instantiation of template class
      'Eigen::CwiseUnaryOpImpl<Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >, Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1>, Eigen::Dense>' requested here
  public CwiseUnaryOpImpl<UnaryOp, XprType, typename internal::traits<XprType>::StorageKind>
         ^
test/classes/eigen_test.cpp:377:65: note: in instantiation of template class 'Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<bertini::complex,
      boost::multiprecision::number<boost::multiprecision::backends::mpfr_float_backend<0, allocate_dynamic>,
      boost::multiprecision::expression_template_option::et_off> >, Eigen::Matrix<bertini::complex, -1, -1, 0, -1, -1> >' requested here
                Eigen::Matrix<data_type, Eigen::Dynamic, Eigen::Dynamic> B = a*A;


this is true when using the current dev 3.3-beta1, the unstable dev, and 3.2. Hence, the scalar_product_traits are necessary.
--

Thanks a ton for reading, and helping!
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
Indeed, so far Eigen provides scalar operator* only with Scalar and Real types, so in your case "complex<mp>" and "mp". I guess that both your complex class and "mp" have a ctor taking an integer? Then I do understand the ambiguity, but I don't understand why it does not show up with 3.2! Prototypes are the same...

Actually, fixing the second issue (conversion) would also fix the former, so maybe we could look at it first. One possibility would be to declare a generic operator* with SFINAE to enable it based on scalar_product_traits, for instance:

Code: Select all
template <typename T>
EIGEN_DEVICE_FUNC inline friend
typename internal::enable_if<internal::scalar_product_traits<T,Scalar>::Defined,
                             const CwiseUnaryOp<internal::scalar_multiple2_op<Scalar,T>, const Derived> >::type
operator*(const T& scalar, const StorageBaseType& matrix)
{
  return CwiseUnaryOp<internal::scalar_multiple2_op<Scalar,T>, const Derived>(
            matrix.derived(), internal::scalar_multiple2_op<Scalar,T>(scalar) );
}


You can try it by adding it in plugins/CommonCwiseUnaryOps.h (and doing the same for the non-friend version, that is left multiply).

Let me known how this is working.
User avatar
ggael
Moderator
Posts
3447
Karma
19
OS
I've pushed such changes: https://bitbucket.org/eigen/eigen/commits/d896d06d500d, so you can simply "hg pull -u", and everything should be fine!
danielbrake
Registered Member
Posts
15
Karma
0
the mentioned commit https://bitbucket.org/eigen/eigen/commits/d896d06d500d does indeed solve the ambiguous operator problem. Thanks!
branwik
Registered Member
Posts
1
Karma
0
Rules of multiplication can be found on Aztekium site - http://www.aztekium.pl/master.
Very helpful, really :)


Bookmarks



Who is online

Registered users: Bing [Bot], blue_bullet, Google [Bot], rockscient, Yahoo [Bot]