Re: [eigen] Implementing a custom-matrix product |

[ Thread Index | Date 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 19:43:47 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

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

**References**:**[eigen] Implementing a custom-matrix product***From:*Manoj Rajagopalan

**Re: [eigen] Implementing a custom-matrix product***From:*Manoj Rajagopalan

**Re: [eigen] Implementing a custom-matrix product***From:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness** - Next by Date:
**Re: [eigen] discrepancy in (triangular view + transpose) and self-adjointness** - Previous by thread:
**Re: [eigen] Implementing a custom-matrix product** - Next by thread:
**[eigen] Arbitrary precision arithmetic support in Eigen**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |