[eigen] about band/skyline matrices

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi all,

I've started to implement a generic BandMatrix class with storage which aims to efficiently and compactly represents matrices where only a few sub- and super-diagonals near the main diagonal are non zero. See:
http://bitbucket.org/bjacob/eigen2_special_matrices/src/tip/Eigen/src/Core/BandMatrix.h

My question is about the storage scheme. Currently I've simply packed the diagonal in a single vector, e.g.:

[11 12 13 00]
[21 22 23 24]
[00 32 33 34] 
[00 00 43 44]

is stored:

size = 4
nb supers = 2
nb subs = 1
data = "" 22 33 44 12 23 34 13 24 21 32 43]

The main advantages are:
 - the diagonal elements are sequentially stored in memory
 - minimal memory consumption
 - it is easy to extract any diagonals
Drawbacks:
 - it is not possible to extract a row or a column
 - => it is not possible to implement, e.g., an efficient matrix product
 - limited to square matrices


So another possibility would be to implement the BLAS storage scheme where the key idea is to "move" the columns up and down such that the diagonals become rows, and then you prune the empty rows. The previous matrix would be stored in a column major matrix like this:

[__ __ 13 24]
[__ 12 23 34]
[11 22 33 44]
[21 32 43 __]


The main advantages are:
 - it is easy to extract the non zero part of any column,
 - the columns are sequentially stored in memory => products can be efficiently implemented
 - it is still possible to extract any row using something similar that DiagonalCoeffs
 - the diagonals are still easy to extract
 - compatible with BLAS/LAPACK
Drawbacks:
 - a little waste of memory (negligible if the number of sub and supers diagonal is small compared to the size of the matrix)
 - the diagonals are not sequentially stored in memory that is not very good as a representation for tridiagonal matrices.

Note that the last drawback can be overcame using a row major matrix. In that case the storage become compatible with the BLAS triadiagonal storage scheme but not with the band storage scheme which impose a colmajor matrix.

So here is what I suggest:

- Implement BandMatrix using the BLAS/LAPACK storage scheme with the additional possiblity to use a rowmajor matrix for the underlying storage
- Add a TridiagonalMatrix<Scalar,Size,Options> class inheriting with simplified template parameters and simplified ctor, e.g.:

template<typename Scalar, int Size, int Supers, int Subs, int Options = ColMajor>
class BandMatrix;

Options = ColMajor or RowMajor | SelfAdjoint

template<typename Scalar, int Size, int Options>
class TridiagonalMatrix : public BandMatrix<Scalar,Size,1,1,Options | RowMajor>;

Another possiblity would be to automatically fall back to a RowMajor storage if Supers<=1 and Subs<=1 and the storage order is left unspecified by the user, e.g:

BandMatrix<float,Dynamic,1,1>               => creates a BLAS compatible tridiagonal matrix
BandMatrix<float,Dynamic,1,1,ColMajor> => creates a BLAS compatible band matrix


Any opinions ?

cheers,
Gael.