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;
46 p = precond.solve(residual);
48 VectorType z(n), tmp(n);
49 RealScalar absNew = internal::real(residual.dot(p));
50 RealScalar rhsNorm2 = rhs.squaredNorm();
51 RealScalar residualNorm2 = 0;
52 RealScalar threshold = tol*tol*rhsNorm2;
56 tmp.noalias() = mat * p;
58 Scalar alpha = absNew / p.dot(tmp);
60 residual -= alpha * tmp;
62 residualNorm2 = residual.squaredNorm();
63 if(residualNorm2 < threshold)
66 z = precond.solve(residual);
68 RealScalar absOld = absNew;
69 absNew = internal::real(residual.dot(z));
70 RealScalar beta = absNew / absOld;
74 tol_error = sqrt(residualNorm2 / rhsNorm2);
80 template<
typename _MatrixType,
int _UpLo=
Lower,
81 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
86 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
89 typedef _MatrixType MatrixType;
90 typedef _Preconditioner Preconditioner;
143 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
144 class ConjugateGradient :
public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
146 typedef IterativeSolverBase<ConjugateGradient> Base;
147 using Base::mp_matrix;
149 using Base::m_iterations;
151 using Base::m_isInitialized;
153 typedef _MatrixType MatrixType;
154 typedef typename MatrixType::Scalar Scalar;
155 typedef typename MatrixType::Index Index;
156 typedef typename MatrixType::RealScalar RealScalar;
157 typedef _Preconditioner Preconditioner;
187 template<
typename Rhs,
typename Guess>
188 inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
191 eigen_assert(m_isInitialized &&
"ConjugateGradient is not initialized.");
192 eigen_assert(Base::rows()==b.rows()
193 &&
"ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
194 return internal::solve_retval_with_guess
199 template<
typename Rhs,
typename Dest>
200 void _solveWithGuess(
const Rhs& b, Dest& x)
const
203 m_error = Base::m_tolerance;
205 for(
int j=0; j<b.cols(); ++j)
208 m_error = Base::m_tolerance;
210 typename Dest::ColXpr xj(x,j);
211 internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
212 Base::m_preconditioner, m_iterations, m_error);
215 m_isInitialized =
true;
220 template<
typename Rhs,
typename Dest>
221 void _solve(
const Rhs& b, Dest& x)
const
224 _solveWithGuess(b,x);
234 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner,
typename Rhs>
235 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
236 : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
238 typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
239 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
241 template<typename Dest>
void evalTo(Dest& dst)
const
243 dec()._solve(rhs(),dst);
251 #endif // EIGEN_CONJUGATE_GRADIENT_H