|Re: [eigen] Implementing a custom-matrix product|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Implementing a custom-matrix product
- From: Manoj Rajagopalan <rmanoj@xxxxxxxxx>
- Date: Wed, 14 Jul 2010 15:22:01 -0400
- Organization: EECS Dept., University of Michigan, Ann Arbor, MI, USA
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.
2. I tried looking at DiagonalProduct. This inherits MatrixBase and declares
EIGEN_DENSE_PUBLIC_INTERFACE. Is this the route I could follow?
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.
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
> 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,
> 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,
> On Wed, Jul 14, 2010 at 12:26 AM, Manoj Rajagopalan <rmanoj@xxxxxxxxx>
> > 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