Home | History | Annotate | Download | only in IterativeLinearSolvers
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_ITERATIVE_SOLVER_BASE_H
     11 #define EIGEN_ITERATIVE_SOLVER_BASE_H
     12 
     13 namespace Eigen {
     14 
     15 /** \ingroup IterativeLinearSolvers_Module
     16   * \brief Base class for linear iterative solvers
     17   *
     18   * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
     19   */
     20 template< typename Derived>
     21 class IterativeSolverBase : internal::noncopyable
     22 {
     23 public:
     24   typedef typename internal::traits<Derived>::MatrixType MatrixType;
     25   typedef typename internal::traits<Derived>::Preconditioner Preconditioner;
     26   typedef typename MatrixType::Scalar Scalar;
     27   typedef typename MatrixType::Index Index;
     28   typedef typename MatrixType::RealScalar RealScalar;
     29 
     30 public:
     31 
     32   Derived& derived() { return *static_cast<Derived*>(this); }
     33   const Derived& derived() const { return *static_cast<const Derived*>(this); }
     34 
     35   /** Default constructor. */
     36   IterativeSolverBase()
     37     : mp_matrix(0)
     38   {
     39     init();
     40   }
     41 
     42   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
     43     *
     44     * This constructor is a shortcut for the default constructor followed
     45     * by a call to compute().
     46     *
     47     * \warning this class stores a reference to the matrix A as well as some
     48     * precomputed values that depend on it. Therefore, if \a A is changed
     49     * this class becomes invalid. Call compute() to update it with the new
     50     * matrix A, or modify a copy of A.
     51     */
     52   IterativeSolverBase(const MatrixType& A)
     53   {
     54     init();
     55     compute(A);
     56   }
     57 
     58   ~IterativeSolverBase() {}
     59 
     60   /** Initializes the iterative solver for the sparcity pattern of the matrix \a A for further solving \c Ax=b problems.
     61     *
     62     * Currently, this function mostly call analyzePattern on the preconditioner. In the future
     63     * we might, for instance, implement column reodering for faster matrix vector products.
     64     */
     65   Derived& analyzePattern(const MatrixType& A)
     66   {
     67     m_preconditioner.analyzePattern(A);
     68     m_isInitialized = true;
     69     m_analysisIsOk = true;
     70     m_info = Success;
     71     return derived();
     72   }
     73 
     74   /** Initializes the iterative solver with the numerical values of the matrix \a A for further solving \c Ax=b problems.
     75     *
     76     * Currently, this function mostly call factorize on the preconditioner.
     77     *
     78     * \warning this class stores a reference to the matrix A as well as some
     79     * precomputed values that depend on it. Therefore, if \a A is changed
     80     * this class becomes invalid. Call compute() to update it with the new
     81     * matrix A, or modify a copy of A.
     82     */
     83   Derived& factorize(const MatrixType& A)
     84   {
     85     eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
     86     mp_matrix = &A;
     87     m_preconditioner.factorize(A);
     88     m_factorizationIsOk = true;
     89     m_info = Success;
     90     return derived();
     91   }
     92 
     93   /** Initializes the iterative solver with the matrix \a A for further solving \c Ax=b problems.
     94     *
     95     * Currently, this function mostly initialized/compute the preconditioner. In the future
     96     * we might, for instance, implement column reodering for faster matrix vector products.
     97     *
     98     * \warning this class stores a reference to the matrix A as well as some
     99     * precomputed values that depend on it. Therefore, if \a A is changed
    100     * this class becomes invalid. Call compute() to update it with the new
    101     * matrix A, or modify a copy of A.
    102     */
    103   Derived& compute(const MatrixType& A)
    104   {
    105     mp_matrix = &A;
    106     m_preconditioner.compute(A);
    107     m_isInitialized = true;
    108     m_analysisIsOk = true;
    109     m_factorizationIsOk = true;
    110     m_info = Success;
    111     return derived();
    112   }
    113 
    114   /** \internal */
    115   Index rows() const { return mp_matrix ? mp_matrix->rows() : 0; }
    116   /** \internal */
    117   Index cols() const { return mp_matrix ? mp_matrix->cols() : 0; }
    118 
    119   /** \returns the tolerance threshold used by the stopping criteria */
    120   RealScalar tolerance() const { return m_tolerance; }
    121 
    122   /** Sets the tolerance threshold used by the stopping criteria */
    123   Derived& setTolerance(const RealScalar& tolerance)
    124   {
    125     m_tolerance = tolerance;
    126     return derived();
    127   }
    128 
    129   /** \returns a read-write reference to the preconditioner for custom configuration. */
    130   Preconditioner& preconditioner() { return m_preconditioner; }
    131 
    132   /** \returns a read-only reference to the preconditioner. */
    133   const Preconditioner& preconditioner() const { return m_preconditioner; }
    134 
    135   /** \returns the max number of iterations */
    136   int maxIterations() const
    137   {
    138     return (mp_matrix && m_maxIterations<0) ? mp_matrix->cols() : m_maxIterations;
    139   }
    140 
    141   /** Sets the max number of iterations */
    142   Derived& setMaxIterations(int maxIters)
    143   {
    144     m_maxIterations = maxIters;
    145     return derived();
    146   }
    147 
    148   /** \returns the number of iterations performed during the last solve */
    149   int iterations() const
    150   {
    151     eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
    152     return m_iterations;
    153   }
    154 
    155   /** \returns the tolerance error reached during the last solve */
    156   RealScalar error() const
    157   {
    158     eigen_assert(m_isInitialized && "ConjugateGradient is not initialized.");
    159     return m_error;
    160   }
    161 
    162   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
    163     *
    164     * \sa compute()
    165     */
    166   template<typename Rhs> inline const internal::solve_retval<Derived, Rhs>
    167   solve(const MatrixBase<Rhs>& b) const
    168   {
    169     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
    170     eigen_assert(rows()==b.rows()
    171               && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
    172     return internal::solve_retval<Derived, Rhs>(derived(), b.derived());
    173   }
    174 
    175   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
    176     *
    177     * \sa compute()
    178     */
    179   template<typename Rhs>
    180   inline const internal::sparse_solve_retval<IterativeSolverBase, Rhs>
    181   solve(const SparseMatrixBase<Rhs>& b) const
    182   {
    183     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
    184     eigen_assert(rows()==b.rows()
    185               && "IterativeSolverBase::solve(): invalid number of rows of the right hand side matrix b");
    186     return internal::sparse_solve_retval<IterativeSolverBase, Rhs>(*this, b.derived());
    187   }
    188 
    189   /** \returns Success if the iterations converged, and NoConvergence otherwise. */
    190   ComputationInfo info() const
    191   {
    192     eigen_assert(m_isInitialized && "IterativeSolverBase is not initialized.");
    193     return m_info;
    194   }
    195 
    196   /** \internal */
    197   template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex>
    198   void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
    199   {
    200     eigen_assert(rows()==b.rows());
    201 
    202     int rhsCols = b.cols();
    203     int size = b.rows();
    204     Eigen::Matrix<DestScalar,Dynamic,1> tb(size);
    205     Eigen::Matrix<DestScalar,Dynamic,1> tx(size);
    206     for(int k=0; k<rhsCols; ++k)
    207     {
    208       tb = b.col(k);
    209       tx = derived().solve(tb);
    210       dest.col(k) = tx.sparseView(0);
    211     }
    212   }
    213 
    214 protected:
    215   void init()
    216   {
    217     m_isInitialized = false;
    218     m_analysisIsOk = false;
    219     m_factorizationIsOk = false;
    220     m_maxIterations = -1;
    221     m_tolerance = NumTraits<Scalar>::epsilon();
    222   }
    223   const MatrixType* mp_matrix;
    224   Preconditioner m_preconditioner;
    225 
    226   int m_maxIterations;
    227   RealScalar m_tolerance;
    228 
    229   mutable RealScalar m_error;
    230   mutable int m_iterations;
    231   mutable ComputationInfo m_info;
    232   mutable bool m_isInitialized, m_analysisIsOk, m_factorizationIsOk;
    233 };
    234 
    235 namespace internal {
    236 
    237 template<typename Derived, typename Rhs>
    238 struct sparse_solve_retval<IterativeSolverBase<Derived>, Rhs>
    239   : sparse_solve_retval_base<IterativeSolverBase<Derived>, Rhs>
    240 {
    241   typedef IterativeSolverBase<Derived> Dec;
    242   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
    243 
    244   template<typename Dest> void evalTo(Dest& dst) const
    245   {
    246     dec().derived()._solve_sparse(rhs(),dst);
    247   }
    248 };
    249 
    250 } // end namespace internal
    251 
    252 } // end namespace Eigen
    253 
    254 #endif // EIGEN_ITERATIVE_SOLVER_BASE_H
    255