ConjugateGradient.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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; //initial residual
00059   VectorType p(n);
00060 
00061   p = precond.solve(residual);      //initial search direction
00062 
00063   VectorType z(n), tmp(n);
00064   RealScalar absNew = internal::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
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;              // the bottleneck of the algorithm
00072 
00073     Scalar alpha = absNew / p.dot(tmp);   // the amount we travel on dir
00074     x += alpha * p;                       // update solution
00075     residual -= alpha * tmp;              // update residue
00076     
00077     residualNorm2 = residual.squaredNorm();
00078     if(residualNorm2 < threshold)
00079       break;
00080     
00081     z = precond.solve(residual);          // approximately solve for "A z = residual"
00082 
00083     RealScalar absOld = absNew;
00084     absNew = internal::real(residual.dot(z));     // update the absolute value of r
00085     RealScalar beta = absNew / absOld;            // calculate the Gram-Schmidt value used to create the new search direction
00086     p = z + beta * p;                             // update search direction
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 } // end namespace internal
00263 
00264 } // end namespace Eigen
00265 
00266 #endif // EIGEN_CONJUGATE_GRADIENT_H