owneroverlapcopy.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_OWNEROVERLAPCOPY_HH
#define DUNE_ISTL_OWNEROVERLAPCOPY_HH
 
#include <new>
#include <iostream>
#include <vector>
#include <list>
#include <map>
#include <set>
#include <tuple>
 
#include "cmath"
 
// MPI header
#if HAVE_MPI
#include <mpi.h>
#endif
 
#include <dune/common/enumset.hh>
 
#if HAVE_MPI
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/communicator.hh>
#include <dune/common/parallel/remoteindices.hh>
#include <dune/common/parallel/mpicommunication.hh>
#endif
 
#include "solvercategory.hh"
#include "istlexception.hh"
#include <dune/common/parallel/communication.hh>
#include <dune/istl/matrixmarket.hh>
 
template<int dim, template<class,class> class Comm>
void testRedistributed(int s);
 
 
namespace Dune {
 
  struct OwnerOverlapCopyAttributeSet
  {
    enum AttributeSet {
      owner=1, overlap=2, copy=3
    };
  };
 
  template <class G, class L>
  class IndexInfoFromGrid
  {
  public:
    typedef G GlobalIdType;
 
    typedef L LocalIdType;
 
    typedef std::tuple<GlobalIdType,LocalIdType,int> IndexTripel;
    typedef std::tuple<int,GlobalIdType,int> RemoteIndexTripel;
 
    void addLocalIndex (const IndexTripel& x)
    {
      if (std::get<2>(x)!=OwnerOverlapCopyAttributeSet::owner &&
          std::get<2>(x)!=OwnerOverlapCopyAttributeSet::overlap &&
          std::get<2>(x)!=OwnerOverlapCopyAttributeSet::copy)
        DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
      localindices.insert(x);
    }
 
    void addRemoteIndex (const RemoteIndexTripel& x)
    {
      if (std::get<2>(x)!=OwnerOverlapCopyAttributeSet::owner &&
          std::get<2>(x)!=OwnerOverlapCopyAttributeSet::overlap &&
          std::get<2>(x)!=OwnerOverlapCopyAttributeSet::copy)
        DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
      remoteindices.insert(x);
    }
 
    const std::set<IndexTripel>& localIndices () const
    {
      return localindices;
    }
 
    const std::set<RemoteIndexTripel>& remoteIndices () const
    {
      return remoteindices;
    }
 
    void clear ()
    {
      localindices.clear();
      remoteindices.clear();
    }
 
  private:
    std::set<IndexTripel> localindices;
    std::set<RemoteIndexTripel> remoteindices;
  };
 
 
#if HAVE_MPI
 
  template <class GlobalIdType, class LocalIdType=int>
  class OwnerOverlapCopyCommunication
  {
    template<typename M, typename G, typename L>
    friend void loadMatrixMarket(M&,
                                 const std::string&,
                                 OwnerOverlapCopyCommunication<G,L>&,
                                 bool);
    // used types
    typedef typename IndexInfoFromGrid<GlobalIdType,LocalIdType>::IndexTripel IndexTripel;
    typedef typename IndexInfoFromGrid<GlobalIdType,LocalIdType>::RemoteIndexTripel RemoteIndexTripel;
    typedef typename std::set<IndexTripel>::const_iterator localindex_iterator;
    typedef typename std::set<RemoteIndexTripel>::const_iterator remoteindex_iterator;
    typedef typename OwnerOverlapCopyAttributeSet::AttributeSet AttributeSet;
    typedef Dune::ParallelLocalIndex<AttributeSet> LI;
  public:
    typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> PIS;
    typedef Dune::RemoteIndices<PIS> RI;
    typedef Dune::RemoteIndexListModifier<PIS,typename RI::Allocator,false> RILM;
    typedef typename RI::RemoteIndex RX;
    typedef Dune::BufferedCommunicator BC;
    typedef Dune::Interface IF;
    typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner> OwnerSet;
    typedef EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopySet;
    typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet;
    typedef Dune::AllSet<AttributeSet> AllSet;
  protected:
 
 
    template<typename T>
    struct CopyGatherScatter
    {
      typedef typename CommPolicy<T>::IndexedType V;
 
      static V gather(const T& a, std::size_t i)
      {
        return a[i];
      }
 
      static void scatter(T& a, V v, std::size_t i)
      {
        a[i] = v;
      }
    };
    template<typename T>
    struct AddGatherScatter
    {
      typedef typename CommPolicy<T>::IndexedType V;
 
      static V gather(const T& a, std::size_t i)
      {
        return a[i];
      }
 
      static void scatter(T& a, V v, std::size_t i)
      {
        a[i] += v;
      }
    };
 
    void buildOwnerOverlapToAllInterface () const
    {
      if (OwnerOverlapToAllInterfaceBuilt)
        OwnerOverlapToAllInterface.free();
      typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> OwnerOverlapSet;
      typedef Combine<OwnerOverlapSet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> AllSet;
      OwnerOverlapSet sourceFlags;
      AllSet destFlags;
      OwnerOverlapToAllInterface.build(ri,sourceFlags,destFlags);
      OwnerOverlapToAllInterfaceBuilt = true;
    }
 
    void buildOwnerToAllInterface () const
    {
      if (OwnerToAllInterfaceBuilt)
        OwnerToAllInterface.free();
      OwnerSet sourceFlags;
      AllSet destFlags;
      OwnerToAllInterface.build(ri,sourceFlags,destFlags);
      OwnerToAllInterfaceBuilt = true;
    }
 
    void buildOwnerCopyToAllInterface () const
    {
      if (OwnerCopyToAllInterfaceBuilt)
        OwnerCopyToAllInterface.free();
      typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet;
      typedef Combine<OwnerCopySet,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::overlap>,AttributeSet> AllSet;
      OwnerCopySet sourceFlags;
      AllSet destFlags;
      OwnerCopyToAllInterface.build(ri,sourceFlags,destFlags);
      OwnerCopyToAllInterfaceBuilt = true;
    }
 
    void buildOwnerCopyToOwnerCopyInterface () const
    {
      if (OwnerCopyToOwnerCopyInterfaceBuilt)
        OwnerCopyToOwnerCopyInterface.free();
      typedef Combine<EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::owner>,EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy>,AttributeSet> OwnerCopySet;
      OwnerCopySet sourceFlags;
      OwnerCopySet destFlags;
      OwnerCopyToOwnerCopyInterface.build(ri,sourceFlags,destFlags);
      OwnerCopyToOwnerCopyInterfaceBuilt = true;
    }
 
    void buildCopyToAllInterface () const
    {
      if (CopyToAllInterfaceBuilt)
        CopyToAllInterface.free();
      CopySet sourceFlags;
      AllSet destFlags;
      CopyToAllInterface.build(ri,sourceFlags,destFlags);
      CopyToAllInterfaceBuilt = true;
    }
 
  public:
 
    void
    DUNE_DEPRECATED_MSG("The solver category can only be set in the constructor. This method is deprecated and will be removed after Dune 2.7")
    setSolverCategory (SolverCategory::Category set) {
      category_ = set;
    }
 
    SolverCategory::Category
    DUNE_DEPRECATED_MSG("This method is deprecated and will be removed after Dune 2.7, use category() instead.")
    getSolverCategory () const {
      return category_;
    }
 
    SolverCategory::Category category () const {
      return category_;
    }
 
    const CollectiveCommunication<MPI_Comm>& communicator() const
    {
      return cc;
    }
 
    template<class T>
    void copyOwnerToAll (const T& source, T& dest) const
    {
      if (!OwnerToAllInterfaceBuilt)
        buildOwnerToAllInterface ();
      BC communicator;
      communicator.template build<T>(OwnerToAllInterface);
      communicator.template forward<CopyGatherScatter<T> >(source,dest);
      communicator.free();
    }
 
    template<class T>
    void copyCopyToAll (const T& source, T& dest) const
    {
      if (!CopyToAllInterfaceBuilt)
        buildCopyToAllInterface ();
      BC communicator;
      communicator.template build<T>(CopyToAllInterface);
      communicator.template forward<CopyGatherScatter<T> >(source,dest);
      communicator.free();
    }
 
    template<class T>
    void addOwnerOverlapToAll (const T& source, T& dest) const
    {
      if (!OwnerOverlapToAllInterfaceBuilt)
        buildOwnerOverlapToAllInterface ();
      BC communicator;
      communicator.template build<T>(OwnerOverlapToAllInterface);
      communicator.template forward<AddGatherScatter<T> >(source,dest);
      communicator.free();
    }
 
    template<class T>
    void addOwnerCopyToAll (const T& source, T& dest) const
    {
      if (!OwnerCopyToAllInterfaceBuilt)
        buildOwnerCopyToAllInterface ();
      BC communicator;
      communicator.template build<T>(OwnerCopyToAllInterface);
      communicator.template forward<AddGatherScatter<T> >(source,dest);
      communicator.free();
    }
 
    template<class T>
    void addOwnerCopyToOwnerCopy (const T& source, T& dest) const
    {
      if (!OwnerCopyToOwnerCopyInterfaceBuilt)
        buildOwnerCopyToOwnerCopyInterface ();
      BC communicator;
      communicator.template build<T>(OwnerCopyToOwnerCopyInterface);
      communicator.template forward<AddGatherScatter<T> >(source,dest);
      communicator.free();
    }
 
 
    template<class T1, class T2>
    void dot (const T1& x, const T1& y, T2& result) const
    {
      // set up mask vector
      if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
      {
        mask.resize(x.size());
        for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
          mask[i] = 1;
        for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
          if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
            mask[i->local().local()] = 0;
      }
      result = T2(0.0);
 
      for (typename T1::size_type i=0; i<x.size(); i++)
        result += x[i]*(y[i])*mask[i];
      result = cc.sum(result);
    }
 
    template<class T1>
    typename FieldTraits<typename T1::field_type>::real_type norm (const T1& x) const
    {
      using real_type = typename FieldTraits<typename T1::field_type>::real_type;
 
      // set up mask vector
      if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size()))
      {
        mask.resize(x.size());
        for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
          mask[i] = 1;
        for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
          if (i->local().attribute()!=OwnerOverlapCopyAttributeSet::owner)
            mask[i->local().local()] = 0;
      }
      auto result = real_type(0.0);
      for (typename T1::size_type i=0; i<x.size(); i++)
        result += Impl::asVector(x[i]).two_norm2()*mask[i];
      return sqrt(cc.sum(result));
    }
 
    typedef Dune::EnumItem<AttributeSet,OwnerOverlapCopyAttributeSet::copy> CopyFlags;
 
    typedef Dune::ParallelIndexSet<GlobalIdType,LI,512> ParallelIndexSet;
 
    typedef Dune::RemoteIndices<PIS> RemoteIndices;
 
    typedef Dune::GlobalLookupIndexSet<ParallelIndexSet> GlobalLookupIndexSet;
 
    const ParallelIndexSet& indexSet() const
    {
      return pis;
    }
 
    const RemoteIndices& remoteIndices() const
    {
      return ri;
    }
 
    ParallelIndexSet& indexSet()
    {
      return pis;
    }
 
 
    RemoteIndices& remoteIndices()
    {
      return ri;
    }
 
    void buildGlobalLookup()
    {
      if(globalLookup_) {
        if(pis.seqNo()==oldseqNo)
          // Nothing changed!
          return;
        delete globalLookup_;
      }
 
      globalLookup_ = new GlobalLookupIndexSet(pis);
      oldseqNo = pis.seqNo();
    }
 
    void buildGlobalLookup(std::size_t size)
    {
      if(globalLookup_) {
        if(pis.seqNo()==oldseqNo)
          // Nothing changed!
          return;
        delete globalLookup_;
      }
      globalLookup_ = new GlobalLookupIndexSet(pis, size);
      oldseqNo = pis.seqNo();
    }
 
    void freeGlobalLookup()
    {
      delete globalLookup_;
      globalLookup_=0;
    }
 
    const GlobalLookupIndexSet& globalLookup() const
    {
      assert(globalLookup_ != 0);
      return *globalLookup_;
    }
 
    template<class T1>
    void project (T1& x) const
    {
      for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
        if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
          x[i->local().local()] = 0;
    }
 
    OwnerOverlapCopyCommunication (MPI_Comm comm_,
                                   SolverCategory::Category cat_ = SolverCategory::overlapping,
                                   bool freecomm_ = false)
      : comm(comm_), cc(comm_), pis(), ri(pis,pis,comm_),
        OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
        OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
        CopyToAllInterfaceBuilt(false), globalLookup_(0), category_(cat_),
        freecomm(freecomm_)
    {}
 
    OwnerOverlapCopyCommunication (SolverCategory::Category cat_ = SolverCategory::overlapping)
      : comm(MPI_COMM_WORLD), cc(MPI_COMM_WORLD), pis(), ri(pis,pis,MPI_COMM_WORLD),
        OwnerToAllInterfaceBuilt(false), OwnerOverlapToAllInterfaceBuilt(false),
        OwnerCopyToAllInterfaceBuilt(false), OwnerCopyToOwnerCopyInterfaceBuilt(false),
        CopyToAllInterfaceBuilt(false), globalLookup_(0), category_(cat_), freecomm(false)
    {}
 
    OwnerOverlapCopyCommunication (const IndexInfoFromGrid<GlobalIdType, LocalIdType>& indexinfo,
                                   MPI_Comm comm_,
                                   SolverCategory::Category cat_ = SolverCategory::overlapping,
                                   bool freecomm_ = false)
      : comm(comm_), cc(comm_), OwnerToAllInterfaceBuilt(false),
        OwnerOverlapToAllInterfaceBuilt(false), OwnerCopyToAllInterfaceBuilt(false),
        OwnerCopyToOwnerCopyInterfaceBuilt(false), CopyToAllInterfaceBuilt(false),
        globalLookup_(0), category_(cat_), freecomm(freecomm_)
    {
      // set up an ISTL index set
      pis.beginResize();
      for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i)
      {
        if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::owner)
          pis.add(std::get<0>(*i),LI(std::get<1>(*i),OwnerOverlapCopyAttributeSet::owner,true));
        if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::overlap)
          pis.add(std::get<0>(*i),LI(std::get<1>(*i),OwnerOverlapCopyAttributeSet::overlap,true));
        if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::copy)
          pis.add(std::get<0>(*i),LI(std::get<1>(*i),OwnerOverlapCopyAttributeSet::copy,true));
        //                std::cout << cc.rank() << ": adding index " << std::get<0>(*i) << " " << std::get<1>(*i) << " " << std::get<2>(*i) << std::endl;
      }
      pis.endResize();
 
      // build remote indices WITHOUT communication
      //          std::cout << cc.rank() << ": build remote indices" << std::endl;
      ri.setIndexSets(pis,pis,cc);
      if (indexinfo.remoteIndices().size()>0)
      {
        remoteindex_iterator i=indexinfo.remoteIndices().begin();
        int p = std::get<0>(*i);
        RILM modifier = ri.template getModifier<false,true>(p);
        typename PIS::const_iterator pi=pis.begin();
        for ( ; i!=indexinfo.remoteIndices().end(); ++i)
        {
          // handle processor change
          if (p!=std::get<0>(*i))
          {
            p = std::get<0>(*i);
            modifier = ri.template getModifier<false,true>(p);
            pi=pis.begin();
          }
 
          // position to correct entry in parallel index set
          while (pi->global()!=std::get<1>(*i) && pi!=pis.end())
            ++pi;
          if (pi==pis.end())
            DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set");
 
          // insert entry
          //                      std::cout << cc.rank() << ": adding remote index " << std::get<0>(*i) << " " << std::get<1>(*i) << " " << std::get<2>(*i) << std::endl;
          if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::owner)
            modifier.insert(RX(OwnerOverlapCopyAttributeSet::owner,&(*pi)));
          if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::overlap)
            modifier.insert(RX(OwnerOverlapCopyAttributeSet::overlap,&(*pi)));
          if (std::get<2>(*i)==OwnerOverlapCopyAttributeSet::copy)
            modifier.insert(RX(OwnerOverlapCopyAttributeSet::copy,&(*pi)));
        }
      }else{
        // Force remote indices to be synced!
        ri.template getModifier<false,true>(0);
      }
    }
 
    // destructor: free memory in some objects
    ~OwnerOverlapCopyCommunication ()
    {
      ri.free();
      if (OwnerToAllInterfaceBuilt) OwnerToAllInterface.free();
      if (OwnerOverlapToAllInterfaceBuilt) OwnerOverlapToAllInterface.free();
      if (OwnerCopyToAllInterfaceBuilt) OwnerCopyToAllInterface.free();
      if (OwnerCopyToOwnerCopyInterfaceBuilt) OwnerCopyToOwnerCopyInterface.free();
      if (CopyToAllInterfaceBuilt) CopyToAllInterface.free();
      if (globalLookup_) delete globalLookup_;
      if (freecomm==true)
        if(comm!=MPI_COMM_NULL)
        {
#ifdef MPI_2
          // If it is possible to query whether MPI_Finalize
          // was called, only free the communicator before
          // calling MPI_Finalize.
          int wasFinalized = 0;
          MPI_Finalized( &wasFinalized );
          if(!wasFinalized)
#endif
          MPI_Comm_free(&comm);
        }
    }
 
  private:
    OwnerOverlapCopyCommunication (const OwnerOverlapCopyCommunication&)
    {}
    MPI_Comm comm;
    CollectiveCommunication<MPI_Comm> cc;
    PIS pis;
    RI ri;
    mutable IF OwnerToAllInterface;
    mutable bool OwnerToAllInterfaceBuilt;
    mutable IF OwnerOverlapToAllInterface;
    mutable bool OwnerOverlapToAllInterfaceBuilt;
    mutable IF OwnerCopyToAllInterface;
    mutable bool OwnerCopyToAllInterfaceBuilt;
    mutable IF OwnerCopyToOwnerCopyInterface;
    mutable bool OwnerCopyToOwnerCopyInterfaceBuilt;
    mutable IF CopyToAllInterface;
    mutable bool CopyToAllInterfaceBuilt;
    mutable std::vector<double> mask;
    int oldseqNo;
    GlobalLookupIndexSet* globalLookup_;
    // re-intruduce const qualifier once deprecated setCategory got removed after Dune 2.7
    /* const */ SolverCategory::Category category_;
    bool freecomm;
  };
 
#endif
 
 
} // end namespace
 
#endif