[eigen] Proposal: Matrices requiring compact storage of triangular part

[Manoj Rajagopalan <rmanoj@xxxxxxxxx>]

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 :-)


Thanks,
Manoj