Re: [eigen] Implementing a custom-matrix product

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Wed, Jul 14, 2010 at 10:45 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:
>
> BTW, what value for Flags should the SymmetricMatrix container and
> SymmetricMatrixMatrix product set? Looking at all the options in
> util/Constants.h, nothing seemed to fit the bill so I'm setting Flags to 0.

for SymmetricMatrixMatrix simply use the default one set for
ProductBase (EvalBeforeNestingBit and EvalBeforeAssignBit are the ones
which matter here).

Regarding SymmetricMatrix, indeed 0 seems to be the most appropriate.

gael

>
> I'm trying out a few more things now ... will get back to you.
>
> And I forgot to answer the last question - yes, I'll be happy to share the
> code for incorporation into Eigen3 when it is mature. But before that I think
> a slight change to Eigen3 design might become necessary. Will start a
> separate thread on that.
>
> Thanks,
> Manoj
>
>
> On Wednesday 14 July 2010 03:42:00 pm Gael Guennebaud wrote:
>> 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
>
>
>
>