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_BASIC_PRECONDITIONERS_H
     11 #define EIGEN_BASIC_PRECONDITIONERS_H
     12 
     13 namespace Eigen {
     14 
     15 /** \ingroup IterativeLinearSolvers_Module
     16   * \brief A preconditioner based on the digonal entries
     17   *
     18   * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.
     19   * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:
     20   * \code
     21   * A.diagonal().asDiagonal() . x = b
     22   * \endcode
     23   *
     24   * \tparam _Scalar the type of the scalar.
     25   *
     26   * This preconditioner is suitable for both selfadjoint and general problems.
     27   * The diagonal entries are pre-inverted and stored into a dense vector.
     28   *
     29   * \note A variant that has yet to be implemented would attempt to preserve the norm of each column.
     30   *
     31   */
     32 template <typename _Scalar>
     33 class DiagonalPreconditioner
     34 {
     35     typedef _Scalar Scalar;
     36     typedef Matrix<Scalar,Dynamic,1> Vector;
     37     typedef typename Vector::Index Index;
     38 
     39   public:
     40     // this typedef is only to export the scalar type and compile-time dimensions to solve_retval
     41     typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;
     42 
     43     DiagonalPreconditioner() : m_isInitialized(false) {}
     44 
     45     template<typename MatType>
     46     DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
     47     {
     48       compute(mat);
     49     }
     50 
     51     Index rows() const { return m_invdiag.size(); }
     52     Index cols() const { return m_invdiag.size(); }
     53 
     54     template<typename MatType>
     55     DiagonalPreconditioner& analyzePattern(const MatType& )
     56     {
     57       return *this;
     58     }
     59 
     60     template<typename MatType>
     61     DiagonalPreconditioner& factorize(const MatType& mat)
     62     {
     63       m_invdiag.resize(mat.cols());
     64       for(int j=0; j<mat.outerSize(); ++j)
     65       {
     66         typename MatType::InnerIterator it(mat,j);
     67         while(it && it.index()!=j) ++it;
     68         if(it && it.index()==j)
     69           m_invdiag(j) = Scalar(1)/it.value();
     70         else
     71           m_invdiag(j) = 0;
     72       }
     73       m_isInitialized = true;
     74       return *this;
     75     }
     76 
     77     template<typename MatType>
     78     DiagonalPreconditioner& compute(const MatType& mat)
     79     {
     80       return factorize(mat);
     81     }
     82 
     83     template<typename Rhs, typename Dest>
     84     void _solve(const Rhs& b, Dest& x) const
     85     {
     86       x = m_invdiag.array() * b.array() ;
     87     }
     88 
     89     template<typename Rhs> inline const internal::solve_retval<DiagonalPreconditioner, Rhs>
     90     solve(const MatrixBase<Rhs>& b) const
     91     {
     92       eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
     93       eigen_assert(m_invdiag.size()==b.rows()
     94                 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
     95       return internal::solve_retval<DiagonalPreconditioner, Rhs>(*this, b.derived());
     96     }
     97 
     98   protected:
     99     Vector m_invdiag;
    100     bool m_isInitialized;
    101 };
    102 
    103 namespace internal {
    104 
    105 template<typename _MatrixType, typename Rhs>
    106 struct solve_retval<DiagonalPreconditioner<_MatrixType>, Rhs>
    107   : solve_retval_base<DiagonalPreconditioner<_MatrixType>, Rhs>
    108 {
    109   typedef DiagonalPreconditioner<_MatrixType> Dec;
    110   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
    111 
    112   template<typename Dest> void evalTo(Dest& dst) const
    113   {
    114     dec()._solve(rhs(),dst);
    115   }
    116 };
    117 
    118 }
    119 
    120 /** \ingroup IterativeLinearSolvers_Module
    121   * \brief A naive preconditioner which approximates any matrix as the identity matrix
    122   *
    123   * \sa class DiagonalPreconditioner
    124   */
    125 class IdentityPreconditioner
    126 {
    127   public:
    128 
    129     IdentityPreconditioner() {}
    130 
    131     template<typename MatrixType>
    132     IdentityPreconditioner(const MatrixType& ) {}
    133 
    134     template<typename MatrixType>
    135     IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
    136 
    137     template<typename MatrixType>
    138     IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
    139 
    140     template<typename MatrixType>
    141     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
    142 
    143     template<typename Rhs>
    144     inline const Rhs& solve(const Rhs& b) const { return b; }
    145 };
    146 
    147 } // end namespace Eigen
    148 
    149 #endif // EIGEN_BASIC_PRECONDITIONERS_H
    150