10 #ifndef EIGEN_CONJUGATE_GRADIENT_H
11 #define EIGEN_CONJUGATE_GRADIENT_H
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,
int& iters,
30 typename Dest::RealScalar& tol_error)
34 typedef typename Dest::RealScalar RealScalar;
35 typedef typename Dest::Scalar Scalar;
36 typedef Matrix<Scalar,Dynamic,1> VectorType;
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);
66 RealScalar absNew = numext::real(residual.dot(p));
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;
83 absNew = numext::real(residual.dot(z));
84 RealScalar beta = absNew / absOld;
88 tol_error = sqrt(residualNorm2 / rhsNorm2);
94 template<
typename _MatrixType,
int _UpLo=
Lower,
95 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
100 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
103 typedef _MatrixType MatrixType;
104 typedef _Preconditioner Preconditioner;
144 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
145 class ConjugateGradient :
public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
147 typedef IterativeSolverBase<ConjugateGradient> Base;
148 using Base::mp_matrix;
150 using Base::m_iterations;
152 using Base::m_isInitialized;
154 typedef _MatrixType MatrixType;
155 typedef typename MatrixType::Scalar Scalar;
156 typedef typename MatrixType::Index Index;
157 typedef typename MatrixType::RealScalar RealScalar;
158 typedef _Preconditioner Preconditioner;
188 template<
typename Rhs,
typename Guess>
189 inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
192 eigen_assert(m_isInitialized &&
"ConjugateGradient is not initialized.");
193 eigen_assert(Base::rows()==b.rows()
194 &&
"ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
195 return internal::solve_retval_with_guess
200 template<
typename Rhs,
typename Dest>
201 void _solveWithGuess(
const Rhs& b, Dest& x)
const
203 typedef typename internal::conditional<UpLo==(
Lower|
Upper),
206 >::type MatrixWrapperType;
208 m_error = Base::m_tolerance;
210 for(
int j=0; j<b.cols(); ++j)
213 m_error = Base::m_tolerance;
215 typename Dest::ColXpr xj(x,j);
216 internal::conjugate_gradient(MatrixWrapperType(*mp_matrix), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
219 m_isInitialized =
true;
224 template<
typename Rhs,
typename Dest>
225 void _solve(
const Rhs& b, Dest& x)
const
228 _solveWithGuess(b,x);
238 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner,
typename Rhs>
239 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
240 : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
242 typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
243 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
245 template<typename Dest>
void evalTo(Dest& dst)
const
247 dec()._solve(rhs(),dst);
255 #endif // EIGEN_CONJUGATE_GRADIENT_H
Definition: Constants.h:167
ConjugateGradient(const MatrixType &A)
Definition: ConjugateGradient.h:179
ConjugateGradient()
Definition: ConjugateGradient.h:167
const internal::solve_retval_with_guess< ConjugateGradient, Rhs, Guess > solveWithGuess(const MatrixBase< Rhs > &b, const Guess &x0) const
Definition: ConjugateGradient.h:190
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:49
A conjugate gradient solver for sparse self-adjoint problems.
Definition: ConjugateGradient.h:96
Definition: Constants.h:169
Definition: Constants.h:380
Definition: Constants.h:376
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:21
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
int maxIterations() const
Definition: IterativeSolverBase.h:136