src/IterativeLinearSolvers/BasicPreconditioners.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_BASIC_PRECONDITIONERS_H
0011 #define EIGEN_BASIC_PRECONDITIONERS_H
0012
0013 namespace RivetEigen {
0014
0015 /** \ingroup IterativeLinearSolvers_Module
0016   * \brief A preconditioner based on the digonal entries
0017   *
0018   * This class allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.
0019   * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:
0020     \code
0021     A.diagonal().asDiagonal() . x = b
0022     \endcode
0023   *
0024   * \tparam _Scalar the type of the scalar.
0025   *
0026   * \implsparsesolverconcept
0027   *
0028   * This preconditioner is suitable for both selfadjoint and general problems.
0029   * The diagonal entries are pre-inverted and stored into a dense vector.
0030   *
0031   * \note A variant that has yet to be implemented would attempt to preserve the norm of each column.
0032   *
0033   * \sa class LeastSquareDiagonalPreconditioner, class ConjugateGradient
0034   */
0035 template <typename _Scalar>
0036 class DiagonalPreconditioner
0037 {
0038     typedef _Scalar Scalar;
0039     typedef Matrix<Scalar,Dynamic,1> Vector;
0040   public:
0041     typedef typename Vector::StorageIndex StorageIndex;
0042     enum {
0043       ColsAtCompileTime = Dynamic,
0044       MaxColsAtCompileTime = Dynamic
0045     };
0046
0047     DiagonalPreconditioner() : m_isInitialized(false) {}
0048
0049     template<typename MatType>
0050     explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
0051     {
0052       compute(mat);
0053     }
0054
0055     EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_invdiag.size(); }
0056     EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_invdiag.size(); }
0057
0058     template<typename MatType>
0059     DiagonalPreconditioner& analyzePattern(const MatType& )
0060     {
0061       return *this;
0062     }
0063
0064     template<typename MatType>
0065     DiagonalPreconditioner& factorize(const MatType& mat)
0066     {
0067       m_invdiag.resize(mat.cols());
0068       for(int j=0; j<mat.outerSize(); ++j)
0069       {
0070         typename MatType::InnerIterator it(mat,j);
0071         while(it && it.index()!=j) ++it;
0072         if(it && it.index()==j && it.value()!=Scalar(0))
0073           m_invdiag(j) = Scalar(1)/it.value();
0074         else
0075           m_invdiag(j) = Scalar(1);
0076       }
0077       m_isInitialized = true;
0078       return *this;
0079     }
0080
0081     template<typename MatType>
0082     DiagonalPreconditioner& compute(const MatType& mat)
0083     {
0084       return factorize(mat);
0085     }
0086
0087     /** \internal */
0088     template<typename Rhs, typename Dest>
0089     void _solve_impl(const Rhs& b, Dest& x) const
0090     {
0091       x = m_invdiag.array() * b.array() ;
0092     }
0093
0094     template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
0095     solve(const MatrixBase<Rhs>& b) const
0096     {
0097       eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
0098       eigen_assert(m_invdiag.size()==b.rows()
0099                 && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
0100       return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
0101     }
0102
0103     ComputationInfo info() { return Success; }
0104
0105   protected:
0106     Vector m_invdiag;
0107     bool m_isInitialized;
0108 };
0109
0110 /** \ingroup IterativeLinearSolvers_Module
0111   * \brief Jacobi preconditioner for LeastSquaresConjugateGradient
0112   *
0113   * This class allows to approximately solve for A' A x  = A' b problems assuming A' A is a diagonal matrix.
0114   * In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:
0115     \code
0116     (A.adjoint() * A).diagonal().asDiagonal() * x = b
0117     \endcode
0118   *
0119   * \tparam _Scalar the type of the scalar.
0120   *
0121   * \implsparsesolverconcept
0122   *
0123   * The diagonal entries are pre-inverted and stored into a dense vector.
0124   *
0125   * \sa class LeastSquaresConjugateGradient, class DiagonalPreconditioner
0126   */
0127 template <typename _Scalar>
0128 class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar>
0129 {
0130     typedef _Scalar Scalar;
0131     typedef typename NumTraits<Scalar>::Real RealScalar;
0132     typedef DiagonalPreconditioner<_Scalar> Base;
0133     using Base::m_invdiag;
0134   public:
0135
0136     LeastSquareDiagonalPreconditioner() : Base() {}
0137
0138     template<typename MatType>
0139     explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
0140     {
0141       compute(mat);
0142     }
0143
0144     template<typename MatType>
0145     LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& )
0146     {
0147       return *this;
0148     }
0149
0150     template<typename MatType>
0151     LeastSquareDiagonalPreconditioner& factorize(const MatType& mat)
0152     {
0153       // Compute the inverse squared-norm of each column of mat
0154       m_invdiag.resize(mat.cols());
0155       if(MatType::IsRowMajor)
0156       {
0157         m_invdiag.setZero();
0158         for(Index j=0; j<mat.outerSize(); ++j)
0159         {
0160           for(typename MatType::InnerIterator it(mat,j); it; ++it)
0161             m_invdiag(it.index()) += numext::abs2(it.value());
0162         }
0163         for(Index j=0; j<mat.cols(); ++j)
0164           if(numext::real(m_invdiag(j))>RealScalar(0))
0165             m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j));
0166       }
0167       else
0168       {
0169         for(Index j=0; j<mat.outerSize(); ++j)
0170         {
0171           RealScalar sum = mat.col(j).squaredNorm();
0172           if(sum>RealScalar(0))
0173             m_invdiag(j) = RealScalar(1)/sum;
0174           else
0175             m_invdiag(j) = RealScalar(1);
0176         }
0177       }
0178       Base::m_isInitialized = true;
0179       return *this;
0180     }
0181
0182     template<typename MatType>
0183     LeastSquareDiagonalPreconditioner& compute(const MatType& mat)
0184     {
0185       return factorize(mat);
0186     }
0187
0188     ComputationInfo info() { return Success; }
0189
0190   protected:
0191 };
0192
0193 /** \ingroup IterativeLinearSolvers_Module
0194   * \brief A naive preconditioner which approximates any matrix as the identity matrix
0195   *
0196   * \implsparsesolverconcept
0197   *
0198   * \sa class DiagonalPreconditioner
0199   */
0200 class IdentityPreconditioner
0201 {
0202   public:
0203
0204     IdentityPreconditioner() {}
0205
0206     template<typename MatrixType>
0207     explicit IdentityPreconditioner(const MatrixType& ) {}
0208
0209     template<typename MatrixType>
0210     IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
0211
0212     template<typename MatrixType>
0213     IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
0214
0215     template<typename MatrixType>
0216     IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
0217
0218     template<typename Rhs>
0219     inline const Rhs& solve(const Rhs& b) const { return b; }
0220
0221     ComputationInfo info() { return Success; }
0222 };
0223
0224 } // end namespace RivetEigen
0225
0226 #endif // EIGEN_BASIC_PRECONDITIONERS_H