scalarproducts.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_SCALARPRODUCTS_HH
#define DUNE_ISTL_SCALARPRODUCTS_HH
 
#include <cmath>
#include <complex>
#include <iostream>
#include <iomanip>
#include <string>
 
#include <dune/common/exceptions.hh>
 
#include "bvector.hh"
#include "solvercategory.hh"
 
 
namespace Dune {
  template<class X>
  class ScalarProduct {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
 
    virtual field_type dot (const X& x, const X& y) = 0;
 
    virtual real_type norm (const X& x) = 0;
 
    virtual SolverCategory::Category category() const
#if DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE
    {
      DUNE_THROW(Dune::Exception,"It is necessary to implement the category method in a derived classes, in the future this method will pure virtual.");
    }
#else
    = 0;
#endif
 
    virtual ~ScalarProduct () {}
  };
 
  //=====================================================================
  // Implementation for ISTL-matrix based operator
  //=====================================================================
 
  template<class X>
  class SeqScalarProduct : public ScalarProduct<X>
  {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
 
    virtual field_type dot (const X& x, const X& y)
    {
      return x.dot(y);
    }
 
    virtual real_type norm (const X& x)
    {
      return x.two_norm();
    }
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::sequential;
    }
 
  };
 
  template<class X, class C>
  class NonoverlappingSchwarzScalarProduct : public ScalarProduct<X>
  {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
    typedef C communication_type;
 
    NonoverlappingSchwarzScalarProduct (const communication_type& com)
      : communication(com)
    {}
 
    virtual field_type dot (const X& x, const X& y)
    {
      field_type result(0);
      communication.dot(x,y,result);
      return result;
    }
 
    virtual real_type norm (const X& x)
    {
      return communication.norm(x);
    }
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::nonoverlapping;
    }
 
    void make_consistent (X& x) const
    {
      communication.copyOwnerToAll(x,x);
    }
 
  private:
    const communication_type& communication;
  };
 
  template<class X, class C>
  class OverlappingSchwarzScalarProduct : public ScalarProduct<X>
  {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
    typedef C communication_type;
 
    OverlappingSchwarzScalarProduct (const communication_type& com)
      : communication(com)
    {}
 
    virtual field_type dot (const X& x, const X& y)
    {
      field_type result(0);
      communication.dot(x,y,result);
      return result;
    }
 
    virtual real_type norm (const X& x)
    {
      return communication.norm(x);
    }
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::overlapping;
    }
 
  private:
    const communication_type& communication;
  };
 
  template<class X, class Comm>
  std::shared_ptr<ScalarProduct<X>> createScalarProduct(const Comm& comm, SolverCategory::Category category)
  {
    switch(category)
    {
      case SolverCategory::sequential:
        return
          std::make_shared<SeqScalarProduct<X>>();
      case SolverCategory::nonoverlapping:
        return
          std::make_shared<NonoverlappingSchwarzScalarProduct<X,Comm>>(comm);
      case SolverCategory::overlapping:
        return
          std::make_shared<OverlappingSchwarzScalarProduct<X,Comm>>(comm);
      default:
        DUNE_THROW(InvalidStateException, "unknown solver category");
    }
  }
 
} // end namespace Dune
 
#endif