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 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  *               Upper, or Lower|Upper in which the full matrix entries will be considered. 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.
141  *
142  * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
143  */
144template< typename _MatrixType, int _UpLo, typename _Preconditioner>
145class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
146{
147  typedef IterativeSolverBase<ConjugateGradient> Base;
148  using Base::mp_matrix;
149  using Base::m_error;
150  using Base::m_iterations;
151  using Base::m_info;
152  using Base::m_isInitialized;
153public:
154  typedef _MatrixType MatrixType;
155  typedef typename MatrixType::Scalar Scalar;
156  typedef typename MatrixType::Index Index;
157  typedef typename MatrixType::RealScalar RealScalar;
158  typedef _Preconditioner Preconditioner;
159
160  enum {
161    UpLo = _UpLo
162  };
163
164public:
165
166  /** Default constructor. */
167  ConjugateGradient() : Base() {}
168
169  /** Initialize the solver with matrix \a A for further \c Ax=b solving.
170    *
171    * This constructor is a shortcut for the default constructor followed
172    * by a call to compute().
173    *
174    * \warning this class stores a reference to the matrix A as well as some
175    * precomputed values that depend on it. Therefore, if \a A is changed
176    * this class becomes invalid. Call compute() to update it with the new
177    * matrix A, or modify a copy of A.
178    */
179  ConjugateGradient(const MatrixType& A) : Base(A) {}
180
181  ~ConjugateGradient() {}
182
183  /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
184    * \a x0 as an initial solution.
185    *
186    * \sa compute()
187    */
188  template<typename Rhs,typename Guess>
189  inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
190  solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
191  {
192    eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
193    eigen_assert(Base::rows()==b.rows()
194              && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
195    return internal::solve_retval_with_guess
196            <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
197  }
198
199  /** \internal */
200  template<typename Rhs,typename Dest>
201  void _solveWithGuess(const Rhs& b, Dest& x) const
202  {
203    typedef typename internal::conditional<UpLo==(Lower|Upper),
204                                           const MatrixType&,
205                                           SparseSelfAdjointView<const MatrixType, UpLo>
206                                          >::type MatrixWrapperType;
207    m_iterations = Base::maxIterations();
208    m_error = Base::m_tolerance;
209
210    for(int j=0; j<b.cols(); ++j)
211    {
212      m_iterations = Base::maxIterations();
213      m_error = Base::m_tolerance;
214
215      typename Dest::ColXpr xj(x,j);
216      internal::conjugate_gradient(MatrixWrapperType(*mp_matrix), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
217    }
218
219    m_isInitialized = true;
220    m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
221  }
222
223  /** \internal */
224  template<typename Rhs,typename Dest>
225  void _solve(const Rhs& b, Dest& x) const
226  {
227    x.setZero();
228    _solveWithGuess(b,x);
229  }
230
231protected:
232
233};
234
235
236namespace internal {
237
238template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
239struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
240  : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
241{
242  typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
243  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
244
245  template<typename Dest> void evalTo(Dest& dst) const
246  {
247    dec()._solve(rhs(),dst);
248  }
249};
250
251} // end namespace internal
252
253} // end namespace Eigen
254
255#endif // EIGEN_CONJUGATE_GRADIENT_H
256