I am attempting to use Eigen as a "drop in" replacement for an existing matrix library. As such, I have defined EIGEN_MATRIXBASE_PLUGIN and am extending the methods in MatrixBase. In particular, I'd like to add a "dotstar" method. I can do so by adding

Code: Select all

template<typename OtherDerived>
inline CwiseBinaryOp<internal::scalar_product_op<Scalar,Scalar>, const ArrayWrapper<Derived>, const ArrayWrapper<OtherDerived> > dotstar(const MatrixBase<OtherDerived> &other) const {
     return this->array() * other.array();
}

(I know that the "Scalar,Scalar" part should be replaced with calls to some traits thing, but getting that syntax correct was too much trouble, and this is general enough for what I want). While this doesn't allow for "chained" operations (like (m1.dotstar(m2)).dotstar(m3)), that's okay.

However, the returned value does not act like a matrix, which is problematic. So, I tried the following

Code: Select all

template<typename OtherDerived>
inline MatrixWrapper<CwiseBinaryOp<internal::scalar_product_op<Scalar,Scalar>, const ArrayWrapper<Derived>, const ArrayWrapper<OtherDerived> > > dotstar(const MatrixBase<OtherDerived> &other) const {
     return (this->array() * other.array()).matrix();
}

which does not compile because MatrixWrapper is not forward-declared. So, I forward declared MatrixWrapper (like ArrayWrapper is). Then, I get an "invalidd use of undefined type" when the template is instantiated. While I understand a bit of why the compiler (g++ 4.1.2) is complaining, I'm wondering why this wasn't a problem in MatrixBase for the method array() and if there is any way of fixing this. Or maybe there is a generally better way that I'm not seeing?

What I don't want to do is have "dotstar" return a MatrixBase or similar and thereby lose the compile-time optimizations provided by the Eigen framework.

Thanks,

Christian

hi you can directly return a

CwiseBinaryOp<internal::scalar_product_op<Scalar,Scalar>, const Derived, const OtherDerived>

by returning, e.g,

this->cwiseProduct(other);

Have a look at the impl of cwiseProduct in src/plugins/MatrixCwiseBinaryOps.h for a complete example. In fact you simply want to create an alias to this later.

Thanks for the quick reply.

I can, certainly, change the top code to this->cwiseProduct(other). However, if I return a CwiseBinaryOp, the return value is in "array mode" and multiplications and the like will be treated as arrays. I want it to look like a matrix (so that the code (m1.dotstar(m2))*m3 will result in m1 component-wise multiplied by m2 and the result then multiplied by m3 "matrix-wise."

So, I figured if I "typecast" back to a matrix using ".matrix()" that would do the trick. The return type then changes to a MatrixWrapper which I'm having trouble returning as it is not forward declared. However, even if I forward declare it, the compiler (g++) fails at template instantiation time as MatrixWrapper is not "defined:"

Code: Select all

error: invalid use of undefined type ‘struct MatrixWrapper<Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, const Eigen::ArrayWrapper<Eigen::Matrix<double, -0x000000001, -0x000000001, 0, -0x000000001, -0x000000001> >, const Eigen::ArrayWrapper<Eigen::Matrix<double, -0x000000001, -0x000000001, 0, -0x000000001, -0x000000001> > > >’

error: return type ‘struct MatrixWrapper<Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, const Eigen::ArrayWrapper<Eigen::Matrix<double, -0x000000001, -0x000000001, 0, -0x000000001, -0x000000001> >, const Eigen::ArrayWrapper<Eigen::Matrix<double, -0x000000001, -0x000000001, 0, -0x000000001, -0x000000001> > > >’ is incomplete

Is there a way around this?

More broadly, is there a reason the "array-type" and "matrix-type" operations were separated? Other than the confusion over operator* and operator/, it would seem that other operations have the same interpretation in either mode and the conversions back and forth are a bit of a pain (although I see that they do not cause any run-time slowdowns). I feel like there must be another implication to the separation that I am missing.

Thanks.

-Christian

no no, with the solution I posted you still get a matrix, not an array expression.

Linear algebra matrices and simple array of scalars are fundamentally different, and many operators would have completely different meanings (scalar substraction/addition, product, comparisons, exp, log, abs, cos, sin, etc.). That's why they are better separated. Also not that .array() and .matrix() are 100% free, they have zero cost.