Custom scalar types and matrix products

Discussion:

Brad Bell

2013-12-24 13:36:20 UTC

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
)

Gael Guennebaud

2013-12-24 16:06:02 UTC

Permalink

Hi,

I don't know why there is an ambiguous call here since I cannot reproduce
with other custom scalar types, but you can workaround with an explicit
evaluation:

A2 * (B2 * A2).eval();

No overhead since matrix products are evaluated into temporaries anyway.

gael

Post by Brad Bell
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.
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
1. svn checkout https://projects.coin-or.org/svn/CppAD/trunkcppad_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
=============
#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
================
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;
^
--------------------------------------------------------------------------
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
--------------------------------------------------------------------------
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
)

Brad Bell

2013-12-25 12:17:29 UTC

Permalink

I have modified the example
https://projects.coin-or.org/CppAD/browser/trunk/bug/eigen_mul.sh
to demonstrate that this does indeed avoid the warning.

In addition, I fixed the fact that echo_eval was missing from the example.

Thanks.

Post by Gael Guennebaud
Hi,
I don't know why there is an ambiguous call here since I cannot
reproduce with other custom scalar types, but you can workaround with
A2 * (B2 * A2).eval();
No overhead since matrix products are evaluated into temporaries anyway.
gael
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.
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
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
=============
#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
================
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;
^
--------------------------------------------------------------------------
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
--------------------------------------------------------------------------
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
)

Brad Bell

2013-12-27 14:56:38 UTC

Permalink

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;
}

Post by Brad Bell
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.
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
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
=============
#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
================
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;
^
--------------------------------------------------------------------------
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
--------------------------------------------------------------------------
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
)

Gael Guennebaud

2013-12-27 16:04:31 UTC

Permalink

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

Post by Brad Bell
I have isolated the problem to a case that only includes Eigen/Core.
==============================
# include <iostream>
# include <Eigen/Core>
// following is 0 or 1. If 0, there is no warning
# define DEFINE_TEMPLATED_CONSTRUCTOR 1
// ------------------------------------------------------------
------------
class myscalar {
// 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;
}

Continue reading on narkive:

Search results for 'Custom scalar types and matrix products' (Questions and Answers)

replies

What's the difference between rules and norms?

started 2008-04-02 21:54:29 UTC

words & wordplay

replies

describ what causes sickle cell anaemia at the DNA and protein level?

started 2006-05-01 12:53:57 UTC

medicine