Re: [eigen] General design issue

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Sun, Nov 15, 2009 at 5:00 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

Just one thing though: I am still not convinced that we want
arithmetic operators (let alone _expression_ templates) for all of the
special matrix classes. It sure would be nice, but I am not convinced
that it's worth the extra work and complexity.

Well, yes and no.

Let's assume we agree that supporting matrix products is important. The problem is that in most cases you don't only want simple product, but you also want to be able to take the adjoint of the factors, do some scaling and maybe accumulate the result at the same time. This is not an hazard if the BLAS GEMM routine does all that in a single call. To make thing clear let me recall that GEMM allow to do:

C = C + alpha * op(A) * op(B)

However to support this kind of high level products we need to be able to build an _expression_ tree containing Transpose, CwiseUnaryOp, and Product. Whence the "no".

Now let's consider an arbitrary fancy storage Foo. A funny thing with my proposal is that if you implement the product for the Foo storage as we did for the dense path (i.e., via a kind of GEMM routine), then you even don't need to implement the specialization of {Transpose|CwiseUnaryOp}Impl for the storage Foo ! Whence the "yes".

Nevertheless, both Transpose and CwiseUnaryOp are usually trivial to implement. The most difficult ones are Block,  CwiseBinary, and partial reductions which are less useful.

Gael.
 

I agree about all the rest. And it's true that your changes, which are
useful anyway, would also be the first thing to do if one wanted
arithmetic operators for special matrices. Finally, such operators and
_expression_ templates can be added at a later date.

Benoit

2009/11/13 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>
>
> On Fri, Nov 13, 2009 at 9:46 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> wrote:
>>
>> 2009/11/13 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>> >
>> >
>> > On Fri, Nov 13, 2009 at 8:58 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
>> > wrote:
>> >>
>> >> Hi,
>> >>
>> >> Wow, your mail took me a lot of effort to understand, I think it's
>> >> mostly ok now...
>> >>
>> >
>> > sorry for that, I should have been more explicit.
>> >
>> >>
>> >> First of all, I agree with your summary of current limitations. Yes, I
>> >> agree we must take care of that before beta1.
>> >>
>> >> Let me rephrase some stuff in my own baby words to see if i got it
>> >> right. The 2 basic problems you want to address are:
>> >> * code sharing between dense and sparse xpr
>> >> * in xpr trees, only the leaves should be affected by what the leaves
>> >> actually are.
>> >>
>> >> The "naive" approach to that would probably be to let e.g. Transpose
>> >> and SparseTranspose inherit a common base AnyTranspose, making it a
>> >> curious base so AnyTranspose<Derived>... then .transpose() methods
>> >> should return objects of AnyTranspose<Derived> which is ugly and may
>> >> (out of empirical experience) give compiler warnings about strict
>> >> aliasing when casting the returned object to the derived type.
>> >>
>> >> So your solution, if I understand well, consists instead in letting
>> >> Transpose itself be a template in one more integer parameter, which
>> >> you call StorageType, which could have special values Dense, Sparse,
>> >> Band...
>> >>
>> >> At this point, I think I have understood how you solve the "the rest
>> >> of the tree doesn't completely change if the leaves are different"
>> >> problem.
>> >>
>> >> But can you enlighten me as to how you solve the "code sharing"
>> >> problem? If Transpose<Dense> and Transpose<Sparse> can share code,
>> >> where do you put that? In a common base class TransposeBase? So that
>> >> problem and solution is actually completely orthogonal to the above?
>> >
>> > What I suggest is a bit different. Let's pick an example:
>> >
>> > template<MatrixType> Transpose : public TransposeImpl<MatrixType,
>> > ei_traits<MatrixType>::StorageType>
>> > {
>> >  int rows();
>> >  int cols();
>> >  MatrixType::Nested m_matrix;
>> >  // etc.
>> > }
>> >
>> > then TranposeImpl<MatrixType,Dense> would inherit MatrixBase<Transpose<>
>> > >
>> > and implements the coeff() methods while
>> > TransposeImpl<MatrixType,Sparse>
>> > would inherit SparseMatrixBase<> and implement the InnerIterators
>> > mechanism.
>>
>> ok so actually in your scheme, the code sharing happens in the derived
>> class Transpose, while the non-shared implementation happens in the
>> base class TransposeImpl.... clever!
>>
>> What I like the most here is that the class that is used as a node in
>> the xpr tree, Transpose, doesn't get any additional template
>> parameter: our xpr trees were already complex enough like that...
>>
>> Are you going to set up a new fork for that? And add it to the
>> todo-3.0 page? It seems quite a long project to complete.
>
> yep, this is definitely not a trivial change, so of course I'll start a fork
> for that and we'll see!
>
> gael
>
>>
>> Benoit
>>
>>
>
>
>
>
>