preconditioners.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- 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 "dilu.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 typename FieldTraits<scalar_field_type>::real_type real_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");
    }

    void pre([[maybe_unused]] domain_type&,[[maybe_unused]] range_type&) override
    {}

    void apply(domain_type& v, const range_type& d) override
    {
      InverseOperatorResult res;
      range_type copy(d);
      inverse_operator_.apply(v, copy, res);
    }

    void post([[maybe_unused]] domain_type&) override
    {}

    SolverCategory::Category category() const override
    {
      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;
    typedef typename FieldTraits<scalar_field_type>::real_type real_field_type;

    SeqSSOR (const M& A, int n, real_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<real_field_type>("relaxation",1.0))
    {}

    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    void apply (X& v, const Y& d) override
    {
      for (int i=0; i<_n; i++) {
        bsorf(_A_,v,d,_w,BL<l>());
        bsorb(_A_,v,d,_w,BL<l>());
      }
    }

    void post ([[maybe_unused]] X& x) override
    {}

    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

  private:
    const M& _A_;
    int _n;
    real_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;
    typedef typename FieldTraits<scalar_field_type>::real_type real_field_type;

    SeqSOR (const M& A, int n, real_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<real_field_type>("relaxation",1.0))
    {}

    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    void apply (X& v, const Y& d) override
    {
      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>());
        }
    }

    void post ([[maybe_unused]] X& x) override
    {}

    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

  private:
    const M& _A_;
    int _n;
    real_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;
    typedef typename FieldTraits<scalar_field_type>::real_type real_field_type;

    SeqJac (const M& A, int n, real_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<real_field_type>("relaxation",1.0))
    {}

    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    void apply (X& v, const Y& d) override
    {
      for (int i=0; i<_n; i++) {
        dbjac(_A_,v,d,_w,BL<l>());
      }
    }

    void post ([[maybe_unused]] X& x) override
    {}

    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

  private:
    const M& _A_;
    int _n;
    real_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("jac", defaultPreconditionerBlockLevelCreator<Dune::SeqJac>());

  template <class M, class X, class Y, int l = 1>
  class SeqDILU : public Preconditioner<X, Y>
  {
  public:
    using matrix_type = M;
    using block_type = typename matrix_type::block_type;
    using domain_type = X;
    using range_type = Y;

    using field_type = typename X::field_type;

    using scalar_field_type = Simd::Scalar<field_type>;
    using real_field_type = typename FieldTraits<scalar_field_type>::real_type;

    SeqDILU(const M &A, real_field_type w)
        : _A_(A),
          _w(w),
          wNotIdentity_([w]
                        {using std::abs; return abs(w - real_field_type(1)) > 1e-15; }())
    {
      Dinv_.resize(_A_.N());
      CheckIfDiagonalPresent<M, l>::check(_A_);
      DILU::blockDILUDecomposition(_A_, Dinv_);
    }

    SeqDILU(const std::shared_ptr<const AssembledLinearOperator<M, X, Y>> &A, const ParameterTree &configuration)
        : SeqDILU(A->getmat(), configuration)
    {
    }

    SeqDILU(const M &A, const ParameterTree &config)
        : SeqDILU(A, config.get<real_field_type>("relaxation", 1.0))
    {
    }

    void pre([[maybe_unused]] X &x, [[maybe_unused]] Y &b) override
    {
    }

    void apply(X &v, const Y &d) override
    {

      DILU::blockDILUBacksolve(_A_, Dinv_, v, d);

      if (wNotIdentity_)
      {
        v *= _w;
      }
    }

    void post([[maybe_unused]] X &x) override
    {
    }

    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

  protected:
    std::vector<block_type> Dinv_;
    const M &_A_;
    const real_field_type _w;
    const bool wNotIdentity_;
  };
  DUNE_REGISTER_PRECONDITIONER("dilu", defaultPreconditionerBlockLevelCreator<Dune::SeqDILU>());

  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 FieldTraits<scalar_field_type>::real_type real_field_type;

    typedef typename ILU::CRS< block_type , typename M::allocator_type> CRS;

    SeqILU (const M& A, real_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<real_field_type>("relaxation", 1.0),
               config.get("resort", false))
    {}

    SeqILU (const M& A, int n, real_field_type w, const bool resort = false )
      : ILU_(),
        lower_(),
        upper_(),
        inv_(),
        w_(w),
        wNotIdentity_([w]{using std::abs; return abs(w - real_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();
      }
    }

    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    void apply (X& v, const Y& d) override
    {
      if( ILU_ )
      {
        ILU::blockILUBacksolve( *ILU_, v, d);
      }
      else
      {
        ILU::blockILUBacksolve(lower_, upper_, inv_, v, d);
      }

      if( wNotIdentity_ )
      {
        v *= w_;
      }
    }

    void post ([[maybe_unused]] X& x) override
    {}

    SolverCategory::Category category() const override
    {
      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 real_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;
    typedef typename FieldTraits<scalar_field_type>::real_type real_field_type;

    Richardson (real_field_type w=1.0) :
      _w(w)
    {}

    Richardson (const ParameterTree& configuration)
      : Richardson(configuration.get<real_field_type>("relaxation", 1.0))
    {}

    void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& b) override
    {}

    void apply (X& v, const Y& d) override
    {
      v = d;
      v *= _w;
    }

    void post ([[maybe_unused]] X& x) override
    {}

    SolverCategory::Category category() const override
    {
      return SolverCategory::sequential;
    }

  private:
    real_field_type _w;
  };
  DUNE_REGISTER_PRECONDITIONER("richardson", [](auto opTraits, const auto& /*op*/, const ParameterTree& config){
    using OpTraits = std::decay_t<decltype(opTraits)>;
    using D = typename OpTraits::domain_type;
    using R = typename OpTraits::range_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;
    typedef typename FieldTraits<scalar_field_type>::real_type real_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<real_field_type>("relaxation", 1.0))
    {}

    explicit SeqILDL ( const matrix_type &A, real_field_type relax = real_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_;
    real_field_type relax_;
  };
  DUNE_REGISTER_PRECONDITIONER("ildl", defaultPreconditionerCreator<Dune::SeqILDL>());

} // end namespace


#endif