MatrixFree-Matrix (and Vector) product with Eigen Library

Question

I am creating a C++ software for which I need a wrapper, based on the Eigen library, implementing the operator* similar to the one explained in this official web-page

https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html

With the code from the aforementioned web-page I am able to wrap my templated class in the MatrixReplacement one.

Keeping the implementation of MatrixReplacement from the example, the following (main) code works

int main()
{
    int n = 10;
    Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
    S = S.transpose()*S;
    MatrixReplacement A;
    A.attachMyMatrix(S);
    Eigen::VectorXd b(n,1), x1(n,1), x2(n,1);
    b.setRandom();

    x1.noalias() = A * b;
    x2.noalias() = S * b;
    std::cout<<(x1-x2).colwise().norm()<<std::endl;
}

But if, instead of using vectors, I want to use matrices for b and x the code will not compile complaining for missing members and types

int main()
{
    int n = 10;
    Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
    S = S.transpose()*S;
    MatrixReplacement A;
    A.attachMyMatrix(S);
    Eigen::MatrixXd b(n,3), x1(n,3), x2(n,3);
    b.setRandom();

    x1.noalias() = A * b;   //<<<<<<<<<<< COMPILE ERROR
    x2.noalias() = S * b;
    std::cout<<(x1-x2).colwise().norm()<<std::endl;
}

MY QUESTION IS: what is missing from the example at the web-page https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html in order to work with my second main?

Thanks!

---

EDIT: in the example in the web-page the for loop

for(Index i=0; i<lhs.cols(); ++i)
    dst += rhs(i) * lhs.my_matrix().col(i); 

will need to be changed to something like

for(Index i=0; i<lhs.cols(); ++i)
    for(Index j=0; j<rhs.cols(); ++j)
        dst.col(j) += rhs(i,j) * lhs.my_matrix().col(i);

or simply

dst.noalias() += lhs.my_matrix() * rhs

Answer 1

In the specialization of internal::generic_product_impl, you need to change GemvProduct to GemmProduct.