Re: [eigen] Implementing a custom-matrix product

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Wed, Jul 14, 2010 at 9:22 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:
>
> 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.

Ideally symm * mat should return a
SymmetricMatrixMatrix<SymmetricMatrix,Matrix> that shoud inherits
ProductBase<SymmetricMatrixMatrix<..>,SymmetricMatrix,Matrix>, then
the sizes should be compatible.

(The SymmetricMatrix::rows() and ::cols() functions have to return the
size of the virtual matrix, not the ones of the underlying compact
object.)


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

that's not a good idea.

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

The scaleAndAddTo() mechanism is specific to ProductBase and is allows
to deal with expression like:

C = s*(A*B);

without evaluating A*B into a temporary. The scalar factor given by
scaleAndAddTo() correspond to s.

I don't see any good reason why inheriting ProductBase would not work.
Actually our sparse*dense products also inherit ProductBase, so you
could have a look at the file src/Sparse/SparseDenseProduct.h (line
148) for a simpler example which should be very close to your use
case.

gael

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