Hi all,

I'd like to multiply a sparse matrix with a (row of) another sparse matrix. Doing that directly does not work currently:

Code: Select all

Eigen::SparseMatrix A(n,n);
Eigen::SparseMatrix b(n,1);
Eigen::VectorX x;
x = A*b;


Using an intermediate sparse matrix to store A*x does work, but that requires memory allocations on each an every call.
Currently, I have a work around that looks like:

Code: Select all

template<typename DerivedA, int Order, typename DerivedB>
void prod(Eigen::SparseMatrix<DerivedA, Order>& A,
                     typename Eigen::SparseMatrix<DerivedB>::InnerIterator& b,
                     Eigen::VectorXi& x
                     )
{
    x.resize(A.rows());
    x.setConstant(0);
    for(; b; ++b)
    {
        auto Ac = A.col(b.index());
        if(Ac.nonZeros())
            x += Ac * b.value();
    }
}


Could someone give some hints on how to make this a little more robust ?
Is it possible to return x, instead of accepting it as input, without re-allocating on return ?

Side-note: If A becomes a const input, this code will endlessly loop in .nonZeros(), since that calls itself.
Is that a bug, or just not supported ?

What is even more sad is that when computing A * B with A and B sparse, the intermediate results of each A * B(:,j) is computed in a dense vector. So if the user evaluate A * B into a dense matrix/vector C then, ideally we would pass C to the sparse product kernel to directly use it. This would require some little code refactoring.

Side-note: If A becomes a const input, this code will endlessly loop in .nonZeros(), since that calls itself.
Is that a bug, or just not supported ?

https://bitbucket.org/eigen/eigen/commits/397aae7aa8d7/
Changeset: 397aae7aa8d7
User: ggael
Date: 2013-11-10 15:26:07
Summary: Add missing nonZeros() overload in Block<SparseMatrixBase<>>

https://bitbucket.org/eigen/eigen/commits/5454b177f8b9/
Changeset: 5454b177f8b9
User: ggael
Date: 2013-11-10 16:16:50
Summary: Use the specialization of Block<SparseMatrix> for const matrices too

ggael wrote:This would require some little code refactoring.

I wouldn't mind taking a shot at that, although I'm not sure my c++-template wizardry-skills are up to it.

Do you have some pointers/hints/ideas on how to approach that ?