Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
ConjugateGradient.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011-2014 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 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 namespace internal {
18 
28 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
29 EIGEN_DONT_INLINE
30 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
31  const Preconditioner& precond, Index& iters,
32  typename Dest::RealScalar& tol_error)
33 {
34  using std::sqrt;
35  using std::abs;
36  typedef typename Dest::RealScalar RealScalar;
37  typedef typename Dest::Scalar Scalar;
38  typedef Matrix<Scalar,Dynamic,1> VectorType;
39 
40  RealScalar tol = tol_error;
41  Index maxIters = iters;
42 
43  Index n = mat.cols();
44 
45  VectorType residual = rhs - mat * x; //initial residual
46 
47  RealScalar rhsNorm2 = rhs.squaredNorm();
48  if(rhsNorm2 == 0)
49  {
50  x.setZero();
51  iters = 0;
52  tol_error = 0;
53  return;
54  }
55  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
56  RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero);
57  RealScalar residualNorm2 = residual.squaredNorm();
58  if (residualNorm2 < threshold)
59  {
60  iters = 0;
61  tol_error = sqrt(residualNorm2 / rhsNorm2);
62  return;
63  }
64 
65  VectorType p(n);
66  p = precond.solve(residual); // initial search direction
67 
68  VectorType z(n), tmp(n);
69  RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM
70  Index i = 0;
71  while(i < maxIters)
72  {
73  tmp.noalias() = mat * p; // the bottleneck of the algorithm
74 
75  Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir
76  x += alpha * p; // update solution
77  residual -= alpha * tmp; // update residual
78 
79  residualNorm2 = residual.squaredNorm();
80  if(residualNorm2 < threshold)
81  break;
82 
83  z = precond.solve(residual); // approximately solve for "A z = residual"
84 
85  RealScalar absOld = absNew;
86  absNew = numext::real(residual.dot(z)); // update the absolute value of r
87  RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction
88  p = z + beta * p; // update search direction
89  i++;
90  }
91  tol_error = sqrt(residualNorm2 / rhsNorm2);
92  iters = i;
93 }
94 
95 }
96 
97 template< typename MatrixType_, int UpLo_=Lower,
98  typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
99 class ConjugateGradient;
100 
101 namespace internal {
102 
103 template< typename MatrixType_, int UpLo_, typename Preconditioner_>
104 struct traits<ConjugateGradient<MatrixType_,UpLo_,Preconditioner_> >
105 {
106  typedef MatrixType_ MatrixType;
107  typedef Preconditioner_ Preconditioner;
108 };
109 
110 }
111 
159 template< typename MatrixType_, int UpLo_, typename Preconditioner_>
160 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<MatrixType_,UpLo_,Preconditioner_> >
161 {
163  using Base::matrix;
164  using Base::m_error;
165  using Base::m_iterations;
166  using Base::m_info;
167  using Base::m_isInitialized;
168 public:
169  typedef MatrixType_ MatrixType;
170  typedef typename MatrixType::Scalar Scalar;
171  typedef typename MatrixType::RealScalar RealScalar;
172  typedef Preconditioner_ Preconditioner;
173 
174  enum {
175  UpLo = UpLo_
176  };
177 
178 public:
179 
182 
193  template<typename MatrixDerived>
194  explicit ConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
195 
196  ~ConjugateGradient() {}
197 
199  template<typename Rhs,typename Dest>
200  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
201  {
202  typedef typename Base::MatrixWrapper MatrixWrapper;
203  typedef typename Base::ActualMatrixType ActualMatrixType;
204  enum {
205  TransposeInput = (!MatrixWrapper::MatrixFree)
206  && (UpLo==(Lower|Upper))
207  && (!MatrixType::IsRowMajor)
208  && (!NumTraits<Scalar>::IsComplex)
209  };
210  typedef std::conditional_t<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&> RowMajorWrapper;
211  EIGEN_STATIC_ASSERT(internal::check_implication(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
212  typedef std::conditional_t<UpLo==(Lower|Upper),
213  RowMajorWrapper,
214  typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
215  > SelfAdjointWrapper;
216 
217  m_iterations = Base::maxIterations();
218  m_error = Base::m_tolerance;
219 
220  RowMajorWrapper row_mat(matrix());
221  internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
222  m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
223  }
224 
225 protected:
226 
227 };
228 
229 } // end namespace Eigen
230 
231 #endif // EIGEN_CONJUGATE_GRADIENT_H
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:161
ConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: ConjugateGradient.h:194
ConjugateGradient()
Definition: ConjugateGradient.h:181
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:146
Index maxIterations() const
Definition: IterativeSolverBase.h:286
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ Success
Definition: Constants.h:444
@ NoConvergence
Definition: Constants.h:448
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32