Re: [eigen] Support for true Array

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2009/11/6 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
> Hi,
>
> I now we already have way too much modules/features which have to be
> polished before thinking about new features, but anyway, let me do this
> proposal.
>
> Eigen has a very good support for mathematical matrices. Even though many
> kind of array-like operations can be performed on our matrices via the
> cwise() prefix, we are still lacking a true support for array of scalars. In
> particular, it would be very convenient to have an Array class that you
> could use just like a scalar value, all operations being berformed
> coefficient wise. I often thought about some options to Matrix/MatrixBase
> for that, but this is rather ugly and prone to many  mistakes.

OK, assuming there is a reasonable implementation for that (and there
is, i agree with your proposal) then this is useful to have in Eigen.
Even if we wanted to claim that Eigen is 100% focused on matrices,
this would still be useful for the following reason. Eigen currently
only does "vertical" i.e. "per-object" vectorization. Which is perfect
for large objects, but not perfect for small objects which may not be
efficiently divided into packets. With your proposal, if a have a
million of Vector3f's, i could instead do a Vector<array,3,1> and all
operations on it would get automatically vectorized, which is
supercool as per-object vectorization is hard/impossible in this case.
This advantage of horizontal vectorization is going to be even more
important with future SIMD such as Larrabee as the packet size gets
bigger.

>
> So here the idea is to mimic what is done for AutoDiffScalar,
> DiagonalMatrix, Quaternion, etc.  Basically, the principle is wrap a
> MatrixBase object (a matrix or an expression) into a wrapper class providing
> a different, specialized API. In practice this is implemented by defining an
> ArrayBase class supporting basics operators and the math functions of the
> standard library (no a.sin(), but std::sin(a)). Then we define two derived
> classes providing the actual coefficients:
> - A class Array<Scalar,Rows,Cols> providing storage and convenient ctors. It
> would basically store a Matrix<Scalar,RowsCols>.
> - A class ArrayWrapper<Expression> wrapping any MatrixBase expression as an
> Array (mostly for internal use to enable expression templates and, e.g., to
> view an existing Matrix as an Array, etc.).
>
> To make thing even more clear, here is the implementation of operator* and
> std::sin:
>
> template<typename OtherDerived>
> ArrayWrapper<CwiseBinaryOp<ei_scalar_multiple_op<Scalar>, typename
> Derived::Coefficients, typename OtherDerived::Coefficients>
> ArrayBase<Derived>::operator*(const ArrayBase<OtherDerived>& other) const
> {
>   return matrix() * other.matrix();
> }

OK, I understand. The cool thing in your idea is that we don't have to
code again all the expression templates, instead we use existing
expression classes wrapped in ArrayWrapper.

That makes me wonder, perhaps the exact same idea should be applied
also to Quaternion where a recent thread on this list suggested that
expression templates are useful too (as non-gcc compilers are not able
to understand that everything fits in registers).

>
> namespace std {
>
> template<typename Derived>
> ArrayWrapper<CwiseUnaryOp<ei_scalar_sin_op<Scalar>, typename
> Derived::Coefficients>
> sin(const ArrayBase<Derived>& a)
> {
>   return a.matrix().cwis().sin();
> }
>
> }

I see. I understand how it makes sense to define std::sin here but i
think we should also define ei_sin so that array type can be
seamlessly used as the scalar type with the rest of Eigen.

> Following the AutoDiffScalar implementation, all of this can be done using a
> couple of macros to reduce coding and maintainability efforts.
>
> A good thing is that can be done without any change in the core of Eigen,
> and so one can do it totally outside Eigen.

Yes, this is an extra bonus.
However, since this doesn't seem like too much extra code, perhaps
this should go into the existing Array module. That would also solve a
problem with naming!

>
> What do you think ?

*kisses*

Benoit

>
> gael.
>
>