preconditioners.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_PRECONDITIONERS_HH
#define DUNE_ISTL_PRECONDITIONERS_HH
 
#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip>
#include <memory>
#include <string>
 
#include <dune/common/simd/simd.hh>
#include <dune/common/parametertree.hh>
 
#include <dune/istl/solverregistry.hh>
#include "preconditioner.hh"
#include "solver.hh"
#include "solvercategory.hh"
#include "istlexception.hh"
#include "matrixutils.hh"
#include "gsetc.hh"
#include "ildl.hh"
#include "ilu.hh"
 
 
namespace Dune {
  template<class O, int c = -1>
  class InverseOperator2Preconditioner :
    public Preconditioner<typename O::domain_type, typename O::range_type>
  {
  public:
    typedef typename O::domain_type domain_type;
    typedef typename O::range_type range_type;
    typedef typename range_type::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
    typedef O InverseOperator;
 
    InverseOperator2Preconditioner(InverseOperator& inverse_operator)
    : inverse_operator_(inverse_operator)
    {
      if(c != -1 && SolverCategory::category(inverse_operator_) != c)
        DUNE_THROW(InvalidStateException, "User-supplied solver category does not match that of the given inverse operator");
    }
 
    virtual void pre(domain_type&,range_type&)
    {}
 
    virtual void apply(domain_type& v, const range_type& d)
    {
      InverseOperatorResult res;
      range_type copy(d);
      inverse_operator_.apply(v, copy, res);
    }
 
    virtual void post(domain_type&)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::category(inverse_operator_);
    }
 
  private:
    InverseOperator& inverse_operator_;
  };
 
  //=====================================================================
  // Implementation of this interface for sequential ISTL-preconditioners
  //=====================================================================
 
 
  template<class M, class X, class Y, int l=1>
  class SeqSSOR : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    SeqSSOR (const M& A, int n, scalar_field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    SeqSSOR (const std::shared_ptr<const AssembledLinearOperator<M,X,Y>>& A, const ParameterTree& configuration)
      : SeqSSOR(A->getmat(), configuration)
    {}
 
    SeqSSOR (const M& A, const ParameterTree& configuration)
      : SeqSSOR(A, configuration.get<int>("iterations",1), configuration.get<scalar_field_type>("relaxation",1.0))
    {}
 
    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
    {}
 
    virtual void apply (X& v, const Y& d)
    {
      for (int i=0; i<_n; i++) {
        bsorf(_A_,v,d,_w,BL<l>());
        bsorb(_A_,v,d,_w,BL<l>());
      }
    }
 
    virtual void post ([[maybe_unused]] X& x)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  private:
    const M& _A_;
    int _n;
    scalar_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("ssor", defaultPreconditionerBlockLevelCreator<Dune::SeqSSOR>());
 
 
  template<class M, class X, class Y, int l=1>
  class SeqSOR : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    SeqSOR (const M& A, int n, scalar_field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    SeqSOR (const std::shared_ptr<const AssembledLinearOperator<M,X,Y>>& A, const ParameterTree& configuration)
      : SeqSOR(A->getmat(), configuration)
    {}
 
    SeqSOR (const M& A, const ParameterTree& configuration)
      : SeqSOR(A, configuration.get<int>("iterations",1), configuration.get<scalar_field_type>("relaxation",1.0))
    {}
 
    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
    {}
 
    virtual void apply (X& v, const Y& d)
    {
      this->template apply<true>(v,d);
    }
 
    template<bool forward>
    void apply(X& v, const Y& d)
    {
      if(forward)
        for (int i=0; i<_n; i++) {
          bsorf(_A_,v,d,_w,BL<l>());
        }
      else
        for (int i=0; i<_n; i++) {
          bsorb(_A_,v,d,_w,BL<l>());
        }
    }
 
    virtual void post ([[maybe_unused]] X& x)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  private:
    const M& _A_;
    int _n;
    scalar_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("sor", defaultPreconditionerBlockLevelCreator<Dune::SeqSOR>());
 
 
  template<class M, class X, class Y, int l=1>
  using SeqGS = SeqSOR<M,X,Y,l>;
  DUNE_REGISTER_PRECONDITIONER("gs", defaultPreconditionerBlockLevelCreator<Dune::SeqGS>());
 
  template<class M, class X, class Y, int l=1>
  class SeqJac : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    SeqJac (const M& A, int n, scalar_field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    SeqJac (const std::shared_ptr<const AssembledLinearOperator<M,X,Y>>& A, const ParameterTree& configuration)
      : SeqJac(A->getmat(), configuration)
    {}
 
    SeqJac (const M& A, const ParameterTree& configuration)
      : SeqJac(A, configuration.get<int>("iterations",1), configuration.get<scalar_field_type>("relaxation",1.0))
    {}
 
    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
    {}
 
    virtual void apply (X& v, const Y& d)
    {
      for (int i=0; i<_n; i++) {
        dbjac(_A_,v,d,_w,BL<l>());
      }
    }
 
    virtual void post ([[maybe_unused]] X& x)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  private:
    const M& _A_;
    int _n;
    scalar_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("jac", defaultPreconditionerBlockLevelCreator<Dune::SeqJac>());
 
 
 
  template<class M, class X, class Y, int l=1>
  class SeqILU : public Preconditioner<X,Y> {
  public:
    typedef typename std::remove_const<M>::type matrix_type;
    typedef typename matrix_type :: block_type block_type;
    typedef X domain_type;
    typedef Y range_type;
 
    typedef typename X::field_type field_type;
 
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    typedef typename ILU::CRS< block_type , typename M::allocator_type> CRS;
 
    SeqILU (const M& A, scalar_field_type w, const bool resort = false )
      : SeqILU( A, 0, w, resort ) // construct ILU(0)
    {
    }
 
    SeqILU (const std::shared_ptr<const AssembledLinearOperator<M,X,Y>>& A, const ParameterTree& configuration)
      : SeqILU(A->getmat(), configuration)
    {}
 
    SeqILU(const M& A, const ParameterTree& config)
      : SeqILU(A, config.get("n", 0),
               config.get<scalar_field_type>("relaxation", 1.0),
               config.get("resort", false))
    {}
 
    SeqILU (const M& A, int n, scalar_field_type w, const bool resort = false )
      : ILU_(),
        lower_(),
        upper_(),
        inv_(),
        w_(w),
        wNotIdentity_([w]{using std::abs; return abs(w - scalar_field_type(1)) > 1e-15;}() )
    {
      if( n == 0 )
      {
        // copy A
        ILU_.reset( new matrix_type( A ) );
        // create ILU(0) decomposition
        ILU::blockILU0Decomposition( *ILU_ );
      }
      else
      {
        // create matrix in build mode
        ILU_.reset( new matrix_type(  A.N(), A.M(), matrix_type::row_wise) );
        // create ILU(n) decomposition
        ILU::blockILUDecomposition( A, n, *ILU_ );
      }
 
      if( resort )
      {
        // store ILU in simple CRS format
        ILU::convertToCRS( *ILU_, lower_, upper_, inv_ );
        ILU_.reset();
      }
    }
 
    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
    {}
 
    virtual void apply (X& v, const Y& d)
    {
      if( ILU_ )
      {
        ILU::blockILUBacksolve( *ILU_, v, d);
      }
      else
      {
        ILU::blockILUBacksolve(lower_, upper_, inv_, v, d);
      }
 
      if( wNotIdentity_ )
      {
        v *= w_;
      }
    }
 
    virtual void post ([[maybe_unused]] X& x)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  protected:
    std::unique_ptr< matrix_type > ILU_;
 
    CRS lower_;
    CRS upper_;
    std::vector< block_type, typename matrix_type::allocator_type > inv_;
 
    const scalar_field_type w_;
    const bool wNotIdentity_;
  };
  DUNE_REGISTER_PRECONDITIONER("ilu", defaultPreconditionerBlockLevelCreator<Dune::SeqILU>());
 
 
  template<class X, class Y>
  class Richardson : public Preconditioner<X,Y> {
  public:
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    Richardson (scalar_field_type w=1.0) :
      _w(w)
    {}
 
    Richardson (const ParameterTree& configuration)
      : Richardson(configuration.get<scalar_field_type>("relaxation", 1.0))
    {}
 
    virtual void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b)
    {}
 
    virtual void apply (X& v, const Y& d)
    {
      v = d;
      v *= _w;
    }
 
    virtual void post ([[maybe_unused]] X& x)
    {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  private:
    scalar_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("richardson", [](auto tl, const auto& /* mat */, const ParameterTree& config){
                                               using D = typename Dune::TypeListElement<1, decltype(tl)>::type;
                                               using R = typename Dune::TypeListElement<2, decltype(tl)>::type;
                                               return std::make_shared<Richardson<D,R>>(config);
                                             });
 
 
  template< class M, class X, class Y >
  class SeqILDL
    : public Preconditioner< X, Y >
  {
    typedef SeqILDL< M, X, Y > This;
    typedef Preconditioner< X, Y > Base;
 
  public:
    typedef std::remove_const_t< M > matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef Simd::Scalar<field_type> scalar_field_type;
 
    SeqILDL (const std::shared_ptr<const AssembledLinearOperator<M,X,Y>>& A, const ParameterTree& configuration)
      : SeqILDL(A->getmat(), configuration)
    {}
 
    SeqILDL(const matrix_type& A, const ParameterTree& config)
      : SeqILDL(A, config.get<scalar_field_type>("relaxation", 1.0))
    {}
 
    explicit SeqILDL ( const matrix_type &A, scalar_field_type relax = scalar_field_type( 1 ) )
      : decomposition_( A.N(), A.M(), matrix_type::random ),
        relax_( relax )
    {
      // setup row sizes for lower triangular matrix
      for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
      {
        const auto &A_i = *i;
        const auto ij = A_i.find( i.index() );
        if( ij != A_i.end() )
          decomposition_.setrowsize( i.index(), ij.offset()+1 );
        else
          DUNE_THROW( ISTLError, "diagonal entry missing" );
      }
      decomposition_.endrowsizes();
 
      // setup row indices for lower triangular matrix
      for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
      {
        const auto &A_i = *i;
        for( auto ij = A_i.begin(); ij.index() < i.index() ; ++ij )
          decomposition_.addindex( i.index(), ij.index() );
        decomposition_.addindex( i.index(), i.index() );
      }
      decomposition_.endindices();
 
      // copy values of lower triangular matrix
      auto i = A.begin();
      for( auto row = decomposition_.begin(), rowend = decomposition_.end(); row != rowend; ++row, ++i )
      {
        auto ij = i->begin();
        for( auto col = row->begin(), colend = row->end(); col != colend; ++col, ++ij )
          *col = *ij;
      }
 
      // perform ILDL decomposition
      bildl_decompose( decomposition_ );
    }
 
    void pre ([[maybe_unused]] X &x, [[maybe_unused]] Y &b) override
    {}
 
    void apply ( X &v, const Y &d ) override
    {
      bildl_backsolve( decomposition_, v, d, true );
      v *= relax_;
    }
 
    void post ([[maybe_unused]] X &x) override
    {}
 
    SolverCategory::Category category () const override { return SolverCategory::sequential; }
 
  private:
    matrix_type decomposition_;
    scalar_field_type relax_;
  };
  DUNE_REGISTER_PRECONDITIONER("ildl", defaultPreconditionerCreator<Dune::SeqILDL>());
 
} // end namespace
 
 
#endif