Re: [eigen] General design issue

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

2009/11/15 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.

Matrix product of two special matrices? I really don't see the use
case for taking the matrix product of e.g. two BandMatrix... do you
have something in mind?

Benoit

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