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, Index& 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;
39 Index maxIters = iters;
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;
156 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
157 class ConjugateGradient :
public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
159 typedef IterativeSolverBase<ConjugateGradient> Base;
162 using Base::m_iterations;
164 using Base::m_isInitialized;
166 typedef _MatrixType MatrixType;
167 typedef typename MatrixType::Scalar Scalar;
168 typedef typename MatrixType::RealScalar RealScalar;
169 typedef _Preconditioner Preconditioner;
190 template<
typename MatrixDerived>
196 template<
typename Rhs,
typename Dest>
197 void _solve_with_guess_impl(
const Rhs& b, Dest& x)
const
199 typedef typename Base::MatrixWrapper MatrixWrapper;
200 typedef typename Base::ActualMatrixType ActualMatrixType;
202 TransposeInput = (!MatrixWrapper::MatrixFree)
204 && (!MatrixType::IsRowMajor)
205 && (!NumTraits<Scalar>::IsComplex)
207 typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType
const&>::type RowMajorWrapper;
208 EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(
Lower|
Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
209 typedef typename internal::conditional<UpLo==(
Lower|
Upper),
211 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
212 >::type SelfAdjointWrapper;
214 m_error = Base::m_tolerance;
216 for(Index j=0; j<b.cols(); ++j)
219 m_error = Base::m_tolerance;
221 typename Dest::ColXpr xj(x,j);
222 RowMajorWrapper row_mat(matrix());
223 internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
226 m_isInitialized =
true;
231 using Base::_solve_impl;
232 template<
typename Rhs,
typename Dest>
233 void _solve_impl(
const MatrixBase<Rhs>& b, Dest& x)
const
236 _solve_with_guess_impl(b.derived(),x);
245 #endif // EIGEN_CONJUGATE_GRADIENT_H
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:96
ConjugateGradient()
Definition: ConjugateGradient.h:178
Index maxIterations() const
Definition: IterativeSolverBase.h:281
Definition: EigenBase.h:28
Definition: Constants.h:204
Definition: Constants.h:432
Definition: Eigen_Colamd.h:54
ConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: ConjugateGradient.h:191
Definition: Constants.h:206
Definition: Constants.h:436