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.
-- GitLab Migration Automatic Message --
This bug has been migrated to gitlab.com's GitLab instance and has been closed from further activity.
You can subscribe and participate further through the new bug through this link to our GitLab instance: https://gitlab.com/libeigen/eigen/issues/574.