|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: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Wed, 14 Jul 2010 21:42:00 +0200
- Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:received:mime-version:received:in-reply-to :references:from:date:message-id:subject:to:content-type :content-transfer-encoding; bh=XaDYoAog6Y2GaNjvAUom7p3iCV+qJbVV9Rzea8TgkmY=; b=oyWj+b8IOe8KgnrZP6staBJy33V3TbX6b512x3gaidf/ET0rlWlInxBwac5horo72m WxAfdJzCixPaX6aLC2uWZRgOXRASH73amjDKcPgZIyKa+Cvnh4FFvlwtQj7AqSdS2c8j tBEh2w4yG9kpBtwn/k3dMP+5ePS96K/xWXy10=
- Domainkey-signature: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type:content-transfer-encoding; b=X3DQdAH8h/wh4o1jHnG0AMTirdaRGz5R4HG75q+7ojjfa7h5JSW5uKlGytmro+2H5K RslvGcWy+X7Zsucptz4G92OUSZyRuEGg4YVorM/HX4nKDuq8//pMkCgHPu0kPTP+xMeG O/7RV3a+XtLL3QSVC+LTqcizhGVwD/wM632/U=
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
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
> 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
> 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