具有特征库的 MatrixFree-Matrix(和 Vector)乘积
MatrixFree-Matrix (and Vector) product with Eigen Library
我正在创建一个 C++ 软件,我需要一个包装器,它基于 Eigen 库,实现类似于官方网页中解释的运算符*
https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html
使用上述网页中的代码,我可以将模板 class 包装在 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
在internal::generic_product_impl
的专精中,需要将GemvProduct
改为GemmProduct
。
我正在创建一个 C++ 软件,我需要一个包装器,它基于 Eigen 库,实现类似于官方网页中解释的运算符*
https://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html
使用上述网页中的代码,我可以将模板 class 包装在 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
在internal::generic_product_impl
的专精中,需要将GemvProduct
改为GemmProduct
。