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 16:45:24 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

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

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

**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] Arbitrary precision arithmetic support in Eigen** - Next by Date:
**[eigen] Proposal: Matrices requiring compact storage of triangular part** - Previous by thread:
**Re: [eigen] Implementing a custom-matrix product** - Next by thread:
**Re: [eigen] Implementing a custom-matrix product**

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