src/IterativeLinearSolvers/ConjugateGradient.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_CONJUGATE_GRADIENT_H
0011 #define EIGEN_CONJUGATE_GRADIENT_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016
0017 /** \internal Low-level conjugate gradient algorithm
0018   * \param mat The matrix A
0019   * \param rhs The right hand side vector b
0020   * \param x On input and initial solution, on output the computed solution.
0021   * \param precond A preconditioner being able to efficiently solve for an
0022   *                approximation of Ax=b (regardless of b)
0023   * \param iters On input the max number of iteration, on output the number of performed iterations.
0024   * \param tol_error On input the tolerance error, on output an estimation of the relative error.
0025   */
0026 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
0027 EIGEN_DONT_INLINE
0028 void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
0029                         const Preconditioner& precond, Index& iters,
0030                         typename Dest::RealScalar& tol_error)
0031 {
0032   using std::sqrt;
0033   using std::abs;
0034   typedef typename Dest::RealScalar RealScalar;
0035   typedef typename Dest::Scalar Scalar;
0036   typedef Matrix<Scalar,Dynamic,1> VectorType;
0037
0038   RealScalar tol = tol_error;
0039   Index maxIters = iters;
0040
0041   Index n = mat.cols();
0042
0043   VectorType residual = rhs - mat * x; //initial residual
0044
0045   RealScalar rhsNorm2 = rhs.squaredNorm();
0046   if(rhsNorm2 == 0)
0047   {
0048     x.setZero();
0049     iters = 0;
0050     tol_error = 0;
0051     return;
0052   }
0053   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
0054   RealScalar threshold = numext::maxi(RealScalar(tol*tol*rhsNorm2),considerAsZero);
0055   RealScalar residualNorm2 = residual.squaredNorm();
0056   if (residualNorm2 < threshold)
0057   {
0058     iters = 0;
0059     tol_error = sqrt(residualNorm2 / rhsNorm2);
0060     return;
0061   }
0062
0063   VectorType p(n);
0064   p = precond.solve(residual);      // initial search direction
0065
0066   VectorType z(n), tmp(n);
0067   RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
0068   Index i = 0;
0069   while(i < maxIters)
0070   {
0071     tmp.noalias() = mat * p;                    // the bottleneck of the algorithm
0072
0073     Scalar alpha = absNew / p.dot(tmp);         // the amount we travel on dir
0074     x += alpha * p;                             // update solution
0075     residual -= alpha * tmp;                    // update residual
0076
0077     residualNorm2 = residual.squaredNorm();
0078     if(residualNorm2 < threshold)
0079       break;
0080
0081     z = precond.solve(residual);                // approximately solve for "A z = residual"
0082
0083     RealScalar absOld = absNew;
0084     absNew = numext::real(residual.dot(z));     // update the absolute value of r
0085     RealScalar beta = absNew / absOld;          // calculate the Gram-Schmidt value used to create the new search direction
0086     p = z + beta * p;                           // update search direction
0087     i++;
0088   }
0089   tol_error = sqrt(residualNorm2 / rhsNorm2);
0090   iters = i;
0091 }
0092
0093 }
0094
0095 template< typename _MatrixType, int _UpLo=Lower,
0096           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
0097 class ConjugateGradient;
0098
0099 namespace internal {
0100
0101 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
0102 struct traits<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
0103 {
0104   typedef _MatrixType MatrixType;
0105   typedef _Preconditioner Preconditioner;
0106 };
0107
0108 }
0109
0110 /** \ingroup IterativeLinearSolvers_Module
0111   * \brief A conjugate gradient solver for sparse (or dense) self-adjoint problems
0112   *
0113   * This class allows to solve for A.x = b linear problems using an iterative conjugate gradient algorithm.
0114   * The matrix A must be selfadjoint. The matrix A and the vectors x and b can be either dense or sparse.
0115   *
0116   * \tparam _MatrixType the type of the matrix A, can be a dense or a sparse matrix.
0117   * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower,
0118   *               \c Upper, or \c Lower|Upper in which the full matrix entries will be considered.
0119   *               Default is \c Lower, best performance is \c Lower|Upper.
0120   * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
0121   *
0122   * \implsparsesolverconcept
0123   *
0124   * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
0125   * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
0126   * and NumTraits<Scalar>::epsilon() for the tolerance.
0127   *
0128   * The tolerance corresponds to the relative residual error: |Ax-b|/|b|
0129   *
0130   * \b Performance: Even though the default value of \c _UpLo is \c Lower, significantly higher performance is
0131   * achieved when using a complete matrix and \b Lower|Upper as the \a _UpLo template parameter. Moreover, in this
0132   * case multi-threading can be exploited if the user code is compiled with OpenMP enabled.
0133   * See \ref TopicMultiThreading for details.
0134   *
0135   * This class can be used as the direct solver classes. Here is a typical usage example:
0136     \code
0137     int n = 10000;
0138     VectorXd x(n), b(n);
0139     SparseMatrix<double> A(n,n);
0140     // fill A and b
0141     ConjugateGradient<SparseMatrix<double>, Lower|Upper> cg;
0142     cg.compute(A);
0143     x = cg.solve(b);
0144     std::cout << "#iterations:     " << cg.iterations() << std::endl;
0145     std::cout << "estimated error: " << cg.error()      << std::endl;
0146     // update b, and solve again
0147     x = cg.solve(b);
0148     \endcode
0149   *
0150   * By default the iterations start with x=0 as an initial guess of the solution.
0151   * One can control the start using the solveWithGuess() method.
0152   *
0153   * ConjugateGradient can also be used in a matrix-free context, see the following \link MatrixfreeSolverExample example \endlink.
0154   *
0155   * \sa class LeastSquaresConjugateGradient, class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
0156   */
0157 template< typename _MatrixType, int _UpLo, typename _Preconditioner>
0158 class ConjugateGradient : public IterativeSolverBase<ConjugateGradient<_MatrixType,_UpLo,_Preconditioner> >
0159 {
0160   typedef IterativeSolverBase<ConjugateGradient> Base;
0161   using Base::matrix;
0162   using Base::m_error;
0163   using Base::m_iterations;
0164   using Base::m_info;
0165   using Base::m_isInitialized;
0166 public:
0167   typedef _MatrixType MatrixType;
0168   typedef typename MatrixType::Scalar Scalar;
0169   typedef typename MatrixType::RealScalar RealScalar;
0170   typedef _Preconditioner Preconditioner;
0171
0172   enum {
0173     UpLo = _UpLo
0174   };
0175
0176 public:
0177
0178   /** Default constructor. */
0179   ConjugateGradient() : Base() {}
0180
0181   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
0182     *
0183     * This constructor is a shortcut for the default constructor followed
0184     * by a call to compute().
0185     *
0186     * \warning this class stores a reference to the matrix A as well as some
0187     * precomputed values that depend on it. Therefore, if \a A is changed
0188     * this class becomes invalid. Call compute() to update it with the new
0189     * matrix A, or modify a copy of A.
0190     */
0191   template<typename MatrixDerived>
0192   explicit ConjugateGradient(const EigenBase<MatrixDerived>& A) : Base(A.derived()) {}
0193
0194   ~ConjugateGradient() {}
0195
0196   /** \internal */
0197   template<typename Rhs,typename Dest>
0198   void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
0199   {
0200     typedef typename Base::MatrixWrapper MatrixWrapper;
0201     typedef typename Base::ActualMatrixType ActualMatrixType;
0202     enum {
0203       TransposeInput  =   (!MatrixWrapper::MatrixFree)
0204                       &&  (UpLo==(Lower|Upper))
0205                       &&  (!MatrixType::IsRowMajor)
0206                       &&  (!NumTraits<Scalar>::IsComplex)
0207     };
0208     typedef typename internal::conditional<TransposeInput,Transpose<const ActualMatrixType>, ActualMatrixType const&>::type RowMajorWrapper;
0209     EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(MatrixWrapper::MatrixFree,UpLo==(Lower|Upper)),MATRIX_FREE_CONJUGATE_GRADIENT_IS_COMPATIBLE_WITH_UPPER_UNION_LOWER_MODE_ONLY);
0210     typedef typename internal::conditional<UpLo==(Lower|Upper),
0211                                            RowMajorWrapper,
0212                                            typename MatrixWrapper::template ConstSelfAdjointViewReturnType<UpLo>::Type
0213                                           >::type SelfAdjointWrapper;
0214
0215     m_iterations = Base::maxIterations();
0216     m_error = Base::m_tolerance;
0217
0218     RowMajorWrapper row_mat(matrix());
0219     internal::conjugate_gradient(SelfAdjointWrapper(row_mat), b, x, Base::m_preconditioner, m_iterations, m_error);
0220     m_info = m_error <= Base::m_tolerance ? Success : NoConvergence;
0221   }
0222
0223 protected:
0224
0225 };
0226
0227 } // end namespace Eigen
0228
0229 #endif // EIGEN_CONJUGATE_GRADIENT_H