Re: [eigen] Implementing a custom-matrix product

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

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.

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