symmetricOperators.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2013-2019 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef SYMMETRICOPERATORS_HH
#define SYMMETRICOPERATORS_HH
 
#include "dune/istl/preconditioners.hh"
 
namespace Kaskade
{
  template <class X, class Xstar>
  class DualPairing {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    
    virtual field_type operator()(X const& x, Xstar const& y) const = 0;
  };
  
  template <class X, class Xstar>
  class DefaultDualPairing: public DualPairing<X,Xstar> {
    typedef typename DualPairing<X,Xstar>::field_type field_type;
  public:
    virtual field_type operator()(X const& x, Xstar const& y) const { return x*y; }
  };
  
  template <class X, class Xstar>
  class ZeroDualPairing: public DualPairing<X,Xstar> {
    typedef typename DualPairing<X,Xstar>::field_type field_type;
  public:
    virtual field_type operator()(X const& , Xstar const& ) const { return 0; }
  };
  
  //--------------------------------------------------------------------------------------------------
  
  template <class X, class Xstar>
  class SymmetricLinearOperator: public Dune::LinearOperator<X,Xstar>
  {
  public:
    typedef typename X::field_type field_type;
    
    
    virtual field_type applyDp(X const& x, Xstar& y) const
    {
      
      this->apply(x,y);
      return dp(x,y);
    }
    
    virtual field_type dp(X const& x, Xstar const& y) const = 0;
 
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
  };
  
  template <class X, class Matrix>
  class SymmetricMatrixOperator: public SymmetricLinearOperator<X,X>
  {
  public:
    typedef typename X::field_type field_type;
    
    SymmetricMatrixOperator(Matrix const& A_)
    : A(A_)
    {}
    
    virtual void apply (X const& x, X& y) const 
    {
      A.mv(x,y);
    }
    
    virtual void applyscaleadd (field_type alpha, X const& x, X& y) const
    {
      A.usmv(alpha,x,y);
    }
    
    virtual field_type applyDp(X const& x, X& y) const
    {
      A.mv(x,y);
      return dp(x,y);
    }
    
    virtual field_type dp(X const& x, X const& y) const
    {
      return x*y;
    }
    
  private:
    Matrix const& A;
  };
  
  template <class X, class Xstar>
  class SymmetricLinearOperatorWrapper: public SymmetricLinearOperator<X,Xstar>
  {
  public:
    typedef typename X::field_type field_type;
    
    SymmetricLinearOperatorWrapper(Dune::LinearOperator<X,Xstar> const& A_, DualPairing<X,Xstar> const& dualPairing_)
    : A(A_), dualPairing(dualPairing_) {}
    
    virtual void apply(X const& x, Xstar& y) const
    {
      A.apply(x,y);
    }
    
    virtual void applyscaleadd(field_type alpha, X const& x, Xstar& y) const
    {
      A.applyscaleadd(alpha,x,y);
    }
    
    virtual field_type dp(X const& x, Xstar const& y) const
    {
      return dualPairing(x,y);
    }
    
  private:
    Dune::LinearOperator<X,Xstar> const& A;
    DualPairing<X,Xstar> const& dualPairing;
  };
    
  //--------------------------------------------------------------------------------------------------
 
  template <class X, class Xstar>
  class SymmetricPreconditioner: public Dune::Preconditioner<X,Xstar> 
  {
  public:
    typedef typename X::field_type field_type;
    
    virtual void pre(X& , Xstar& ) {}
    
    virtual void post(X& x) {}
    
    virtual field_type applyDp(X& x, Xstar const& y) = 0;
    
    virtual bool requiresInitializedInput() const = 0;
 
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
  };
  
  template <class X>
  class IdentityPreconditioner: public SymmetricPreconditioner<X,X> 
  {
  public:
    typedef typename X::field_type field_type;
    
    virtual void apply (X& v, const X& d) 
    {
      v = d;
    }
    
    virtual field_type applyDp(X& x, X const& y)
    {
      x = y;
      return x.two_norm2();
    }
    
    virtual bool requiresInitializedInput() const
    {
      return false;
    }
  };
  
  template <class X, class Xstar>
  class SymmetricPreconditionerWrapper: public SymmetricPreconditioner<X,Xstar>
  {
  public:
    typedef typename X::field_type field_type;
    
    SymmetricPreconditionerWrapper(Dune::Preconditioner<X,Xstar>& B_, DualPairing<X,Xstar> const& dualPairing_)
    : B(B_), dualPairing(dualPairing_) {}
    
    virtual void pre(X& x, Xstar& y)
    {
      B.pre(x,y);
    }
    
    virtual void post(X& x)
    {
      B.post(x);
    }
    
    virtual void apply(X& x, Xstar const& y)
    {
      B.apply(x,y);
    }
    
    virtual field_type applyDp(X& x, Xstar const& y)
    {
      B.apply(x,y); 
      return dualPairing(x,y);
    }
    
    virtual bool requiresInitializedInput() const { return true; }
 
  private:
    Dune::Preconditioner<X,Xstar>& B;
    DualPairing<X,Xstar> const& dualPairing;
  };
  
} // end of namespace Kaskade
 
#endif