• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
13 namespace Eigen {
14 
15 namespace 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   */
26 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
27 EIGEN_DONT_INLINE
conjugate_gradient(const MatrixType & mat,const Rhs & rhs,Dest & x,const Preconditioner & precond,int & iters,typename Dest::RealScalar & tol_error)28 void 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   VectorType p(n);
45 
46   p = precond.solve(residual);      //initial search direction
47 
48   VectorType z(n), tmp(n);
49   RealScalar absNew = internal::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
50   RealScalar rhsNorm2 = rhs.squaredNorm();
51   RealScalar residualNorm2 = 0;
52   RealScalar threshold = tol*tol*rhsNorm2;
53   int i = 0;
54   while(i < maxIters)
55   {
56     tmp.noalias() = mat * p;              // the bottleneck of the algorithm
57 
58     Scalar alpha = absNew / p.dot(tmp);   // the amount we travel on dir
59     x += alpha * p;                       // update solution
60     residual -= alpha * tmp;              // update residue
61 
62     residualNorm2 = residual.squaredNorm();
63     if(residualNorm2 < threshold)
64       break;
65 
66     z = precond.solve(residual);          // approximately solve for "A z = residual"
67 
68     RealScalar absOld = absNew;
69     absNew = internal::real(residual.dot(z));     // update the absolute value of r
70     RealScalar beta = absNew / absOld;            // calculate the Gram-Schmidt value used to create the new search direction
71     p = z + beta * p;                             // update search direction
72     i++;
73   }
74   tol_error = sqrt(residualNorm2 / rhsNorm2);
75   iters = i;
76 }
77 
78 }
79 
80 template< typename _MatrixType, int _UpLo=Lower,
81           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
82 class ConjugateGradient;
83 
84 namespace internal {
85 
86 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
87 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
88 {
89   typedef _MatrixType MatrixType;
90   typedef _Preconditioner Preconditioner;
91 };
92 
93 }
94 
95 /** \ingroup IterativeLinearSolvers_Module
96   * \brief A conjugate gradient solver for sparse self-adjoint problems
97   *
98   * This class allows to solve for A.x = b sparse linear problems using a conjugate gradient algorithm.
99   * The sparse matrix A must be selfadjoint. The vectors x and b can be either dense or sparse.
100   *
101   * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
102   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
103   *               or Upper. Default is Lower.
104   * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
105   *
106   * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
107   * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
108   * and NumTraits<Scalar>::epsilon() for the tolerance.
109   *
110   * This class can be used as the direct solver classes. Here is a typical usage example:
111   * \code
112   * int n = 10000;
113   * VectorXd x(n), b(n);
114   * SparseMatrix<double> A(n,n);
115   * // fill A and b
116   * ConjugateGradient<SparseMatrix<double> > cg;
117   * cg.compute(A);
118   * x = cg.solve(b);
119   * std::cout << "#iterations:     " << cg.iterations() << std::endl;
120   * std::cout << "estimated error: " << cg.error()      << std::endl;
121   * // update b, and solve again
122   * x = cg.solve(b);
123   * \endcode
124   *
125   * By default the iterations start with x=0 as an initial guess of the solution.
126   * One can control the start using the solveWithGuess() method. Here is a step by
127   * step execution example starting with a random guess and printing the evolution
128   * of the estimated error:
129   * * \code
130   * x = VectorXd::Random(n);
131   * cg.setMaxIterations(1);
132   * int i = 0;
133   * do {
134   *   x = cg.solveWithGuess(b,x);
135   *   std::cout << i << " : " << cg.error() << std::endl;
136   *   ++i;
137   * } while (cg.info()!=Success && i<100);
138   * \endcode
139   * Note that such a step by step excution is slightly slower.
140   *
141   * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
142   */
143 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
144 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
145 {
146   typedef IterativeSolverBase<ConjugateGradient> Base;
147   using Base::mp_matrix;
148   using Base::m_error;
149   using Base::m_iterations;
150   using Base::m_info;
151   using Base::m_isInitialized;
152 public:
153   typedef _MatrixType MatrixType;
154   typedef typename MatrixType::Scalar Scalar;
155   typedef typename MatrixType::Index Index;
156   typedef typename MatrixType::RealScalar RealScalar;
157   typedef _Preconditioner Preconditioner;
158 
159   enum {
160     UpLo = _UpLo
161   };
162 
163 public:
164 
165   /** Default constructor. */
166   ConjugateGradient() : Base() {}
167 
168   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
169     *
170     * This constructor is a shortcut for the default constructor followed
171     * by a call to compute().
172     *
173     * \warning this class stores a reference to the matrix A as well as some
174     * precomputed values that depend on it. Therefore, if \a A is changed
175     * this class becomes invalid. Call compute() to update it with the new
176     * matrix A, or modify a copy of A.
177     */
178   ConjugateGradient(const MatrixType& A) : Base(A) {}
179 
180   ~ConjugateGradient() {}
181 
182   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
183     * \a x0 as an initial solution.
184     *
185     * \sa compute()
186     */
187   template<typename Rhs,typename Guess>
188   inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
189   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
190   {
191     eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
192     eigen_assert(Base::rows()==b.rows()
193               && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
194     return internal::solve_retval_with_guess
195             <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
196   }
197 
198   /** \internal */
199   template<typename Rhs,typename Dest>
200   void _solveWithGuess(const Rhs& b, Dest& x) const
201   {
202     m_iterations = Base::maxIterations();
203     m_error = Base::m_tolerance;
204 
205     for(int j=0; j<b.cols(); ++j)
206     {
207       m_iterations = Base::maxIterations();
208       m_error = Base::m_tolerance;
209 
210       typename Dest::ColXpr xj(x,j);
211       internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
212                                    Base::m_preconditioner, m_iterations, m_error);
213     }
214 
215     m_isInitialized = true;
216     m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
217   }
218 
219   /** \internal */
220   template<typename Rhs,typename Dest>
221   void _solve(const Rhs& b, Dest& x) const
222   {
223     x.setOnes();
224     _solveWithGuess(b,x);
225   }
226 
227 protected:
228 
229 };
230 
231 
232 namespace internal {
233 
234 template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
235 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
236   : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
237 {
238   typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
239   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
240 
241   template<typename Dest> void evalTo(Dest& dst) const
242   {
243     dec()._solve(rhs(),dst);
244   }
245 };
246 
247 } // end namespace internal
248 
249 } // end namespace Eigen
250 
251 #endif // EIGEN_CONJUGATE_GRADIENT_H
252