New user self-registration is disabled due to spam. Please email eigen-core-team @ lists.tuxfamily.org if you need an account.
Before reporting a bug, please make sure that your Eigen version is up-to-date!
Bug 574 - Support for Matrices consisting of fixed-size Blocks
Summary: Support for Matrices consisting of fixed-size Blocks
Status: NEW
Alias: None
Product: Eigen
Classification: Unclassified
Component: Core - general (show other bugs)
Version: 3.2
Hardware: All All
: Normal enhancement
Assignee: Nobody
URL:
Whiteboard:
Keywords:
Depends on:
Blocks:
 
Reported: 2013-03-23 11:37 UTC by Daniel Vollmer
Modified: 2013-04-08 08:57 UTC (History)
3 users (show)



Attachments
Example implementation (does not quite compile) (2.64 KB, application/zip)
2013-03-23 11:37 UTC, Daniel Vollmer
no flags Details

Description Daniel Vollmer 2013-03-23 11:37:58 UTC
Created attachment 320 [details]
Example implementation (does not quite compile)

It would be quite useful to define Eigen-matrices (of all kinds, dense, sparse, diagonal). I am not expecting Eigen's built-in decompositions and advanced to work immediately with such an entity, but supporting the basic operations should be possible.

One approach demonstrated in the attached example, is to treat a fixed size BlockMatrix (e.g. Eigen::Matrix<double, 5, 5>) similarly to a complex number. When one takes care to define the proper scalar_product_traits, both matrix-matrix and matrix-vector products work.
Two problems with the example implementation are:
1) The existing product implementations declare a pre-multiplier alpha of type "const typename Res::Scalar& alpha". This type declaration leads to problems when LHS is SparseMatrix<Matrix<double, 5, 5> > and RHS is
Vector<Matrix<double, 5, 1> >. Then (via scalar_product_traits) the product result scalar is again Vector<Matrix<double, 5, 1> >, which means we try to multiply Matrix<double, 5, 1> * Matrix<double, 5, 1>. In this case alpha really should be a scalar, but another typedef Matrix::DefinitelyForSureScalar is also very ugly. I'm not sure how alpha should be declared to make this work (without introducing *another* typedef to all Matrix-types which I'm not sure is the right approach).

2) Many temporary "Scalar"s of the matrix's type (i.e. variables of type Matrix<double, 5, 5>) are created via their single-argument constructor and initialized for little benefit, for example SparseMatrix::insertUncompressed, insertBackByOuter, inside AmbiVector, etc. and then overwritten be the actual value. This is cheap for double, but costly for a heavier weight scalar-type. When one wants the additive or multiplicative identity for a custom type, one idea could be to not treat this as a Scalar but maybe CwiseNullaryOp<> that could be implemented efficiently and correctly.


Another approach is to define an Adaptor-wrapper (or maybe a "View") that makes an 
Matrix<Matrix<double, 5, 5>, M, N> look like Matrix<double, M*5, N*5>. This increases interoperability with everything else in Eigen that relies on Scalar being a real scalar, but would probably lose all advantages of the block-structure (i.e. the ability to work on tightly-packed sub-blocks).

Maybe both are worthwhile?

Relevant mailing list thread:
http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2013/03/msg00072.html
Possibly relevant / related / helpful bug:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=279 (Be more flexible regarding mixing types)
Comment 1 Daniel Vollmer 2013-03-23 11:40:28 UTC
Sorry, I messed up the introductory paragraph:
It would be quite useful to be able define Eigen-matrices (of all kinds, dense, sparse,
diagonal) to have a fixed-size matrix as their "Scalar" type.

Note You need to log in before you can comment on or make changes to this bug.