ConjugateGradient.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_CONJUGATE_GRADIENT_H
26 #define EIGEN_CONJUGATE_GRADIENT_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
41 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
43 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
44  const Preconditioner& precond, int& iters,
45  typename Dest::RealScalar& tol_error)
46 {
47  using std::sqrt;
48  using std::abs;
49  typedef typename Dest::RealScalar RealScalar;
50  typedef typename Dest::Scalar Scalar;
51  typedef Matrix<Scalar,Dynamic,1> VectorType;
52 
53  RealScalar tol = tol_error;
54  int maxIters = iters;
55 
56  int n = mat.cols();
57 
58  VectorType residual = rhs - mat * x; //initial residual
59  VectorType p(n);
60 
61  p = precond.solve(residual); //initial search direction
62 
63  VectorType z(n), tmp(n);
64  RealScalar absNew = internal::real(residual.dot(p)); // the square of the absolute value of r scaled by invM
65  RealScalar rhsNorm2 = rhs.squaredNorm();
66  RealScalar residualNorm2 = 0;
67  RealScalar threshold = tol*tol*rhsNorm2;
68  int i = 0;
69  while(i < maxIters)
70  {
71  tmp.noalias() = mat * p; // the bottleneck of the algorithm
72 
73  Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir
74  x += alpha * p; // update solution
75  residual -= alpha * tmp; // update residue
76 
77  residualNorm2 = residual.squaredNorm();
78  if(residualNorm2 < threshold)
79  break;
80 
81  z = precond.solve(residual); // approximately solve for "A z = residual"
82 
83  RealScalar absOld = absNew;
84  absNew = internal::real(residual.dot(z)); // update the absolute value of r
85  RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction
86  p = z + beta * p; // update search direction
87  i++;
88  }
89  tol_error = sqrt(residualNorm2 / rhsNorm2);
90  iters = i;
91 }
92 
93 }
94 
95 template< typename _MatrixType, int _UpLo=Lower,
96  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
97 class ConjugateGradient;
98 
99 namespace internal {
100 
101 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
102 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
103 {
104  typedef _MatrixType MatrixType;
105  typedef _Preconditioner Preconditioner;
106 };
107 
108 }
109 
158 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
159 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
160 {
162  using Base::mp_matrix;
163  using Base::m_error;
164  using Base::m_iterations;
165  using Base::m_info;
166  using Base::m_isInitialized;
167 public:
168  typedef _MatrixType MatrixType;
169  typedef typename MatrixType::Scalar Scalar;
170  typedef typename MatrixType::Index Index;
171  typedef typename MatrixType::RealScalar RealScalar;
172  typedef _Preconditioner Preconditioner;
173 
174  enum {
175  UpLo = _UpLo
176  };
177 
178 public:
179 
182 
194 
196 
202  template<typename Rhs,typename Guess>
203  inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
204  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
205  {
206  eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
207  eigen_assert(Base::rows()==b.rows()
208  && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
209  return internal::solve_retval_with_guess
210  <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
211  }
212 
214  template<typename Rhs,typename Dest>
215  void _solveWithGuess(const Rhs& b, Dest& x) const
216  {
219 
220  for(int j=0; j<b.cols(); ++j)
221  {
224 
225  typename Dest::ColXpr xj(x,j);
226  internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
228  }
229 
230  m_isInitialized = true;
232  }
233 
235  template<typename Rhs,typename Dest>
236  void _solve(const Rhs& b, Dest& x) const
237  {
238  x.setOnes();
239  _solveWithGuess(b,x);
240  }
241 
242 protected:
243 
244 };
245 
246 
247 namespace internal {
248 
249 template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
250 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
251  : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
252 {
253  typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
254  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
255 
256  template<typename Dest> void evalTo(Dest& dst) const
257  {
258  dec()._solve(rhs(),dst);
259  }
260 };
261 
262 } // end namespace internal
263 
264 } // end namespace Eigen
265 
266 #endif // EIGEN_CONJUGATE_GRADIENT_H