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 <string>
 
#include <dune/common/unused.hh>
 
#include "preconditioner.hh"
#include "solver.hh"
#include "solvercategory.hh"
#include "istlexception.hh"
#include "matrixutils.hh"
#include "gsetc.hh"
#include "ilu.hh"
 
 
namespace Dune {
  template<class O, int c>
  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 O InverseOperator;
 
    // define the category
    enum {
      category=c
    };
 
    InverseOperator2Preconditioner(InverseOperator& inverse_operator)
    : inverse_operator_(inverse_operator)
    {}
 
    void pre(domain_type&,range_type&)
    {}
 
    void apply(domain_type& v, const range_type& d)
    {
      InverseOperatorResult res;
      range_type copy(d);
      inverse_operator_.apply(v, copy, res);
    }
 
    void post(domain_type&)
    {}
 
  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;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqSSOR (const M& A, int n, field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(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 (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    const M& _A_;
    int _n;
    field_type _w;
  };
 
 
 
  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;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqSOR (const M& A, int n, field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(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 (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    const M& _A_;
    int _n;
    field_type _w;
  };
 
 
  template<class M, class X, class Y, int l=1>
  class SeqGS : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqGS (const M& A, int n, field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(b);
    }
 
    virtual void apply (X& v, const Y& d)
    {
      for (int i=0; i<_n; i++) {
        dbgs(_A_,v,d,_w,BL<l>());
      }
    }
 
    virtual void post (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    const M& _A_;
    int _n;
    field_type _w;
  };
 
 
  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;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqJac (const M& A, int n, field_type w)
      : _A_(A), _n(n), _w(w)
    {
      CheckIfDiagonalPresent<M,l>::check(_A_);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(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 (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    const M& _A_;
    int _n;
    field_type _w;
  };
 
 
 
  template<class M, class X, class Y, int l=1>
  class SeqILU0 : public Preconditioner<X,Y> {
  public:
    typedef typename Dune::remove_const<M>::type matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqILU0 (const M& A, field_type w)
      : ILU(A) // copy A
    {
      _w =w;
      bilu0_decomposition(ILU);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(b);
    }
 
    virtual void apply (X& v, const Y& d)
    {
      bilu_backsolve(ILU,v,d);
      v *= _w;
    }
 
    virtual void post (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    field_type _w;
    matrix_type ILU;
  };
 
 
  template<class M, class X, class Y, int l=1>
  class SeqILUn : public Preconditioner<X,Y> {
  public:
    typedef typename Dune::remove_const<M>::type matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    SeqILUn (const M& A, int n, field_type w)
      : ILU(A.N(),A.M(),M::row_wise)
    {
      _n = n;
      _w = w;
      bilu_decomposition(A,n,ILU);
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(b);
    }
 
    virtual void apply (X& v, const Y& d)
    {
      bilu_backsolve(ILU,v,d);
      v *= _w;
    }
 
    virtual void post (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    matrix_type ILU;
    int _n;
    field_type _w;
  };
 
 
 
  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;
 
    // define the category
    enum {
      category=SolverCategory::sequential
    };
 
    Richardson (field_type w=1.0)
    {
      _w = w;
    }
 
    virtual void pre (X& x, Y& b)
    {
      DUNE_UNUSED_PARAMETER(x);
      DUNE_UNUSED_PARAMETER(b);
    }
 
    virtual void apply (X& v, const Y& d)
    {
      v = d;
      v *= _w;
    }
 
    virtual void post (X& x)
    {
      DUNE_UNUSED_PARAMETER(x);
    }
 
  private:
    field_type _w;
  };
 
 
 
} // end namespace
 
#endif