matrixfree_cg.cpp
Go to the documentation of this file.
1 #include <iostream>
2 #include <Eigen/Core>
3 #include <Eigen/Dense>
4 #include <Eigen/IterativeLinearSolvers>
5 #include <unsupported/Eigen/IterativeSolvers>
6 
9 
10 namespace Eigen {
11 namespace 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.
21 class MatrixReplacement : public Eigen::EigenBase<MatrixReplacement> {
22 public:
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>
39  }
40 
41  // Custom API:
42  MatrixReplacement() : mp_mat(0) {}
43 
45  mp_mat = &mat;
46  }
47  const SparseMatrix<double> my_matrix() const { return *mp_mat; }
48 
49 private:
51 };
52 
53 
54 // Implementation of MatrixReplacement * Eigen::DenseVector though a specialization of internal::generic_product_impl:
55 namespace Eigen {
56 namespace 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  {
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");
71 
72  // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs,
73  // but let's do something fancier (and less efficient):
74  for(Index i=0; i<lhs.cols(); ++i)
75  dst += rhs(i) * lhs.my_matrix().col(i);
76  }
77  };
78 
79 }
80 }
81 
82 int main()
83 {
84  int n = 10;
85  Eigen::SparseMatrix<double> S = Eigen::MatrixXd::Random(n,n).sparseView(0.5,1);
86  S = S.transpose()*S;
87 
89  A.attachMyMatrix(S);
90 
91  Eigen::VectorXd b(n), x;
92  b.setRandom();
93 
94  // Solve Ax = b using various iterative solver with matrix-free version:
95  {
97  cg.compute(A);
98  x = cg.solve(b);
99  std::cout << "CG: #iterations: " << cg.iterations() << ", estimated error: " << cg.error() << std::endl;
100  }
101 
102  {
104  bicg.compute(A);
105  x = bicg.solve(b);
106  std::cout << "BiCGSTAB: #iterations: " << bicg.iterations() << ", estimated error: " << bicg.error() << std::endl;
107  }
108 
109  {
111  gmres.compute(A);
112  x = gmres.solve(b);
113  std::cout << "GMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
114  }
115 
116  {
118  gmres.compute(A);
119  x = gmres.solve(b);
120  std::cout << "DGMRES: #iterations: " << gmres.iterations() << ", estimated error: " << gmres.error() << std::endl;
121  }
122 
123  {
125  minres.compute(A);
126  x = minres.solve(b);
127  std::cout << "MINRES: #iterations: " << minres.iterations() << ", estimated error: " << minres.error() << std::endl;
128  }
129 }
SCALAR Scalar
Definition: bench_gemm.cpp:46
Expression of the product of two arbitrary matrices or vectors.
Definition: Product.h:71
Scalar * b
Definition: benchVecAdd.cpp:17
A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse...
Definition: DGMRES.h:19
A versatible sparse matrix representation.
Definition: SparseMatrix.h:96
int n
const Solve< ConjugateGradient< _MatrixType, _UpLo, _Preconditioner >, Rhs > solve(const MatrixBase< Rhs > &b) const
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
void attachMyMatrix(const SparseMatrix< double > &mat)
DiscreteKey S(1, 2)
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ColXpr col(Index i)
ConjugateGradient< _MatrixType, _UpLo, _Preconditioner > & compute(const EigenBase< MatrixDerived > &A)
TransposeReturnType transpose()
A minimal residual solver for sparse symmetric problems.
Definition: MINRES.h:143
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Index cols() const
RealScalar alpha
const SparseMatrix< double > my_matrix() const
Eigen::Product< MatrixReplacement, Rhs, Eigen::AliasFreeProduct > operator*(const Eigen::MatrixBase< Rhs > &x) const
A GMRES solver for sparse square problems.
Definition: GMRES.h:221
Index rows() const
int main()
A bi conjugate gradient stabilized solver for sparse square problems.
Definition: BiCGSTAB.h:113
const int Dynamic
Definition: Constants.h:22
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
static void scaleAndAddTo(Dest &dst, const MatrixReplacement &lhs, const Rhs &rhs, const Scalar &alpha)
const SparseMatrix< double > * mp_mat
#define EIGEN_ONLY_USED_FOR_DEBUG(x)
Definition: Macros.h:1049
EIGEN_DONT_INLINE void minres(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters, typename Dest::RealScalar &tol_error)
Definition: MINRES.h:32
bool gmres(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters, const Index &restart, typename Dest::RealScalar &tol_error)
Definition: GMRES.h:56


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:52