|Re: [eigen] Proposal: Matrices requiring compact storage of triangular part|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
- To: eigen@xxxxxxxxxxxxxxxxxxx
- Subject: Re: [eigen] Proposal: Matrices requiring compact storage of triangular part
- From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
- Date: Mon, 19 Jul 2010 17:03:44 +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=ZqG0cCohEXITm4OCwIX9wVE7RpW6dOdEArtvs+E+jVo=; b=Y976Le7KVQfKzICzdqiY3PFFCvSNZQv2Ml/wvfDv1cCSOC+EfFiKyUriMeZv+BCnBs WZZ40OAezWraIY2yFf40BRochqsQ+gWg7ToEu3yYagKQP63N99jzcmXMw7VEG3mzyfHO pgfV8y5waqnjrpE8ZNlycy9w62WoH4MU5anWE=
- 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=n61VkroP3wAueVc1PGf2p1PmZiXP66ZAJm3HHpqNeVM37ZTCb7UIZ4dheOLx7hbkyk onfTvZ3oHWwypKKyryVj1helNgxI2xnrRQIq0S4hBvijLYqlhoe9IZ46WcQ1m1FPAx0i VOF9ISjxOYmyr8AaMZHrNz2FpgRQpz8leMzKg=
thanks for this detailed email, but as Benoit said this is a bit too
much for now, and most of it can be easily implemented without any
modifications in Eigen. Actually, I think there are only two welcome
changes that are affecting Core:
- rename the current TriangularBase to TriangularViewBase
- add support for symmetric complex matrix (not only selfadjoint).
For the rest, this can be easily implemented in a separated module
providing compact triangular representations and related algorithms.
Also, regarding the design, I'd rather see something like:
rather than a selfadjoint/symmetric view on a CompactTriangularMatrix.
Also, since CompactSelfAdjointMatrix and CompactSymmetricMatrix are
very very similar, they could either be factorized in a unique class
with a template option (but I cannot find a good name for either
symmetric or selfadjoint), or they could inherit a common base class
factoring most of the code.
LAPACK original triangular layout has zero interest IMO, so let's
forget about it.
that was just some quick thoughts ;)
On Wed, Jul 14, 2010 at 11:00 PM, Manoj Rajagopalan <rmanoj@xxxxxxxxx> wrote:
> Hi Eigen developers,
> Gael's recent reference to the RFPF storage layout for triangular (portions
> of) matrices got me experimenting with implementing a SymmetricMatrix class
> and related products. The question of sharing the code with Eigen came up..
> Here are my thoughts on slight changes to Eigen3 that might be required
> before the integration can be performed.
> 1. Define a Triangular_RFPF class (for lack of a better name). I don't know if
> this can be folded into the existing TriangularBase hierarchy (which assumes
> rectangular/dense 'models' being viewed). If not, then I suggest renaming the
> existing TriangularBase to TriangularViewBase and the introducing a
> TriangularBase as the CRTP base-class for all types of triangle-based views.
> 2. Define a selfadjointView() method that returns a hermitian view of the
> stored triangle. Since the Triangular_RFPF class itself will take a Lower|
> Upper mode indicator, this method won't be templated like it is for
> MatrixBase<> derivatives.
> 3. Also define a symmetricView() method that returns a symmetric view of the
> stored triangle.
> 4. The two views will be identical for real Scalar types. For complex the two
> are mutually-exclusive. So we will need to introduce a Symmetric=0x20 into
> the mode enum in util/Constants.h. For real Scalar cases, Flags &
> (SelfAdjoint|Symmetric) should be true. In other words, both bits should be
> set. For the complex case only one should.
> 5. Some kind of ei_meta_if<> will need to be used to return identical wrapper
> classes for symmetricView() and selfadjointView() if Scalar is real, so that
> existing product implementations etc. can work as is. Then, only symmetric
> products will have to be implemented for complex Scalars.
> 6. New product wrappers will have to be written when Triangular_RFPF or its
> (selfadjoint|symmetric)View() are involved but this will involve working with
> blocks of the RFPF storage matrix as discussed in the original paper so this
> part should be easy.
> 7. New triangular solver wrappers will have to be written but this should
> again be easy - the paper discusses the lower-triangular case.
> 8. The SelfAdjointView<> wrapper serves as a validator for LLT and LDLT. So
> symmetricView() for real Scalars should have a SymmetricViewReturnType
> typedef that maps to SelfAdjointView<> for real scalars. This will ensure
> minimal changes to Eigen3. Since there are no specializations of L(D)LT
> implementations for classes other than SelfAdjointView, we are good.
> 9. Since TriangularView<>, Triangular_RFPF<> and a hypothetical vanilla
> TriangularMatrix class (with LAPACK layout) share a common abstraction, this
> begs the question of whether they should inherit from siblings at some level
> since they differ either in ownership of the data or in the layout but are
> all triangle-accessors. In fact, layout and ownership are orthogonal (you can
> have a view with RFPF layout instead of the vanilla LAPACK layout), there is
> scope for more generalization but this might be too much work.
> 10. Complex symmetric matrices arise in certain differential equations
> situations and there are efficient solvers for certain classes that use
> complex orthogonal (as opposed to unitary) transformations. Having complex
> symmetric classes would be desirable at some point.
> 11. I think the complex case can be safely postponed to a later release of
> Eigen depending on use-community's interest. Then SymmetricProduct wrappers
> can be deferred. Common cases will work with the SelfAdjointProduct wrappers
> and Eigen releases won't be stalled :-)