ConjugateGradient.h
Go to the documentation of this file.
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 namespace Eigen {
14 
15 namespace internal {
16 
26 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
28 void 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;
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 
95 template< typename _MatrixType, int _UpLo=Lower,
96  typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
98 
99 namespace internal {
100 
101 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
102 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
103 {
104  typedef _MatrixType MatrixType;
105  typedef _Preconditioner Preconditioner;
106 };
107 
108 }
109 
157 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
158 class 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;
166 public:
167  typedef _MatrixType MatrixType;
168  typedef typename MatrixType::Scalar Scalar;
170  typedef _Preconditioner Preconditioner;
171 
172  enum {
173  UpLo = _UpLo
174  };
175 
176 public:
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)
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 
223 protected:
224 
225 };
226 
227 } // end namespace Eigen
228 
229 #endif // EIGEN_CONJUGATE_GRADIENT_H
Eigen::internal::traits< ConjugateGradient< _MatrixType, _UpLo, _Preconditioner > >::MatrixType
_MatrixType MatrixType
Definition: ConjugateGradient.h:104
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::internal::traits< ConjugateGradient< _MatrixType, _UpLo, _Preconditioner > >::Preconditioner
_Preconditioner Preconditioner
Definition: ConjugateGradient.h:105
Eigen::ConjugateGradient::ConjugateGradient
ConjugateGradient()
Definition: ConjugateGradient.h:179
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
Eigen::ConjugateGradient::RealScalar
MatrixType::RealScalar RealScalar
Definition: ConjugateGradient.h:169
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::EigenBase
Definition: EigenBase.h:29
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Eigen::ConjugateGradient::Preconditioner
_Preconditioner Preconditioner
Definition: ConjugateGradient.h:170
Eigen::Upper
@ Upper
Definition: Constants.h:211
Eigen::Success
@ Success
Definition: Constants.h:442
real
float real
Definition: datatypes.h:10
type
Definition: pytypes.h:1525
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
beta
double beta(double a, double b)
Definition: beta.c:61
Eigen::ConjugateGradient::UpLo
@ UpLo
Definition: ConjugateGradient.h:173
n
int n
Definition: BiCGSTAB_simple.cpp:1
Eigen::NoConvergence
@ NoConvergence
Definition: Constants.h:446
A
Definition: test_numpy_dtypes.cpp:298
Eigen::Architecture::Type
Type
Definition: Constants.h:471
Eigen::IterativeSolverBase::maxIterations
Index maxIterations() const
Definition: IterativeSolverBase.h:281
Eigen::ConjugateGradient
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:97
Eigen::IterativeSolverBase::m_info
ComputationInfo m_info
Definition: IterativeSolverBase.h:438
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
Eigen::ConjugateGradient::Scalar
MatrixType::Scalar Scalar
Definition: ConjugateGradient.h:168
Eigen::Lower
@ Lower
Definition: Constants.h:209
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::ConjugateGradient::MatrixType
_MatrixType MatrixType
Definition: ConjugateGradient.h:167
Eigen::IterativeSolverBase
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:143
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::IterativeSolverBase::m_tolerance
RealScalar m_tolerance
Definition: IterativeSolverBase.h:434
Eigen::internal::traits
Definition: ForwardDeclarations.h:17
EIGEN_STATIC_ASSERT
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
Eigen::internal::Rhs
@ Rhs
Definition: TensorContractionMapper.h:18
Eigen::ConjugateGradient::~ConjugateGradient
~ConjugateGradient()
Definition: ConjugateGradient.h:194
Eigen::internal::conditional
Definition: Meta.h:109
Eigen::IterativeSolverBase::m_iterations
Index m_iterations
Definition: IterativeSolverBase.h:437
p
float * p
Definition: Tutorial_Map_using.cpp:9
Eigen::IterativeSolverBase::matrix
const ActualMatrixType & matrix() const
Definition: IterativeSolverBase.h:419
min
#define min(a, b)
Definition: datatypes.h:19
gtsam::tol
const G double tol
Definition: Group.h:79
Eigen::Matrix< Scalar, Dynamic, 1 >
abs
#define abs(x)
Definition: datatypes.h:17
Eigen::ConjugateGradient::Base
IterativeSolverBase< ConjugateGradient > Base
Definition: ConjugateGradient.h:160
internal
Definition: BandTriangularSolver.h:13
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
VectorType
Definition: FFTW.cpp:65
EIGEN_IMPLIES
#define EIGEN_IMPLIES(a, b)
Definition: Macros.h:1315
Eigen::MatrixWrapper
Expression of an array as a mathematical vector or matrix.
Definition: ArrayBase.h:15
Eigen::internal::generic_matrix_wrapper< MatrixType >
Eigen::IterativeSolverBase::ActualMatrixType
MatrixWrapper::ActualMatrixType ActualMatrixType
Definition: IterativeSolverBase.h:417
Eigen::IterativeSolverBase::m_preconditioner
Preconditioner m_preconditioner
Definition: IterativeSolverBase.h:431
Eigen::internal::conjugate_gradient
EIGEN_DONT_INLINE void conjugate_gradient(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters, typename Dest::RealScalar &tol_error)
Definition: ConjugateGradient.h:28
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::ConjugateGradient::ConjugateGradient
ConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: ConjugateGradient.h:192
Eigen::ConjugateGradient::_solve_vector_with_guess_impl
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: ConjugateGradient.h:198
Eigen::IterativeSolverBase::m_error
RealScalar m_error
Definition: IterativeSolverBase.h:436
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DONT_INLINE
#define EIGEN_DONT_INLINE
Definition: Macros.h:940
Eigen::SparseSolverBase::m_isInitialized
bool m_isInitialized
Definition: SparseSolverBase.h:119


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:10