Bugzilla – Bug 553
Matrix-matrix product with sparse self-adjoint matrix doesn't compile
Last modified: 2013-06-28 23:00:32 UTC
I have a self-adjoint sparse matrix and want to multiply it with another sparse matrix. Compiling the following code produces a compiler error on gcc 4.5.4:
using namespace Eigen;
z = x * y; // works
z = x.selfadjointView<Lower>() * y; // compiler error
The error is the following:
testspmatrix.cpp:14:38: error: no match for 'operator*' in 'Eigen::SparseMatrixBase<Derived>::selfadjointView() [with unsigned int UpLo = 1u, Derived = Eigen::SparseMatrix<double>]() * y'
Eigen/src/SparseCore/../plugins/CommonCwiseUnaryOps.h:80:1: note: candidates are: const Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple2_op<double, std::complex<double> >, const Eigen::SparseMatrix<double> > Eigen::operator*(const std::complex<double>&, const Eigen::SparseMatrixBase<Eigen::SparseMatrix<double> >::StorageBaseType&)
Eigen/src/SparseCore/../plugins/CommonCwiseUnaryOps.h:76:1: note: const Eigen::SparseMatrixBase<Eigen::SparseMatrix<double> >::ScalarMultipleReturnType Eigen::operator*(const Eigen::SparseMatrixBase<Eigen::SparseMatrix<double, 0, int> >::Scalar&, const Eigen::SparseMatrixBase<Eigen::SparseMatrix<double> >::StorageBaseType&)
This seems like very basic functionality to have and is an absolute blocker for me. Am I doing something wrong?
"very basic functionality" -> I'm sure as I'm not aware of any lib implementing this. Actually, the only way I think of to implement this would be to first evaluate x.selfadjointView<Lower>() into a full SparseMatrix:
x2 = x.selfadjointView<Lower>();
z = x2 * y;
(needs the devel branch to work)
Nevertheless, the API could support that.
That's what I ended up doing in my code. If you tell me that there is no better way to implement this, then I guess this can be closed.
I would maybe suggest, if possible, to add one of your UPPERCASE_ERROR_MESSAGES for this use case, as well as to update the documentation for sparse matrices (which is right now very, uh, *sparse*) to explicitly state which operations are supported and which not. Otherwise it's very hard to figure out what's going on.
I forgot to add that I think the API should *not* create the temporary matrix transparently, because I think that would be a huge performance trap. Better tell people directly that this doesn't work. Especially in tight loops you usually want to avoid any unnecessary allocations.
sparse-sparse products already have to introduce temporaries unless the two factors and result matrices have all the same storage order, so there is no reason for not supporting sparse selfadjoint time sparse products. The introduction of temporaries is documented.
Date: 2013-06-28 22:27:45
Summary: Fiw bug 553: add support for sparse matrix time sparse self-adjoint view products
Date: 2013-06-28 22:56:26
Summary: Add missing sparse matrix constructor from sparse self-adjoint views, and add documentation for sparse time selfadjoint matrix