使用Eigen库进行MatrixFree矩阵（和向量）乘法

Question

我正在创建一个C++软件，需要一个基于Eigen库的包装器，实现类似于官方网页中解释的operator*操作符。

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

使用上述网页的代码，我能够将我的模板类包装在MatrixReplacement中。

保留示例中MatrixReplacement的实现，以下（主要）代码可用：

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;
}

但是，如果我想使用矩阵作为b和x的数据类型，而不是向量，代码将无法编译，因为缺少成员和类型。

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;
}

我的问题是：在 https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html 的示例中缺少什么，以便与我的第二个主要程序一起使用？

谢谢！

---

编辑： 在网页示例中的for循环

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

将需要更改为类似于以下内容

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);

或者简单地

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

- FlavioG

Answer 1

在internal::generic_product_impl的专业领域中，您需要将GemvProduct更改为GemmProduct。

- ggael