10 #ifndef EIGEN_CONJUGATE_GRADIENT_H 11 #define EIGEN_CONJUGATE_GRADIENT_H 26 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
29 const Preconditioner& precond,
int& iters,
30 typename Dest::RealScalar& tol_error)
34 typedef typename Dest::RealScalar RealScalar;
35 typedef typename Dest::Scalar Scalar;
38 RealScalar tol = tol_error;
43 VectorType residual = rhs - mat * x;
45 RealScalar rhsNorm2 = rhs.squaredNorm();
53 RealScalar threshold = tol*tol*rhsNorm2;
54 RealScalar residualNorm2 = residual.squaredNorm();
55 if (residualNorm2 < threshold)
58 tol_error =
sqrt(residualNorm2 / rhsNorm2);
63 p = precond.solve(residual);
65 VectorType z(n), tmp(n);
70 tmp.noalias() = mat * p;
72 Scalar alpha = absNew / p.dot(tmp);
74 residual -= alpha * tmp;
76 residualNorm2 = residual.squaredNorm();
77 if(residualNorm2 < threshold)
80 z = precond.solve(residual);
82 RealScalar absOld = absNew;
84 RealScalar beta = absNew / absOld;
88 tol_error =
sqrt(residualNorm2 / rhsNorm2);
94 template<
typename _MatrixType,
int _UpLo=
Lower,
100 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
157 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
161 using Base::mp_matrix;
163 using Base::m_iterations;
165 using Base::m_isInitialized;
168 typedef typename MatrixType::Scalar
Scalar;
169 typedef typename MatrixType::Index
Index;
201 template<
typename Rhs,
typename Guess>
205 eigen_assert(m_isInitialized &&
"ConjugateGradient is not initialized.");
207 &&
"ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
213 template<
typename Rhs,
typename Dest>
216 m_iterations = Base::maxIterations();
217 m_error = Base::m_tolerance;
219 for(
int j=0; j<b.cols(); ++j)
221 m_iterations = Base::maxIterations();
222 m_error = Base::m_tolerance;
226 Base::m_preconditioner, m_iterations, m_error);
229 m_isInitialized =
true;
234 template<
typename Rhs,
typename Dest>
238 _solveWithGuess(b,x);
248 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner,
typename Rhs>
255 template<typename Dest>
void evalTo(Dest& dst)
const 257 dec()._solve(
rhs(),dst);
265 #endif // EIGEN_CONJUGATE_GRADIENT_H A preconditioner based on the digonal entries.
IntermediateState sqrt(const Expression &arg)
ConjugateGradient< _MatrixType, _UpLo, _Preconditioner > Dec
EIGEN_DONT_INLINE void conjugate_gradient(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, int &iters, typename Dest::RealScalar &tol_error)
ConjugateGradient(const MatrixType &A)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Block< Derived, internal::traits< Derived >::RowsAtCompileTime, 1,!IsRowMajor > ColXpr
A conjugate gradient solver for sparse self-adjoint problems.
const internal::solve_retval_with_guess< ConjugateGradient, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
RealReturnType real() const
_Preconditioner Preconditioner
void _solveWithGuess(const Rhs &b, Dest &x) const
void _solve(const Rhs &b, Dest &x) const
void rhs(const real_t *x, real_t *f)
MatrixType::RealScalar RealScalar
MatrixType::Scalar Scalar
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
IterativeSolverBase< ConjugateGradient > Base
#define EIGEN_DONT_INLINE
EIGEN_STRONG_INLINE Index cols() const
Base class for linear iterative solvers.
Base class for all dense matrices, vectors, and expressions.
_Preconditioner Preconditioner