00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_CONJUGATE_GRADIENT_H
00026 #define EIGEN_CONJUGATE_GRADIENT_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00041 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
00042 EIGEN_DONT_INLINE
00043 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
00044 const Preconditioner& precond, int& iters,
00045 typename Dest::RealScalar& tol_error)
00046 {
00047 using std::sqrt;
00048 using std::abs;
00049 typedef typename Dest::RealScalar RealScalar;
00050 typedef typename Dest::Scalar Scalar;
00051 typedef Matrix<Scalar,Dynamic,1> VectorType;
00052
00053 RealScalar tol = tol_error;
00054 int maxIters = iters;
00055
00056 int n = mat.cols();
00057
00058 VectorType residual = rhs - mat * x;
00059 VectorType p(n);
00060
00061 p = precond.solve(residual);
00062
00063 VectorType z(n), tmp(n);
00064 RealScalar absNew = internal::real(residual.dot(p));
00065 RealScalar rhsNorm2 = rhs.squaredNorm();
00066 RealScalar residualNorm2 = 0;
00067 RealScalar threshold = tol*tol*rhsNorm2;
00068 int i = 0;
00069 while(i < maxIters)
00070 {
00071 tmp.noalias() = mat * p;
00072
00073 Scalar alpha = absNew / p.dot(tmp);
00074 x += alpha * p;
00075 residual -= alpha * tmp;
00076
00077 residualNorm2 = residual.squaredNorm();
00078 if(residualNorm2 < threshold)
00079 break;
00080
00081 z = precond.solve(residual);
00082
00083 RealScalar absOld = absNew;
00084 absNew = internal::real(residual.dot(z));
00085 RealScalar beta = absNew / absOld;
00086 p = z + beta * p;
00087 i++;
00088 }
00089 tol_error = sqrt(residualNorm2 / rhsNorm2);
00090 iters = i;
00091 }
00092
00093 }
00094
00095 template< typename _MatrixType, int _UpLo=Lower,
00096 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
00097 class ConjugateGradient;
00098
00099 namespace internal {
00100
00101 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
00102 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
00103 {
00104 typedef _MatrixType MatrixType;
00105 typedef _Preconditioner Preconditioner;
00106 };
00107
00108 }
00109
00158 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
00159 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
00160 {
00161 typedef IterativeSolverBase<ConjugateGradient> Base;
00162 using Base::mp_matrix;
00163 using Base::m_error;
00164 using Base::m_iterations;
00165 using Base::m_info;
00166 using Base::m_isInitialized;
00167 public:
00168 typedef _MatrixType MatrixType;
00169 typedef typename MatrixType::Scalar Scalar;
00170 typedef typename MatrixType::Index Index;
00171 typedef typename MatrixType::RealScalar RealScalar;
00172 typedef _Preconditioner Preconditioner;
00173
00174 enum {
00175 UpLo = _UpLo
00176 };
00177
00178 public:
00179
00181 ConjugateGradient() : Base() {}
00182
00193 ConjugateGradient(const MatrixType& A) : Base(A) {}
00194
00195 ~ConjugateGradient() {}
00196
00202 template<typename Rhs,typename Guess>
00203 inline const internal::solve_retval_with_guess<ConjugateGradient, Rhs, Guess>
00204 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
00205 {
00206 eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
00207 eigen_assert(Base::rows()==b.rows()
00208 && "ConjugateGradient::solve(): invalid number of rows of the right hand side matrix b");
00209 return internal::solve_retval_with_guess
00210 <ConjugateGradient, Rhs, Guess>(*this, b.derived(), x0);
00211 }
00212
00214 template<typename Rhs,typename Dest>
00215 void _solveWithGuess(const Rhs& b, Dest& x) const
00216 {
00217 m_iterations = Base::maxIterations();
00218 m_error = Base::m_tolerance;
00219
00220 for(int j=0; j<b.cols(); ++j)
00221 {
00222 m_iterations = Base::maxIterations();
00223 m_error = Base::m_tolerance;
00224
00225 typename Dest::ColXpr xj(x,j);
00226 internal::conjugate_gradient(mp_matrix->template selfadjointView<UpLo>(), b.col(j), xj,
00227 Base::m_preconditioner, m_iterations, m_error);
00228 }
00229
00230 m_isInitialized = true;
00231 m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
00232 }
00233
00235 template<typename Rhs,typename Dest>
00236 void _solve(const Rhs& b, Dest& x) const
00237 {
00238 x.setOnes();
00239 _solveWithGuess(b,x);
00240 }
00241
00242 protected:
00243
00244 };
00245
00246
00247 namespace internal {
00248
00249 template<typename _MatrixType, int _UpLo, typename _Preconditioner, typename Rhs>
00250 struct solve_retval<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
00251 : solve_retval_base<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner>, Rhs>
00252 {
00253 typedef ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> Dec;
00254 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00255
00256 template<typename Dest> void evalTo(Dest& dst) const
00257 {
00258 dec()._solve(rhs(),dst);
00259 }
00260 };
00261
00262 }
00263
00264 }
00265
00266 #endif // EIGEN_CONJUGATE_GRADIENT_H