[eigen] slow sparse dense product

[Sebastian Birk <birk@xxxxxxxxxxxxxxxxxxxxx>]

Hi,

again I encountered an unexpected slow product of a sparse matrix A and
a dense matrix X. As can be seen in the attached file, computing A*X or
A^H*X is about 50% slower for X being col-major than X being row-major
and even computing A*c for all columns c of X. Obviously, the col-major
storage order is non-optimal, but 50% slower seems to be too much.

Sebastian

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

#include <boost/timer.hpp>

using namespace std;
using namespace Eigen;

int main()
{
    const int size = 100000;
    const int nnzpc = 90;
    const int iter = 30;
    const int cols = 4;
    boost::timer time;

    MatrixXcd X(size,cols);
    static MatrixXcd Y(size,cols);
    Matrix<complex<double>,Dynamic,Dynamic,RowMajor> Xrm(size,cols);
    static Matrix<complex<double>,Dynamic,Dynamic,RowMajor> Yrm(size,cols);
    SparseMatrix<complex<double>,ColMajor> A(size,size);

    A.reserve(VectorXi::Constant(size,nnzpc));
    for (int n = 0; n < size; ++n)
        for (int m = 0; m < nnzpc; ++m)
            A.insert((m*nnzpc+n)%size,n) = complex<double>(m,n);
    A.makeCompressed();

    time.restart();
    for (int i = 0; i < iter; ++i) {
        Y.noalias() = A*X;
        Y.noalias() = A.adjoint()*X;
    }
    cout << "time Y=AX,Y=A^HX: " << time.elapsed() << endl;

    time.restart();
    for (int i = 0; i < iter; ++i) {
        Yrm.noalias() = A*Xrm;
        Yrm.noalias() = A.adjoint()*Xrm;
    }
    cout << "time Y=AX,Y=A^HX (row-major-matrix): " << time.elapsed() << endl;

    time.restart();
    for (int i = 0; i < iter; ++i)
        for (int c = 0; c < cols; c++) {
            Y.col(c).noalias() = A*X.col(c);
            Y.col(c).noalias() = A.adjoint()*X.col(c);
        }
    cout << "time Y=AX,Y=A^HX (col-wise): " << time.elapsed() << endl;
}

//time Y=AX,Y=A^HX: 13.59
//time Y=AX,Y=A^HX (row-major-matrix): 9.43
//time Y=AX,Y=A^HX (col-wise): 9.27