Eigen  3.2.91
ConjugateGradient.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_CONJUGATE_GRADIENT_H
11 #define EIGEN_CONJUGATE_GRADIENT_H
12 
13 namespace Eigen {
14 
15 namespace internal {
16 
26 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
27 EIGEN_DONT_INLINE
28 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
29  const Preconditioner& precond, Index& iters,
30  typename Dest::RealScalar& tol_error)
31 {
32  using std::sqrt;
33  using std::abs;
34  typedef typename Dest::RealScalar RealScalar;
35  typedef typename Dest::Scalar Scalar;
36  typedef Matrix<Scalar,Dynamic,1> VectorType;
37 
38  RealScalar tol = tol_error;
39  Index maxIters = iters;
40 
41  Index n = mat.cols();
42 
43  VectorType residual = rhs - mat * x; //initial residual
44 
45  RealScalar rhsNorm2 = rhs.squaredNorm();
46  if(rhsNorm2 == 0)
47  {
48  x.setZero();
49  iters = 0;
50  tol_error = 0;
51  return;
52  }
53  RealScalar threshold = tol*tol*rhsNorm2;
54  RealScalar residualNorm2 = residual.squaredNorm();
55  if (residualNorm2 < threshold)
56  {
57  iters = 0;
58  tol_error = sqrt(residualNorm2 / rhsNorm2);
59  return;
60  }
61 
62  VectorType p(n);
63  p = precond.solve(residual); // initial search direction
64 
65  VectorType z(n), tmp(n);
66  RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM
67  Index i = 0;
68  while(i < maxIters)
69  {
70  tmp.noalias() = mat * p; // the bottleneck of the algorithm
71 
72  Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir
73  x += alpha * p; // update solution
74  residual -= alpha * tmp; // update residual
75 
76  residualNorm2 = residual.squaredNorm();
77  if(residualNorm2 < threshold)
78  break;
79 
80  z = precond.solve(residual); // approximately solve for "A z = residual"
81 
82  RealScalar absOld = absNew;
83  absNew = numext::real(residual.dot(z)); // update the absolute value of r
84  RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction
85  p = z + beta * p; // update search direction
86  i++;
87  }
88  tol_error = sqrt(residualNorm2 / rhsNorm2);
89  iters = i;
90 }
91 
92 }
93 
94 template< typename _MatrixType, int _UpLo=Lower,
95  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
97 
98 namespace internal {
99 
100 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
101 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
102 {
103  typedef _MatrixType MatrixType;
104  typedef _Preconditioner Preconditioner;
105 };
106 
107 }
108 
152 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
153 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
154 {
155  typedef IterativeSolverBase<ConjugateGradient> Base;
156  using Base::mp_matrix;
157  using Base::m_error;
158  using Base::m_iterations;
159  using Base::m_info;
160  using Base::m_isInitialized;
161 public:
162  typedef _MatrixType MatrixType;
163  typedef typename MatrixType::Scalar Scalar;
164  typedef typename MatrixType::RealScalar RealScalar;
165  typedef _Preconditioner Preconditioner;
166 
167  enum {
168  UpLo = _UpLo
169  };
170 
171 public:
172 
174  ConjugateGradient() : Base() {}
175 
186  explicit ConjugateGradient(const MatrixType& A) : Base(A) {}
187 
188  ~ConjugateGradient() {}
189 
191  template<typename Rhs,typename Dest>
192  void _solve_with_guess_impl(const Rhs& b, Dest& x) const
193  {
194  typedef Ref<const MatrixType> MatRef;
195  typedef typename internal::conditional<UpLo==(Lower|Upper) && (!MatrixType::IsRowMajor) && (!NumTraits<Scalar>::IsComplex),
196  Transpose<const MatRef>, MatRef const&>::type RowMajorWrapper;
197  typedef typename internal::conditional<UpLo==(Lower|Upper),
198  RowMajorWrapper,
199  typename MatRef::template ConstSelfAdjointViewReturnType<UpLo>::Type
200  >::type SelfAdjointWrapper;
201  m_iterations = Base::maxIterations();
202  m_error = Base::m_tolerance;
203 
204  for(Index j=0; j<b.cols(); ++j)
205  {
206  m_iterations = Base::maxIterations();
207  m_error = Base::m_tolerance;
208 
209  typename Dest::ColXpr xj(x,j);
210  RowMajorWrapper row_mat(mp_matrix);
211  internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
212  }
213 
214  m_isInitialized = true;
215  m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
216  }
217 
219  using Base::_solve_impl;
220  template<typename Rhs,typename Dest>
221  void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const
222  {
223  x.setZero();
224  _solve_with_guess_impl(b.derived(),x);
225  }
226 
227 protected:
228 
229 };
230 
231 } // end namespace Eigen
232 
233 #endif // EIGEN_CONJUGATE_GRADIENT_H
ConjugateGradient(const MatrixType &A)
Definition: ConjugateGradient.h:186
Definition: Constants.h:196
Definition: LDLT.h:16
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:96
ConjugateGradient()
Definition: ConjugateGradient.h:174
Index maxIterations() const
Definition: IterativeSolverBase.h:155
Definition: Constants.h:198
Definition: Constants.h:424
Definition: Eigen_Colamd.h:54
Definition: Constants.h:428