Eigen教程6 - Matrix-free solvers

博客新址: http://blog.xuezhisd.top
邮箱：xuezhisd@126.com

Matrix-free solvers

像ConjugateGradient 和 BiCGSTAB这样的迭代求解器可以用在 matrix free context。为此，用户必须提供一个继承EigenBase<>的封装类，并实现下面的方法：

• Index rows() 和 Index cols()：分别返回行数和列数。
• operator* ：操作数分别是你自己的类型和一个Eigen密集列向量。

Eigen::internal::traits<> 也必须为封装类型专门定制。

下面的例程，对Eigen::SparseMatrix进行了封装：

// 参考链接：http://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
#include <unsupported/Eigen/IterativeSolvers>
class MatrixReplacement; //声明类，定义在下面
using Eigen::SparseMatrix;
namespace Eigen {
	namespace internal {
		// MatrixReplacement和SparseMatrix相似，因此继承它的特性
		template<>
		struct traits<MatrixReplacement> :  public Eigen::internal::traits<Eigen::SparseMatrix<double> >
		{};
	}
}
// matrix-free wrapper的简单例子：将用户类型转换成Eigen兼容的类型
// 为了简单，该例子简单封装了Eigen::SparseMatrix
class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
public:
	// Required typedefs, constants, and method:
	typedef double Scalar;
	typedef double RealScalar;
	typedef int StorageIndex;
	enum {
		ColsAtCompileTime = Eigen::Dynamic,
		MaxColsAtCompileTime = Eigen::Dynamic,
		IsRowMajor = false
	};
	Index rows() const { return mp_mat->rows(); }
	Index cols() const { return mp_mat->cols(); }
	template<typename Rhs>
	Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const {
		return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived());
	}
	// 自定义API:
	MatrixReplacement() : mp_mat(0) {}
	void attachMyMatrix(const SparseMatrix<double> &mat) {
		mp_mat = &mat;
	}
	const SparseMatrix<double> my_matrix() const { return *mp_mat; }
private:
	const SparseMatrix<double> *mp_mat;
};
// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
namespace Eigen {
	namespace internal {
		template<typename Rhs>
		struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector
			: generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> >
		{
			typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar;
			template<typename Dest>
			static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha)
			{
				// This method should implement "dst += alpha * lhs * rhs" inplace,
				// however, for iterative solvers, alpha is always equal to 1, so let's not bother about it.
				assert(alpha==Scalar(1) && "scaling is not implemented");
				// Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
				// but let's do something fancier (and less efficient):
				for(Index i=0; i<lhs.cols(); ++i)
					dst += rhs(i) * lhs.my_matrix().col(i);
			}
		};
	}
}
int main()
{
	int n = 10;
	Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
	S = S.transpose()*S; //S'S 对称正定矩阵
	MatrixReplacement A;
	A.attachMyMatrix(S);
	Eigen::VectorXd b(n), x;
	b.setRandom();
	// 使用不同的迭代求解器，来求解 Ax = b with matrix-free version:
	{
		Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg;
		cg.compute(A);
		x = cg.solve(b);
		std::cout << "CG:       #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
	}
	{
		Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg;
		bicg.compute(A);
		x = bicg.solve(b);
		std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
	}
	{
		Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
		gmres.compute(A);
		x = gmres.solve(b);
		std::cout << "GMRES:    #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
	}
	{
		Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres;
		gmres.compute(A);
		x = gmres.solve(b);
		std::cout << "DGMRES:   #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
	}
	{
		Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres;
		minres.compute(A);
		x = minres.solve(b);
		std::cout << "MINRES:   #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
	}
}

参考

• http://eigen.tuxfamily.org/dox/group__MatrixfreeSolverExample.html