Free cookie consent management tool by TermsFeed Policy Generator

source: branches/HeuristicLab.Problems.GaussianProcessTuning/HeuristicLab.Eigen/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h @ 9562

Last change on this file since 9562 was 9562, checked in by gkronber, 11 years ago

#1967 worked on Gaussian process evolution.

File size: 8.5 KB
Line 
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  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
80template< typename _MatrixType, int _UpLo=Lower,
81          typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
82class ConjugateGradient;
83
84namespace internal {
85
86template< typename _MatrixType, int _UpLo, typename _Preconditioner>
87struct 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  */
143template< typename _MatrixType, int _UpLo, typename _Preconditioner>
144class 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;
152public:
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
163public:
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
227protected:
228
229};
230
231
232namespace internal {
233
234template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
235struct 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
Note: See TracBrowser for help on using the repository browser.