Bug 553 - Matrix-matrix product with sparse self-adjoint matrix doesn't compile
Matrix-matrix product with sparse self-adjoint matrix doesn't compile
Status: RESOLVED FIXED
Product: Eigen
Classification: Unclassified
Component: Sparse
unspecified
All All
: Normal major
Assigned To: Desire NUENTSA
:
Depends on:
Blocks: 3.2
  Show dependency treegraph
 
Reported: 2013-02-14 15:59 UTC by Clemens Hofreither
Modified: 2013-06-28 23:00 UTC (History)
1 user (show)



Attachments

Description Clemens Hofreither 2013-02-14 15:59:50 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:

-------------------------------------

#include <Eigen/Core>
#include <Eigen/Sparse>

using namespace Eigen;

int main()
{
    SparseMatrix<double> x(10,10);
    SparseMatrix<double> y(10,10);
    SparseMatrix<double> z(10,10);

    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?
Comment 1 Gael Guennebaud 2013-02-25 19:03:02 UTC
"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:

SparseMatrix<double>x2;
x2 = x.selfadjointView<Lower>();
z = x2 * y;

(needs the devel branch to work)

Nevertheless, the API could support that.
Comment 2 Clemens Hofreither 2013-02-26 10:16:05 UTC
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.
Comment 3 Clemens Hofreither 2013-02-26 14:39:35 UTC
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.
Comment 4 Gael Guennebaud 2013-06-28 23:00:32 UTC
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.

https://bitbucket.org/eigen/eigen/commits/39429a449b68/
Changeset:   39429a449b68
User:        dnuentsa
Date:        2013-06-28 22:27:45
Summary:     Fiw bug 553: add support for sparse matrix time sparse self-adjoint view products

https://bitbucket.org/eigen/eigen/commits/743776e65442/
Changeset:   743776e65442
User:        ggael
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

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