Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
13namespace Eigen {
14
15namespace internal {
16
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, Index& 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 Index maxIters = iters;
40
41 Index 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 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
54 RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero);
55 RealScalar residualNorm2 = residual.squaredNorm();
56 if (residualNorm2 < threshold)
57 {
58 iters = 0;
59 tol_error = sqrt(residualNorm2 / rhsNorm2);
60 return;
61 }
62
63 VectorType p(n);
64 p = precond.solve(residual); // initial search direction
65
66 VectorType z(n), tmp(n);
67 RealScalar absNew = numext::real(residual.dot(p)); // the square of the absolute value of r scaled by invM
68 Index i = 0;
69 while(i < maxIters)
70 {
71 tmp.noalias() = mat * p; // the bottleneck of the algorithm
72
73 Scalar alpha = absNew / p.dot(tmp); // the amount we travel on dir
74 x += alpha * p; // update solution
75 residual -= alpha * tmp; // update residual
76
77 residualNorm2 = residual.squaredNorm();
78 if(residualNorm2 < threshold)
79 break;
80
81 z = precond.solve(residual); // approximately solve for "A z = residual"
82
83 RealScalar absOld = absNew;
84 absNew = numext::real(residual.dot(z)); // update the absolute value of r
85 RealScalar beta = absNew / absOld; // calculate the Gram-Schmidt value used to create the new search direction
86 p = z + beta * p; // update search direction
87 i++;
88 }
89 tol_error = sqrt(residualNorm2 / rhsNorm2);
90 iters = i;
91}
92
93}
94
95template< typename _MatrixType, int _UpLo=Lower,
96 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
97class ConjugateGradient;
98
99namespace internal {
100
101template< typename _MatrixType, int _UpLo, typename _Preconditioner>
102struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
103{
104 typedef _MatrixType MatrixType;
105 typedef _Preconditioner Preconditioner;
106};
107
108}
109
157template< typename _MatrixType, int _UpLo, typename _Preconditioner>
158class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
159{
161 using Base::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::RealScalar RealScalar;
170 typedef _Preconditioner Preconditioner;
171
172 enum {
173 UpLo = _UpLo
174 };
175
176public:
177
180
191 template<typename MatrixDerived>
192 explicit ConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
193
195
197 template<typename Rhs,typename Dest>
198 void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
199 {
200 typedef typename Base::MatrixWrapper MatrixWrapper;
201 typedef typename Base::ActualMatrixType ActualMatrixType;
202 enum {
203 TransposeInput = (!MatrixWrapper::MatrixFree)
204 && (UpLo==(Lower|Upper))
205 && (!MatrixType::IsRowMajor)
206 && (!NumTraits<Scalar>::IsComplex)
207 };
208 typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
209 EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
210 typedef typename internal::conditional<UpLo==(Lower|Upper),
211 RowMajorWrapper,
212 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
213 >::type SelfAdjointWrapper;
214
215 m_iterations = Base::maxIterations();
216 m_error = Base::m_tolerance;
217
218 RowMajorWrapper row_mat(matrix());
219 internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
220 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
221 }
222
223protected:
224
225};
226
227} // end namespace Eigen
228
229#endif // EIGEN_CONJUGATE_GRADIENT_H
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:159
ConjugateGradient()
Definition: ConjugateGradient.h:179
ConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: ConjugateGradient.h:192
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:144
Index maxIterations() const
Definition: IterativeSolverBase.h:281
@ Lower
Definition: Constants.h:209
@ Upper
Definition: Constants.h:211
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30