Re: [eigen] Custom scalar types and matrix products

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

ok, thanks for the finding. In my opinion such a constructor should be declared explicit because otherwise it is quite obvious such a constructor is a source for troubles. I don't have any workaround for the 3.2 branch (except calling eval explicitly) but this precise issue is already fixed by the evaluators that should make it for the 3.3 release (currently developed in the eigen-evaluators fork).

gael

On Fri, Dec 27, 2013 at 3:56 PM, Brad Bell <bradbell@xxxxxxxxxx> wrote:

I have  isolated the problem to a case that only includes Eigen/Core.
The source code is included below:
==============================

# include <iostream>
# include <Eigen/Core>

// following is 0 or 1. If 0, there is no warning
# define DEFINE_TEMPLATED_CONSTRUCTOR 1

// ------------------------------------------------------------------------
class myscalar {
public:
     // data
     double value_;
     // constructors
     myscalar(void)
     {    value_ = 0.0; }
     myscalar(const double& value)
     {    value_ = value; }

# if DEFINE_TEMPLATED_CONSTRUCTOR
     template <class T> myscalar(const T& value)
     {    value_ = double(value); }
# endif

     // binary operators
     myscalar operator+(const myscalar& right) const
     {    myscalar result = value_ + right.value_;
          return result;
     }
     myscalar operator*(const myscalar& right) const
     {    myscalar result = value_ * right.value_;
          return result;
     }
     bool operator!=(const myscalar& right) const
     {    bool result = value_ != right.value_;
          return result;
     }

     // assignments
     myscalar& operator=(const double& value)
     {    value_ = value;
          return *this;
     }

     // computed assignment operators
     myscalar& operator +=(const myscalar& right)
     {    value_ += right.value_;
          return *this;
     }

};
// ------------------------------------------------------------------------

int main() {
     using Eigen::Matrix;
     using Eigen::Dynamic;

     Matrix<myscalar, Dynamic, Dynamic> A(1,1);
     Matrix<myscalar, Dynamic, Dynamic> B(1,1);
     Matrix<myscalar, Dynamic, Dynamic> C(1,1);
     Matrix<myscalar, Dynamic, Dynamic> D(1,1);

     A(0,0) = 1.0;
     B(0,0) = 2.0;
     C(0,0) = 3.0;
     D      = A * (B * C);

     if( D(0,0) != myscalar(6.0) )
     {    std::cout << "eigen_mul: Error" << std::endl;
          return 1;
     }
     std::cout << "eigen_mul: OK" << std::endl;
     return 0;

}

On 12/24/2013 5:36 AM, Brad Bell wrote:

I am having trouble using the templated scalar type CppAD::AD<Base> with eigen's matrix multiply.

To be specific, the code eigen_mul.cpp (included below)
is generating the error message eigen_mul.stderr (included below).
I have reformatted the error message output so that it is easier to read.

A script that:
generates the file build/eigen_mul.cpp,
runs the g++ compiler,
and generates the error message,
can be found at
https://projects.coin-or.org/CppAD/browser/trunk/bug/eigen_mul.sh

You can run this script using the following steps:
    1. svn checkout https://projects.coin-or.org/svn/CppAD/trunk cppad_trunk
    2. cd cppad_trunk
    3. mkdir build
    4. cd build
    5. cmake ..
    6. cd ../bug
Then modify the script eigen_mul.sh line
     -I$HOME/prefix/eigen/include
to be the location of the eigen include directory on your machine.
Then execute the command
    ./eigen_mul.sh


eigen_mul.cpp:
=============
#include <Eigen/Core>
#include <cppad/cppad.hpp>

int main() {
     using Eigen::Matrix;
     using Eigen::Dynamic;
     typedef CppAD::AD<double> AScalar;

     Matrix<AScalar, Dynamic, Dynamic> A2(1,1);
     Matrix<AScalar, Dynamic, Dynamic> B2(1,1);

     A2(0,0) = 1.0;
     B2(0,0) = 2.0;

     A2 * B2 * A2;

     return 0;
}

eigen_mul.stderr
================

eigen_mul.cpp:15:12:
warning: ISO C++ says that these are ambiguous,
even though the worst conversion for the first is better than the
worst conversion for the second: [enabled by default]
  A2 * B2 * A2;
            ^
--------------------------------------------------------------------------
/home/bradbell/prefix/eigen/include/Eigen/src/Core/GeneralProduct.h:571:1:
note: candidate 1: const typename

Eigen::ProductReturnType<Derived, OtherDerived>::Type
Eigen::MatrixBase<Derived>::operator*(
     const Eigen::MatrixBase<OtherDerived>&) const

OtherDerived = Eigen::Matrix<CppAD::AD<double>, -1, -1>;
Derived      = Eigen::GeneralProduct<
     Eigen::Matrix<CppAD::AD<double>, -1, -1>,
     Eigen::Matrix<CppAD::AD<double>, -1, -1>, 5>;

typename Eigen::ProductReturnType<Derived, OtherDerived>::Type =
Eigen::GeneralProduct<
     Eigen::GeneralProduct<
          Eigen::Matrix<CppAD::AD<double>, -1, -1>,
          Eigen::Matrix<CppAD::AD<double>, -1, -1>, 5>,
     Eigen::Matrix<CppAD::AD<double>, -1, -1>, 5>

MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const


--------------------------------------------------------------------------
/home/bradbell/prefix/eigen/include/Eigen/src/Core/ProductBase.h:198:1:
note: candidate 2: const

Eigen::ScaledProduct<Derived>
Eigen::operator*(
     const Eigen::ProductBase<Derived, _Lhs, _Rhs>&,
     const typename Derived::Scalar&)

Derived = Eigen::GeneralProduct<
     Eigen::Matrix<CppAD::AD<double>, -1, -1>,
     Eigen::Matrix<CppAD::AD<double>, -1, -1>, 5>;
Lhs = Eigen::Matrix<CppAD::AD<double>, -1, -1>;
Rhs = Eigen::Matrix<CppAD::AD<double>, -1, -1>;

typename Derived::Scalar = CppAD::AD<double>

operator*(
     const ProductBase<Derived,Lhs,Rhs>& prod,
     const typename Derived::Scalar& x
)