10 #ifndef EIGEN_CONJUGATE_GRADIENT_H
11 #define EIGEN_CONJUGATE_GRADIENT_H
13 #include "./InternalHeaderCheck.h"
28 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
30 void conjugate_gradient(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
31 const Preconditioner& precond,
Index& iters,
32 typename Dest::RealScalar& tol_error)
36 typedef typename Dest::RealScalar RealScalar;
37 typedef typename Dest::Scalar Scalar;
38 typedef Matrix<Scalar,Dynamic,1> VectorType;
40 RealScalar tol = tol_error;
41 Index maxIters = iters;
45 VectorType residual = rhs - mat * x;
47 RealScalar rhsNorm2 = rhs.squaredNorm();
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)
61 tol_error =
sqrt(residualNorm2 / rhsNorm2);
66 p = precond.solve(residual);
68 VectorType z(n), tmp(n);
69 RealScalar absNew = numext::real(residual.dot(p));
73 tmp.noalias() = mat * p;
75 Scalar alpha = absNew / p.dot(tmp);
77 residual -= alpha * tmp;
79 residualNorm2 = residual.squaredNorm();
80 if(residualNorm2 < threshold)
83 z = precond.solve(residual);
85 RealScalar absOld = absNew;
86 absNew = numext::real(residual.dot(z));
87 RealScalar beta = absNew / absOld;
91 tol_error =
sqrt(residualNorm2 / rhsNorm2);
97 template<
typename MatrixType_,
int UpLo_=
Lower,
98 typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
99 class ConjugateGradient;
103 template<
typename MatrixType_,
int UpLo_,
typename Preconditioner_>
104 struct traits<ConjugateGradient<MatrixType_,UpLo_,Preconditioner_> >
106 typedef MatrixType_ MatrixType;
107 typedef Preconditioner_ Preconditioner;
159 template<
typename MatrixType_,
int UpLo_,
typename Preconditioner_>
165 using Base::m_iterations;
167 using Base::m_isInitialized;
169 typedef MatrixType_ MatrixType;
170 typedef typename MatrixType::Scalar Scalar;
171 typedef typename MatrixType::RealScalar RealScalar;
172 typedef Preconditioner_ Preconditioner;
193 template<
typename MatrixDerived>
199 template<
typename Rhs,
typename Dest>
200 void _solve_vector_with_guess_impl(
const Rhs& b, Dest& x)
const
202 typedef typename Base::MatrixWrapper MatrixWrapper;
203 typedef typename Base::ActualMatrixType ActualMatrixType;
205 TransposeInput = (!MatrixWrapper::MatrixFree)
207 && (!MatrixType::IsRowMajor)
208 && (!NumTraits<Scalar>::IsComplex)
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),
214 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
215 > SelfAdjointWrapper;
218 m_error = Base::m_tolerance;
220 RowMajorWrapper row_mat(matrix());
221 internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
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