Re: [eigen] General design issue

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

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

oh, and that is exactly what your e-mail was about, right? that's why
you didn't start discussing the implementation of xprs ?

Benoit


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