10 #ifndef EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 11 #define EIGEN_LEAST_SQUARE_CONJUGATE_GRADIENT_H 26 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
28 void least_square_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;
41 Index m = mat.rows(), n = mat.cols();
43 VectorType residual = rhs - mat * x;
44 VectorType normal_residual = mat.adjoint() * residual;
46 RealScalar rhsNorm2 = (mat.adjoint()*rhs).squaredNorm();
54 RealScalar threshold = tol*tol*rhsNorm2;
55 RealScalar residualNorm2 = normal_residual.squaredNorm();
56 if (residualNorm2 < threshold)
59 tol_error =
sqrt(residualNorm2 / rhsNorm2);
64 p = precond.solve(normal_residual);
66 VectorType z(n), tmp(m);
67 RealScalar absNew = numext::real(normal_residual.dot(p));
71 tmp.noalias() = mat * p;
73 Scalar alpha = absNew / tmp.squaredNorm();
75 residual -= alpha * tmp;
76 normal_residual = mat.adjoint() * residual;
78 residualNorm2 = normal_residual.squaredNorm();
79 if(residualNorm2 < threshold)
82 z = precond.solve(normal_residual);
84 RealScalar absOld = absNew;
85 absNew = numext::real(normal_residual.dot(z));
86 RealScalar beta = absNew / absOld;
90 tol_error =
sqrt(residualNorm2 / rhsNorm2);
96 template<
typename _MatrixType,
97 typename _Preconditioner = LeastSquareDiagonalPreconditioner<typename _MatrixType::Scalar> >
102 template<
typename _MatrixType,
typename _Preconditioner>
105 typedef _MatrixType MatrixType;
106 typedef _Preconditioner Preconditioner;
148 template<
typename _MatrixType,
typename _Preconditioner>
154 using Base::m_iterations;
156 using Base::m_isInitialized;
158 typedef _MatrixType MatrixType;
159 typedef typename MatrixType::Scalar Scalar;
160 typedef typename MatrixType::RealScalar RealScalar;
161 typedef _Preconditioner Preconditioner;
178 template<
typename MatrixDerived>
184 template<
typename Rhs,
typename Dest>
185 void _solve_with_guess_impl(
const Rhs& b, Dest& x)
const 187 m_iterations = Base::maxIterations();
188 m_error = Base::m_tolerance;
190 for(
Index j=0; j<b.cols(); ++j)
192 m_iterations = Base::maxIterations();
193 m_error = Base::m_tolerance;
195 typename Dest::ColXpr xj(x,j);
196 internal::least_square_conjugate_gradient(matrix(), b.col(j), xj, Base::m_preconditioner, m_iterations, m_error);
199 m_isInitialized =
true;
204 using Base::_solve_impl;
205 template<
typename Rhs,
typename Dest>
209 _solve_with_guess_impl(b.
derived(),x);
216 #endif // EIGEN_LEAST_SQUARE_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:287
A conjugate gradient solver for sparse (or dense) least-square problems.
Definition: LeastSquareConjugateGradient.h:98
Derived & derived()
Definition: EigenBase.h:45
Definition: EigenBase.h:29
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Definition: Constants.h:432
Definition: Eigen_Colamd.h:50
LeastSquaresConjugateGradient()
Definition: LeastSquareConjugateGradient.h:166
LeastSquaresConjugateGradient(const EigenBase< MatrixDerived > &A)
Definition: LeastSquareConjugateGradient.h:179
Base class for linear iterative solvers.
Definition: IterativeSolverBase.h:143
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: Constants.h:436