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// 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
13namespace Eigen {
14
15namespace internal {
16
17/** \internal Low-level conjugate gradient algorithm
18  * \param mat The matrix A
19  * \param rhs The right hand side vector b
20  * \param x On input and initial solution, on output the computed solution.
21  * \param precond A preconditioner being able to efficiently solve for an
22  *                approximation of Ax=b (regardless of b)
23  * \param iters On input the max number of iteration, on output the number of performed iterations.
24  * \param tol_error On input the tolerance error, on output an estimation of the relative error.
25  */
26template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
27EIGEN_DONT_INLINE
28void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
29                        const Preconditioner& precond, int& 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  int maxIters = iters;
40
41  int 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  int 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 residue
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
94template< typename _MatrixType, int _UpLo=Lower,
95          typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
96class ConjugateGradient;
97
98namespace internal {
99
100template< typename _MatrixType, int _UpLo, typename _Preconditioner>
101struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
102{
103  typedef _MatrixType MatrixType;
104  typedef _Preconditioner Preconditioner;
105};
106
107}
108
109/** \ingroup IterativeLinearSolvers_Module
110  * \brief A conjugate gradient solver for sparse self-adjoint problems
111  *
112  * This class allows to solve for A.x = b sparse linear problems using a conjugate gradient algorithm.
113  * The sparse matrix A must be selfadjoint. The vectors x and b can be either dense or sparse.
114  *
115  * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
116  * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
117  *               or Upper. Default is Lower.
118  * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
119  *
120  * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
121  * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
122  * and NumTraits<Scalar>::epsilon() for the tolerance.
123  *
124  * This class can be used as the direct solver classes. Here is a typical usage example:
125  * \code
126  * int n = 10000;
127  * VectorXd x(n), b(n);
128  * SparseMatrix<double> A(n,n);
129  * // fill A and b
130  * ConjugateGradient<SparseMatrix<double> > cg;
131  * cg.compute(A);
132  * x = cg.solve(b);
133  * std::cout << "#iterations:     " << cg.iterations() << std::endl;
134  * std::cout << "estimated error: " << cg.error()      << std::endl;
135  * // update b, and solve again
136  * x = cg.solve(b);
137  * \endcode
138  *
139  * By default the iterations start with x=0 as an initial guess of the solution.
140  * One can control the start using the solveWithGuess() method. Here is a step by
141  * step execution example starting with a random guess and printing the evolution
142  * of the estimated error:
143  * * \code
144  * x = VectorXd::Random(n);
145  * cg.setMaxIterations(1);
146  * int i = 0;
147  * do {
148  *   x = cg.solveWithGuess(b,x);
149  *   std::cout << i << " : " << cg.error() << std::endl;
150  *   ++i;
151  * } while (cg.info()!=Success && i<100);
152  * \endcode
153  * Note that such a step by step excution is slightly slower.
154  *
155  * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
156  */
157template< typename _MatrixType, int _UpLo, typename _Preconditioner>
158class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
159{
160  typedef IterativeSolverBase<ConjugateGradient> Base;
161  using Base::mp_matrix;
162  using Base::m_error;
163  using Base::m_iterations;
164  using Base::m_info;
165  using Base::m_isInitialized;
166public:
167  typedef _MatrixType MatrixType;
168  typedef typename MatrixType::Scalar Scalar;
169  typedef typename MatrixType::Index Index;
170  typedef typename MatrixType::RealScalar RealScalar;
171  typedef _Preconditioner Preconditioner;
172
173  enum {
174    UpLo = _UpLo
175  };
176
177public:
178
179  /** Default constructor. */
180  ConjugateGradient() : Base() {}
181
182  /** Initialize the solver with matrix \a A for further \c Ax=b solving.
183    *
184    * This constructor is a shortcut for the default constructor followed
185    * by a call to compute().
186    *
187    * \warning this class stores a reference to the matrix A as well as some
188    * precomputed values that depend on it. Therefore, if \a A is changed
189    * this class becomes invalid. Call compute() to update it with the new
190    * matrix A, or modify a copy of A.
191    */
192  ConjugateGradient(const MatrixType& A) : Base(A) {}
193
194  ~ConjugateGradient() {}
195
196  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
197    * \a x0 as an initial solution.
198    *
199    * \sa compute()
200    */
201  template<typename Rhs,typename Guess>
202  inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
203  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
204  {
205    eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
206    eigen_assert(Base::rows()==b.rows()
207              && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
208    return internal::solve_retval_with_guess
209            <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
210  }
211
212  /** \internal */
213  template<typename Rhs,typename Dest>
214  void _solveWithGuess(const Rhs& b, Dest& x) const
215  {
216    m_iterations = Base::maxIterations();
217    m_error = Base::m_tolerance;
218
219    for(int j=0; j<b.cols(); ++j)
220    {
221      m_iterations = Base::maxIterations();
222      m_error = Base::m_tolerance;
223
224      typename Dest::ColXpr xj(x,j);
225      internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
226                                   Base::m_preconditioner, m_iterations, m_error);
227    }
228
229    m_isInitialized = true;
230    m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
231  }
232
233  /** \internal */
234  template<typename Rhs,typename Dest>
235  void _solve(const Rhs& b, Dest& x) const
236  {
237    x.setOnes();
238    _solveWithGuess(b,x);
239  }
240
241protected:
242
243};
244
245
246namespace internal {
247
248template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
249struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
250  : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
251{
252  typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
253  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
254
255  template<typename Dest> void evalTo(Dest& dst) const
256  {
257    dec()._solve(rhs(),dst);
258  }
259};
260
261} // end namespace internal
262
263} // end namespace Eigen
264
265#endif // EIGEN_CONJUGATE_GRADIENT_H
266