Re: [eigen] Implementing a custom-matrix product

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

Hi Gael,

   I ran into some problems integrating my block arithmetic into the Eigen3 
framework. I tried following the SelfadjointMatrixMatrix example.

1. My product wrapper class is named SymmetricMatrixMatrix (for now). Should 
this inherit ProductBase like SelfadjointMatrixMatrix? In doing so, I get an 
assertion failure in the ProductBase CTOR because it is expecting the nested 
storage matrices to be dimensionally multiplication-compatible whereas the 
RFPF format violates this because of re-packaging the layout.


2. I tried looking at DiagonalProduct. This inherits MatrixBase and declares 
EIGEN_DENSE_PUBLIC_INTERFACE. Is this the route I could follow?

3. I noticed there is a EIGEN_GENERIC_PUBLIC_INTERFACE also. Should I just not 
inherit from any previous baseclass and implement the multiplier as a 
top-level class? In this case, I don't see an function-call chain from 
DenseBase<>::operator=() to SymmetricMatrixMatrix::scaleAndAddTo(). I am 
guessing the various base-classes in Eigen take care of bridging the 
expressions to their evaluations.

thanks,
Manoj


On Wednesday 14 July 2010 04:52:18 am Gael Guennebaud wrote:
> Hi Manoj,
>
> this is simple. The first thing to do is transform high level code like:
>
> C += A * B;
>
> where A is a CompactSymmetricMatrix, to something like:
>
> compact_symm(A, B, C);
>
> so that we have both the product operands and the destination matrix at
> hand.
>
> Then, as you said, implementing compact_symm is a matter of three
> lines of code, maybe something like that:
> C.block() += A.internalStorage().block().triangularView<..>() * B.block();
> C.block() += A.internalStorage().block().triangularView<..>().transpose()
> * B.block();
> C.block() += A.internalStorage().block() * B.block();
>
> where A.internalStorage() corresponds to dense Matrix used to store
> the triangular part in a compact way. (internalStorage is not a good
> name, but anyway...)
>
> So it remains to details the first step. For this I suggest you to
> mimic what is done for the SelfadjointView matrix product:
> file src/Core/SelfadjointView.h, lines 111,119
> file src/Core/product/SelfadjointMatrixMatrix.h, lines 385 and below...
>
> Actually this is example is much more complicated that what you have
> to achieve. Typically, you can implement the the above three lines of
> code directly in the scaleAndAddTo() method instead of calling another
> function, you don't have to distinguish between matrix and vector,
> etc.
>
> Are you willing (able) to share your code? I'm asking because this is
> definitively something we'll add in eigen eventually, when all other
> high priority tasks will be completed ;)
>
> good luck,
>
> gael
>
> On Wed, Jul 14, 2010 at 12:26 AM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> 
wrote:
> > Hi,
> >
> >   I have a basic implementation of a SymmetricMatrix class that stores
> > only the lower triangle in Rectangular Full Packed Format (RFPF) from a
> > LAWN reference that Gael mentioned a little while ago.
> >
> >   I would like to implement a SymmetricMatrix * Matrix product. How do I
> > go about it? Here are my preliminary issues:
> >
> > - classes like DiagonalWrapper and TriangularView return a custom class
> > like DiagonalProduct and TriangularProduct when the pre-multiplication
> > operator * (member function) is called. But then, someone mentioned
> > ReturnByValue in a different context but that also seems to serve a
> > similar purpose. I am a little confused. How do these concepts relate?
> > Are they mutually exclusive?
> >
> > - I am not able to identify which member of the dense matrix class
> > hierarchy actually performs the evaluation of the product into a matrix
> > with dense-storage. Even DenseStorageBase seems to be passing this off to
> > its base class operator=() or some _set() method.
> >
> > - There are coeff() methods that return Scalar, and there are packet()
> > methods that return PacketScalar. I am guessing the latter is used when
> > the expression is vectorizable. How and where is this decided?
> >
> > - The product that I am interested in shouldn't have to worry about
> > PacketScalar etc. because it reduces to a sequence of high-level Eigen
> > TRMM and GEMM calls. This is due to the layout - the triangular part is
> > stored as two sub-triangles (one is adjoint-ed) and one rectangular
> > block, all packed so that they form a rectangular array of total size
> > N*(N+1)/2. When the product evaluation is triggered, the most efficient
> > operation would be with these triangular/rectangular blocks as opposed to
> > coeff() or packet() access which would be expensive. How can this be
> > ensured?
> >
> > Thanks,
> > Manoj