25 #ifndef EIGEN_CONJUGATE_GRADIENT_H
26 #define EIGEN_CONJUGATE_GRADIENT_H
41 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
44 const Preconditioner& precond,
int& iters,
45 typename Dest::RealScalar& tol_error)
49 typedef typename Dest::RealScalar RealScalar;
50 typedef typename Dest::Scalar Scalar;
53 RealScalar tol = tol_error;
58 VectorType residual = rhs - mat * x;
61 p = precond.solve(residual);
63 VectorType z(n), tmp(n);
65 RealScalar rhsNorm2 = rhs.squaredNorm();
66 RealScalar residualNorm2 = 0;
67 RealScalar threshold = tol*tol*rhsNorm2;
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;
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> >
97 class ConjugateGradient;
101 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
102 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
104 typedef _MatrixType MatrixType;
105 typedef _Preconditioner Preconditioner;
158 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner>
169 typedef typename MatrixType::Scalar
Scalar;
170 typedef typename MatrixType::Index
Index;
202 template<
typename Rhs,
typename Guess>
203 inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
208 &&
"ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
209 return internal::solve_retval_with_guess
214 template<
typename Rhs,
typename Dest>
220 for(
int j=0; j<b.cols(); ++j)
225 typename Dest::ColXpr xj(x,j);
235 template<
typename Rhs,
typename Dest>
249 template<
typename _MatrixType,
int _UpLo,
typename _Preconditioner,
typename Rhs>
250 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
251 : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
253 typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
256 template<typename Dest>
void evalTo(Dest& dst)
const
258 dec()._solve(rhs(),dst);
266 #endif // EIGEN_CONJUGATE_GRADIENT_H