BasicPreconditioners.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
#ifndef EIGEN_BASIC_PRECONDITIONERS_H
#define EIGEN_BASIC_PRECONDITIONERS_H
 
namespace Eigen { 
 
template <typename _Scalar>
class DiagonalPreconditioner
{
    typedef _Scalar Scalar;
    typedef Matrix<Scalar,Dynamic,1> Vector;
  public:
    typedef typename Vector::StorageIndex StorageIndex;
    enum {
      ColsAtCompileTime = Dynamic,
      MaxColsAtCompileTime = Dynamic
    };
 
    DiagonalPreconditioner() : m_isInitialized(false) {}
 
    template<typename MatType>
    explicit DiagonalPreconditioner(const MatType& mat) : m_invdiag(mat.cols())
    {
      compute(mat);
    }
 
    Index rows() const { return m_invdiag.size(); }
    Index cols() const { return m_invdiag.size(); }
    
    template<typename MatType>
    DiagonalPreconditioner& analyzePattern(const MatType& )
    {
      return *this;
    }
    
    template<typename MatType>
    DiagonalPreconditioner& factorize(const MatType& mat)
    {
      m_invdiag.resize(mat.cols());
      for(int j=0; j<mat.outerSize(); ++j)
      {
        typename MatType::InnerIterator it(mat,j);
        while(it && it.index()!=j) ++it;
        if(it && it.index()==j && it.value()!=Scalar(0))
          m_invdiag(j) = Scalar(1)/it.value();
        else
          m_invdiag(j) = Scalar(1);
      }
      m_isInitialized = true;
      return *this;
    }
    
    template<typename MatType>
    DiagonalPreconditioner& compute(const MatType& mat)
    {
      return factorize(mat);
    }
 
    template<typename Rhs, typename Dest>
    void _solve_impl(const Rhs& b, Dest& x) const
    {
      x = m_invdiag.array() * b.array() ;
    }
 
    template<typename Rhs> inline const Solve<DiagonalPreconditioner, Rhs>
    solve(const MatrixBase<Rhs>& b) const
    {
      eigen_assert(m_isInitialized && "DiagonalPreconditioner is not initialized.");
      eigen_assert(m_invdiag.size()==b.rows()
                && "DiagonalPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
      return Solve<DiagonalPreconditioner, Rhs>(*this, b.derived());
    }
    
    ComputationInfo info() { return Success; }
 
  protected:
    Vector m_invdiag;
    bool m_isInitialized;
};
 
template <typename _Scalar>
class LeastSquareDiagonalPreconditioner : public DiagonalPreconditioner<_Scalar>
{
    typedef _Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef DiagonalPreconditioner<_Scalar> Base;
    using Base::m_invdiag;
  public:
 
    LeastSquareDiagonalPreconditioner() : Base() {}
 
    template<typename MatType>
    explicit LeastSquareDiagonalPreconditioner(const MatType& mat) : Base()
    {
      compute(mat);
    }
 
    template<typename MatType>
    LeastSquareDiagonalPreconditioner& analyzePattern(const MatType& )
    {
      return *this;
    }
    
    template<typename MatType>
    LeastSquareDiagonalPreconditioner& factorize(const MatType& mat)
    {
      // Compute the inverse squared-norm of each column of mat
      m_invdiag.resize(mat.cols());
      if(MatType::IsRowMajor)
      {
        m_invdiag.setZero();
        for(Index j=0; j<mat.outerSize(); ++j)
        {
          for(typename MatType::InnerIterator it(mat,j); it; ++it)
            m_invdiag(it.index()) += numext::abs2(it.value());
        }
        for(Index j=0; j<mat.cols(); ++j)
          if(numext::real(m_invdiag(j))>RealScalar(0))
            m_invdiag(j) = RealScalar(1)/numext::real(m_invdiag(j));
      }
      else
      {
        for(Index j=0; j<mat.outerSize(); ++j)
        {
          RealScalar sum = mat.col(j).squaredNorm();
          if(sum>RealScalar(0))
            m_invdiag(j) = RealScalar(1)/sum;
          else
            m_invdiag(j) = RealScalar(1);
        }
      }
      Base::m_isInitialized = true;
      return *this;
    }
    
    template<typename MatType>
    LeastSquareDiagonalPreconditioner& compute(const MatType& mat)
    {
      return factorize(mat);
    }
    
    ComputationInfo info() { return Success; }
 
  protected:
};
 
class IdentityPreconditioner
{
  public:
 
    IdentityPreconditioner() {}
 
    template<typename MatrixType>
    explicit IdentityPreconditioner(const MatrixType& ) {}
    
    template<typename MatrixType>
    IdentityPreconditioner& analyzePattern(const MatrixType& ) { return *this; }
    
    template<typename MatrixType>
    IdentityPreconditioner& factorize(const MatrixType& ) { return *this; }
 
    template<typename MatrixType>
    IdentityPreconditioner& compute(const MatrixType& ) { return *this; }
    
    template<typename Rhs>
    inline const Rhs& solve(const Rhs& b) const { return b; }
    
    ComputationInfo info() { return Success; }
};
 
} // end namespace Eigen
 
#endif // EIGEN_BASIC_PRECONDITIONERS_H