Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace Eigen {
16
17namespace internal {
18
28template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
29EIGEN_DONT_INLINE
30void 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
97template< typename MatrixType_, int UpLo_=Lower,
98 typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
99class ConjugateGradient;
100
101namespace internal {
102
103template< typename MatrixType_, int UpLo_, typename Preconditioner_>
104struct traits<ConjugateGradient<MatrixType_,UpLo_,Preconditioner_> >
105{
106 typedef MatrixType_ MatrixType;
107 typedef Preconditioner_ Preconditioner;
108};
109
110}
111
159template< typename MatrixType_, int UpLo_, typename Preconditioner_>
160class 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;
168public:
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
178public:
179
182
193 template<typename MatrixDerived>
194 explicit ConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
195
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 typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type 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 typename internal::conditional<UpLo==(Lower|Upper),
213 RowMajorWrapper,
214 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
215 >::type 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
225protected:
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:283
@ 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: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: EigenBase.h:32