matrixfree_cg.cpp revision 2b8756b6f1de65d3f8bffab45be6c44ceb7411fc
1#include <iostream> 2#include <Eigen/Core> 3#include <Eigen/Dense> 4#include <Eigen/IterativeLinearSolvers> 5#include <unsupported/Eigen/IterativeSolvers> 6 7class MatrixReplacement; 8using Eigen::SparseMatrix; 9 10namespace Eigen { 11namespace internal { 12 // MatrixReplacement looks-like a SparseMatrix, so let's inherits its traits: 13 template<> 14 struct traits<MatrixReplacement> : public Eigen::internal::traits<Eigen::SparseMatrix<double> > 15 {}; 16} 17} 18 19// Example of a matrix-free wrapper from a user type to Eigen's compatible type 20// For the sake of simplicity, this example simply wrap a Eigen::SparseMatrix. 21class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> { 22public: 23 // Required typedefs, constants, and method: 24 typedef double Scalar; 25 typedef double RealScalar; 26 typedef int StorageIndex; 27 enum { 28 ColsAtCompileTime = Eigen::Dynamic, 29 MaxColsAtCompileTime = Eigen::Dynamic, 30 IsRowMajor = false 31 }; 32 33 Index rows() const { return mp_mat->rows(); } 34 Index cols() const { return mp_mat->cols(); } 35 36 template<typename Rhs> 37 Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct> operator*(const Eigen::MatrixBase<Rhs>& x) const { 38 return Eigen::Product<MatrixReplacement,Rhs,Eigen::AliasFreeProduct>(*this, x.derived()); 39 } 40 41 // Custom API: 42 MatrixReplacement() : mp_mat(0) {} 43 44 void attachMyMatrix(const SparseMatrix<double> &mat) { 45 mp_mat = &mat; 46 } 47 const SparseMatrix<double> my_matrix() const { return *mp_mat; } 48 49private: 50 const SparseMatrix<double> *mp_mat; 51}; 52 53 54// Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl: 55namespace Eigen { 56namespace internal { 57 58 template<typename Rhs> 59 struct generic_product_impl<MatrixReplacement, Rhs, SparseShape, DenseShape, GemvProduct> // GEMV stands for matrix-vector 60 : generic_product_impl_base<MatrixReplacement,Rhs,generic_product_impl<MatrixReplacement,Rhs> > 61 { 62 typedef typename Product<MatrixReplacement,Rhs>::Scalar Scalar; 63 64 template<typename Dest> 65 static void scaleAndAddTo(Dest& dst, const MatrixReplacement& lhs, const Rhs& rhs, const Scalar& alpha) 66 { 67 // This method should implement "dst += alpha * lhs * rhs" inplace, 68 // however, for iterative solvers, alpha is always equal to 1, so let's not bother about it. 69 assert(alpha==Scalar(1) && "scaling is not implemented"); 70 71 // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, 72 // but let's do something fancier (and less efficient): 73 for(Index i=0; i<lhs.cols(); ++i) 74 dst += rhs(i) * lhs.my_matrix().col(i); 75 } 76 }; 77 78} 79} 80 81int main() 82{ 83 int n = 10; 84 Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1); 85 S = S.transpose()*S; 86 87 MatrixReplacement A; 88 A.attachMyMatrix(S); 89 90 Eigen::VectorXd b(n), x; 91 b.setRandom(); 92 93 // Solve Ax = b using various iterative solver with matrix-free version: 94 { 95 Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> cg; 96 cg.compute(A); 97 x = cg.solve(b); 98 std::cout << "CG: #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl; 99 } 100 101 { 102 Eigen::BiCGSTAB<MatrixReplacement, Eigen::IdentityPreconditioner> bicg; 103 bicg.compute(A); 104 x = bicg.solve(b); 105 std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl; 106 } 107 108 { 109 Eigen::GMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres; 110 gmres.compute(A); 111 x = gmres.solve(b); 112 std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; 113 } 114 115 { 116 Eigen::DGMRES<MatrixReplacement, Eigen::IdentityPreconditioner> gmres; 117 gmres.compute(A); 118 x = gmres.solve(b); 119 std::cout << "DGMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl; 120 } 121 122 { 123 Eigen::MINRES<MatrixReplacement, Eigen::Lower|Eigen::Upper, Eigen::IdentityPreconditioner> minres; 124 minres.compute(A); 125 x = minres.solve(b); 126 std::cout << "MINRES: #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl; 127 } 128} 129