[eigen] Recursive Matrices in Eigen

[Darius Rueckert <darius.rueckert@xxxxxx>]

Hello everyone,

for the past months we have been working on making recursive matrices 
work in eigen. Here is an example:

SparseMatrix<Matrix<double,4,3>> A;
Matrix<Matrix<double,4,1>,-1,1> x,b;
x = A * b;

The obvious advantage is that we get support for block matrices on all 
eigen types for "free" ( DiagonalMatrix,TriangularView,LDLT,...). After 
some testing, the performance appeared to be "better than expected" and 
we have decided to write a paper on how to apply recursive matrices on 
non-linear least squares problems. The paper is accepted and will be 
published at HPG in July. Here is the current version:

https://cloud9.cs.fau.de/index.php/s/UINyrgAUtaAiWvR


Unfortunately there are two limitations which I was not able to fix.

1. All inner matrix types require an additional wrapper class, which I 
have called "MatrixScalar". The matrix A from above is then:

SparseMatrix<MatrixScalar<Matrix<double,4,3>>> A;

2. Some eigen kernels do not work directly. For example, some operations 
use the function scaleAndAddTo(..., Scalar(1)), which does not work if 
"Scalar" is a recursive matrix type. To fix this I had to specialize a 
lot of eigen kernels, which creates redundancy with only minor changes.


So here are my questions to you:

What do you think about this idea?

Do you think a deep integration into eigen would be possible? (Without 
recursive specializations and without "MatrixScalar" class)

The code can be found here:
https://github.com/darglein/EigenRecursive

Example for sparse block-matrices:
https://github.com/darglein/EigenRecursive/blob/master/samples/sparseOperators/sparseOperators_main.cpp

With best regards,
Darius

--
Darius Rueckert
Informatik 9 (Graphische Datenverarbeitung) http://lgdv.tf.fau.de
Universitaet Erlangen-Nuernberg, Cauerstraße 11, 91058 Erlangen