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 15:22:01 -0400*Organization*: EECS Dept., University of Michigan, Ann Arbor, MI, USA

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. 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. 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:*Gael Guennebaud

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Arbitrary precision arithmetic support in Eigen** - Next by Date:
**Re: [eigen] Arbitrary precision arithmetic support in Eigen** - 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/ |