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 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
54 RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero);
55 RealScalar residualNorm2 = residual.squaredNorm();
56 if (residualNorm2 < threshold)
59 tol_error =
sqrt(residualNorm2 / rhsNorm2);
64 p = precond.solve(residual);
66 VectorType z(n), tmp(n);
67 RealScalar absNew = numext::real(residual.dot(p));
71 tmp.noalias() = mat * p;
73 Scalar alpha = absNew / p.dot(tmp);
75 residual -= alpha * tmp;
77 residualNorm2 = residual.squaredNorm();
78 if(residualNorm2 < threshold)
81 z = precond.solve(residual);
83 RealScalar absOld = absNew;
84 absNew = numext::real(residual.dot(z));
85 RealScalar beta = absNew / absOld;
89 tol_error =
sqrt(residualNorm2 / rhsNorm2);
95 template<
typename _MatrixType,
int _UpLo=
Lower,
96 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
101 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
104 typedef _MatrixType MatrixType;
105 typedef _Preconditioner Preconditioner;
157 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
163 using Base::m_iterations;
165 using Base::m_isInitialized;
167 typedef _MatrixType MatrixType;
168 typedef typename MatrixType::Scalar Scalar;
169 typedef typename MatrixType::RealScalar RealScalar;
170 typedef _Preconditioner Preconditioner;
191 template<
typename MatrixDerived>
197 template<
typename Rhs,
typename Dest>
198 void _solve_vector_with_guess_impl(
const Rhs& b, Dest& x)
const 200 typedef typename Base::MatrixWrapper MatrixWrapper;
201 typedef typename Base::ActualMatrixType ActualMatrixType;
203 TransposeInput = (!MatrixWrapper::MatrixFree)
205 && (!MatrixType::IsRowMajor)
208 typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType
const&>::type RowMajorWrapper;
209 EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(
Lower|
Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
210 typedef typename internal::conditional<UpLo==(
Lower|
Upper),
212 typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
213 >::type SelfAdjointWrapper;
215 m_iterations = Base::maxIterations();
216 m_error = Base::m_tolerance;
218 RowMajorWrapper row_mat(matrix());
219 internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
229 #endif // EIGEN_CONJUGATE_GRADIENT_H const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
Definition: Core:141
A conjugate gradient solver for sparse (or dense) self-adjoint problems.
Definition: ConjugateGradient.h:97
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
ConjugateGradient()
Definition: ConjugateGradient.h:179
Definition: EigenBase.h:29
Definition: Constants.h:209
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: Constants.h:211
Definition: Constants.h:442
Definition: Eigen_Colamd.h:50
ConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: ConjugateGradient.h:192
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:143
Definition: Constants.h:446