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:
Possibly relevant / related / helpful bug:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=279 (Be more flexible regarding mixing types)
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.