Re: [eigen] General design issue

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

one more remark:

perhaps the most urgent / worthwile effort in this direction would be
to ensure that our architecture / API allows doing that at a later
date. Indeed, i think that we agree that these features are not by
themselves 3.0 blockers, but it would be nice to be able to code (part
of) that in 3.x.

Benoit

2009/11/16 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> 2009/11/16 Gael Guennebaud <gael.guennebaud@xxxxxxxxx>:
>>
>>
>> On Sun, Nov 15, 2009 at 6:29 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
>> wrote:
>>>
>>> 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?
>>
>> well, for many kind of physic problems where it is advantageous to use a
>> band storage rather than a general sparse storage.
>
> Well ok, I'll just trust you that there exist such problems, in which
> one actually needs to do matrix products... so obviously i don't mind
> if you want to work on this, but personnally i still need motivational
> examples if you would like me to work on this myself. In other words,
> i still find it very tempting to say that special matrix classes are
> only for special uses like solving.
>
>>
>> In order to factorize the code implementing the operators and functions (+,
>> *, .transpose(), etc.) over all the different matrix type, I see three
>> options:
>>
>> 1 - Move everything to AnyMatrixBase and find a nice way to generate
>> explicit compilation errors when one operation is not supported by the given
>> storage type.
>>
>> 2 - Split them into a set of "base" matrix classes, one per type of
>> expression (e.g., CwiseUnaryOps, TransposeOps, etc.) and then let the main
>> matrix base class of each kind of storage inherits the set of supported base
>> operations. So the current MatrixBase class for dense storage would inherit
>> CwiseUnaryOps, TransposeOps, CwiseBinaryOps, ProductOps, etc. One issue is
>> that we cannot use multiple inheritence because:
>>  a) MSVC does not like multiple inheritence of empty classes, it generates
>> additional storage cost.
>>  b) operator - and * would be inplemented by various base classes.
>> To overcome this issue, one common meta programming solution is to use a
>> derivation chain, e.g.:
>>
>> template<typename Derived> class MatrixBase : public
>> DerivationChain<Derived,CwiseUnaryOps,TransposeOps, CwiseBinaryOps,
>> ProductOps> {};
>>
>> The DerivationChain helper class will generate that inheritence tree:
>>
>> MatrixBase<Derived> -> CwiseUnaryOps<Derived> -> TransposeOps<Derived> ->
>> CwiseBinaryOps<Derived> -> ProductOps<Derived> -> AnyMatrixBase<Derived>
>>
>> 3 - the third option is to use the preprocessor, not macro, but rather one
>> #include per set of operations, e.g.:
>>
>> template<typename Derived> class MatrixBase
>> {
>>   // ...
>>   #include <CwiseUnaryOps>
>>   #include <TransposeOps>
>>   #include <CwiseBinaryOps>
>>   #include <ProductOps>
>> };
>>
>>
>> Because of the difficulties to generate good documentation and good
>> compilation errors I don't really like solution 1.
>>
>> The solution 2 is the C++ solution. It also splits the documentation which
>> could be an advantage (or not ?). Still regarding the documentation, we
>> could cheat a bit by showing a flat inheritence tree hidding the derivation
>> chain stuff.
>>
>> Finally, the solution 3 is much easier, both for us and for the compiler,
>> but I'm not sure how doxygen handle that. I'll try and let you know.
>>
>
> If solution 3 is easier for both us and the compiler, that dwarves any
> drawbacks it might have for doxygen...
>
> There's something i'm not sure to understand: this solution only takes
> care of defining operators/methods as just returning xpr objects,
> right? The actual implementation of the xprs can't be shared at all
> across base matrix types, as far as I can see.
>
> Then frankly, I'd say don't worry too much about this aspect (of the
> operators returning xprs) because that is only a relatively small part
> of the code and effort... anyhow, the bulk of the work will be to
> write all the xpr classes.
>
> Benoit
>
>> gael.
>>
>>>
>>> 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
>>> >> >>
>>> >> >>
>>> >> >
>>> >> >
>>> >> >
>>> >> >
>>> >> >
>>> >>
>>> >>
>>> >
>>> >
>>> >
>>> >
>>>
>>>
>>
>>
>>
>>
>