Re: [eigen] Implementing a custom-matrix product

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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