Re: [eigen] Implementing a custom-matrix product

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

Success!

I now have the container SymmetricMatrix<Lower> and the product 
SymmetricMatrixMatrix<Lower> working. On to the solver ... but that's a 
different thread! :-)

thanks, Gael and others.
-- Manoj


On Wednesday 14 July 2010 05:39:22 pm Gael Guennebaud wrote:
> 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