indexset.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:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_COMMON_PARALLEL_INDEXSET_HH
#define DUNE_COMMON_PARALLEL_INDEXSET_HH

#include <algorithm>
#include <cstdint> // for uint32_t
#include <iostream>

#include <dune/common/arraylist.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/localindex.hh>
#include <dune/common/parallel/mpitraits.hh>

namespace Dune
{
  // forward declarations

  template<class TG, class TL>
  class IndexPair;

  template<class TG, class TL>
  std::ostream& operator<<(std::ostream& os, const IndexPair<TG,TL>& pair);

  template<class TG, class TL>
  bool operator==(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator!=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator<(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator<=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator >=(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);

  template<class TG, class TL>
  bool operator==(const IndexPair<TG,TL>&, const TG&);

  template<class TG, class TL>
  bool operator!=(const IndexPair<TG,TL>&, const TG&);

  template<class TG, class TL>
  bool operator<(const IndexPair<TG,TL>&, const TG&);

  template<class TG, class TL>
  bool operator>(const IndexPair<TG,TL>&, const TG&);

  template<class TG, class TL>
  bool operator<=(const IndexPair<TG,TL>&, const TG&);

  template<class TG, class TL>
  bool operator >=(const IndexPair<TG,TL>&, const TG&);

  template<typename T>
  struct MPITraits;

  template<class TG, class TL>
  class IndexPair
  {
    friend std::ostream& operator<< <>(std::ostream&, const IndexPair<TG,TL>&);
    friend bool operator==<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator!=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator< <>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator><>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator<=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator>=<>(const IndexPair<TG,TL>&, const IndexPair<TG,TL>&);
    friend bool operator==<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator!=<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator< <>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator> <>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator<=<>(const IndexPair<TG,TL>&, const TG &);
    friend bool operator>=<>(const IndexPair<TG,TL>&, const TG &);
    friend struct MPITraits<IndexPair<TG,TL> >;

  public:
    typedef TG GlobalIndex;

    typedef TL LocalIndex;

    IndexPair(const GlobalIndex& global, const LocalIndex& local);

    IndexPair();
    IndexPair(const GlobalIndex& global);

    inline const GlobalIndex& global() const;

    inline LocalIndex& local();

    inline const LocalIndex& local() const;

    inline void setLocal(int index);
  private:
    GlobalIndex global_;
    LocalIndex local_;
  };

  enum ParallelIndexSetState
  {
    GROUND,
    RESIZE
  };

  class InvalidIndexSetState : public InvalidStateException {};

  // Forward declaration
  template<class I> class GlobalLookupIndexSet;

  template<typename TG, typename TL, int N=100>
  class ParallelIndexSet
  {
    friend class GlobalLookupIndexSet<ParallelIndexSet<TG,TL,N> >;

  public:
    typedef TG GlobalIndex;

    typedef TL LocalIndex;

    typedef Dune::IndexPair<GlobalIndex,LocalIndex> IndexPair;

    constexpr static int arraySize = (N>0) ? N : 1;

    class iterator :
      public ArrayList<IndexPair,N>::iterator
    {
      typedef typename ArrayList<IndexPair,N>::iterator
      Father;
      friend class ParallelIndexSet<GlobalIndex,LocalIndex,N>;
    public:
      iterator(ParallelIndexSet<TG,TL,N>& indexSet, const Father& father)
        : Father(father), indexSet_(&indexSet)
      {}

    private:
      inline void markAsDeleted() const
      {
#ifndef NDEBUG
        if(indexSet_->state_ != RESIZE)
          DUNE_THROW(InvalidIndexSetState, "Indices can only be removed "
                     <<"while in RESIZE state!");
#endif
        Father::operator*().local().setState(DELETED);
      }

      ParallelIndexSet<TG,TL,N>* indexSet_;

    };



    typedef typename
    ArrayList<IndexPair,N>::const_iterator
    const_iterator;

    ParallelIndexSet();

    inline const ParallelIndexSetState& state()
    {
      return state_;
    }

    void beginResize();

    inline void add(const GlobalIndex& global);

    inline void add(const GlobalIndex& global, const LocalIndex& local);

    inline void markAsDeleted(const iterator& position);

    void endResize();

    inline IndexPair&
    operator[](const GlobalIndex& global);

    inline IndexPair&
    at(const GlobalIndex& global);

    inline bool
    exists (const GlobalIndex& global) const;

    inline const IndexPair&
    operator[](const GlobalIndex& global) const;

    inline const IndexPair&
    at(const GlobalIndex& global) const;

    inline iterator begin();

    inline iterator end();

    inline const_iterator begin() const;

    inline const_iterator end() const;

    inline void renumberLocal();

    inline int seqNo() const;

    inline size_t size() const;

  private:
    ArrayList<IndexPair,N> localIndices_;
    ArrayList<IndexPair,N> newIndices_;
    ParallelIndexSetState state_;
    int seqNo_;
    bool deletedEntries_;
    inline void merge();
  };


  template<class TG, class TL, int N>
  std::ostream& operator<<(std::ostream& os, const ParallelIndexSet<TG,TL,N>& indexSet);

  template<class I>
  class GlobalLookupIndexSet
  {
  public:
    typedef I ParallelIndexSet;

    typedef typename ParallelIndexSet::LocalIndex LocalIndex;

    typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;

    typedef typename ParallelIndexSet::const_iterator const_iterator;

    typedef Dune::IndexPair<typename I::GlobalIndex, typename I::LocalIndex> IndexPair;

    GlobalLookupIndexSet(const ParallelIndexSet& indexset, std::size_t size);

    GlobalLookupIndexSet(const ParallelIndexSet& indexset);

    ~GlobalLookupIndexSet();

    inline const IndexPair&
    operator[](const GlobalIndex& global) const;

    inline const IndexPair*
    pair(const std::size_t& local) const;

    inline const_iterator begin() const;

    inline const_iterator end() const;

    inline int seqNo() const;

    inline size_t size() const;
  private:
    const ParallelIndexSet& indexSet_;

    std::size_t size_;

    std::vector<const IndexPair*> indices_;

  };


  template<typename T>
  struct LocalIndexComparator
  {
    static bool compare([[maybe_unused]] const T& t1, [[maybe_unused]] const T& t2)
    {
      return false;
    }
  };

  template<class TG, class TL>
  struct IndexSetSortFunctor
  {
    bool operator()(const IndexPair<TG,TL>& i1, const IndexPair<TG,TL>& i2)
    {
      return i1.global()<i2.global() || (i1.global()==i2.global() &&
                                         LocalIndexComparator<TL>::compare(i1.local(),
                                                                           i2.local()));
    }
  };



  template<class TG, class TL>
  inline std::ostream& operator<<(std::ostream& os, const IndexPair<TG,TL>& pair)
  {
    os<<"{global="<<pair.global_<<", local="<<pair.local_<<"}";
    return os;
  }

  template<class TG, class TL, int N>
  inline std::ostream& operator<<(std::ostream& os, const ParallelIndexSet<TG,TL,N>& indexSet)
  {
    typedef typename ParallelIndexSet<TG,TL,N>::const_iterator Iterator;
    Iterator end = indexSet.end();
    os<<"{";
    for(Iterator index = indexSet.begin(); index != end; ++index)
      os<<*index<<" ";
    os<<"}";
    return os;

  }

  template<class TG, class TL>
  inline bool operator==(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_==b.global_;
  }

  template<class TG, class TL>
  inline bool operator!=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_!=b.global_;
  }

  template<class TG, class TL>
  inline bool operator<(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_<b.global_;
  }

  template<class TG, class TL>
  inline bool operator>(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_>b.global_;
  }

  template<class TG, class TL>
  inline bool operator<=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_<=b.global_;
  }

  template<class TG, class TL>
  inline bool operator >=(const IndexPair<TG,TL>& a, const IndexPair<TG,TL>& b)
  {
    return a.global_>=b.global_;
  }

  template<class TG, class TL>
  inline bool operator==(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_==b;
  }

  template<class TG, class TL>
  inline bool operator!=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_!=b;
  }

  template<class TG, class TL>
  inline bool operator<(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_<b;
  }

  template<class TG, class TL>
  inline bool operator>(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_>b;
  }

  template<class TG, class TL>
  inline bool operator<=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_<=b;
  }

  template<class TG, class TL>
  inline bool operator >=(const IndexPair<TG,TL>& a, const TG& b)
  {
    return a.global_>=b;
  }

#ifndef DOXYGEN

  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair(const TG& global, const TL& local)
    : global_(global), local_(local){}

  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair(const TG& global)
    : global_(global), local_(){}

  template<class TG, class TL>
  IndexPair<TG,TL>::IndexPair()
    : global_(), local_(){}

  template<class TG, class TL>
  inline const TG& IndexPair<TG,TL>::global() const {
    return global_;
  }

  template<class TG, class TL>
  inline TL& IndexPair<TG,TL>::local() {
    return local_;
  }

  template<class TG, class TL>
  inline const TL& IndexPair<TG,TL>::local() const {
    return local_;
  }

  template<class TG, class TL>
  inline void IndexPair<TG,TL>::setLocal(int local){
    local_=local;
  }

  template<class TG, class TL, int N>
  ParallelIndexSet<TG,TL,N>::ParallelIndexSet()
    : state_(GROUND), seqNo_(0), deletedEntries_()
  {}

  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::beginResize()
  {

    // Checks in unproductive code
#ifndef NDEBUG
    if(state_!=GROUND)
      DUNE_THROW(InvalidIndexSetState,
                 "IndexSet has to be in GROUND state, when "
                 << "beginResize() is called!");
#endif

    state_ = RESIZE;
    deletedEntries_ = false;
  }

  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::add(const GlobalIndex& global)
  {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be added "
                 <<"while in RESIZE state!");
#endif
    newIndices_.push_back(IndexPair(global));
  }

  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::add(const TG& global, const TL& local)
  {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be added "
                 <<"while in RESIZE state!");
#endif
    newIndices_.push_back(IndexPair(global,local));
  }

  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::markAsDeleted(const iterator& global)
  {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "Indices can only be removed "
                 <<"while in RESIZE state!");
#endif
    deletedEntries_ = true;

    global.markAsDeleted();
  }

  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::endResize() {
    // Checks in unproductive code
#ifndef NDEBUG
    if(state_ != RESIZE)
      DUNE_THROW(InvalidIndexSetState, "endResize called while not "
                 <<"in RESIZE state!");
#endif

    std::sort(newIndices_.begin(), newIndices_.end(), IndexSetSortFunctor<TG,TL>());
    merge();
    seqNo_++;
    state_ = GROUND;
  }


  template<class TG, class TL, int N>
  inline void ParallelIndexSet<TG,TL,N>::merge(){
    if(localIndices_.size()==0)
    {
      localIndices_=newIndices_;
      newIndices_.clear();
    }
    else if(newIndices_.size()>0 || deletedEntries_)
    {
      ArrayList<IndexPair,N> tempPairs;

      auto old = localIndices_.begin();
      auto added = newIndices_.begin();
      const auto endold = localIndices_.end();
      const auto endadded = newIndices_.end();

      while(old != endold && added!= endadded)
      {
        if(old->local().state()==DELETED) {
          old.eraseToHere();
        }
        else
        {
          if(old->global() < added->global() ||
             (old->global() == added->global()
              && LocalIndexComparator<TL>::compare(old->local(),added->local())))
          {
            tempPairs.push_back(*old);
            old.eraseToHere();
            continue;
          }else
          {
            tempPairs.push_back(*added);
            added.eraseToHere();
          }
        }
      }

      while(old != endold)
      {
        if(old->local().state()!=DELETED) {
          tempPairs.push_back(*old);
        }
        old.eraseToHere();
      }

      while(added!= endadded)
      {
        tempPairs.push_back(*added);
        added.eraseToHere();
      }
      localIndices_ = tempPairs;
    }
  }


  template<class TG, class TL, int N>
  inline const IndexPair<TG,TL>&
  ParallelIndexSet<TG,TL,N>::at(const TG& global) const
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;

    while(low<high)
    {
      probe = (high + low) / 2;
      if(global <= localIndices_[probe].global())
        high = probe;
      else
        low = probe+1;
    }

    if(probe==-1)
      DUNE_THROW(RangeError, "No entries!");

    if( localIndices_[low].global() != global)
      DUNE_THROW(RangeError, "Could not find entry of "<<global);
    else
      return localIndices_[low];
  }

  template<class TG, class TL, int N>
  inline const IndexPair<TG,TL>&
  ParallelIndexSet<TG,TL,N>::operator[](const TG& global) const
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;

    while(low<high)
    {
      probe = (high + low) / 2;
      if(global <= localIndices_[probe].global())
        high = probe;
      else
        low = probe+1;
    }

    return localIndices_[low];
  }
  template<class TG, class TL, int N>
  inline IndexPair<TG,TL>& ParallelIndexSet<TG,TL,N>::at(const TG& global)
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;

    while(low<high)
    {
      probe = (high + low) / 2;
      if(localIndices_[probe].global() >= global)
        high = probe;
      else
        low = probe+1;
    }

    if(probe==-1)
      DUNE_THROW(RangeError, "No entries!");

    if( localIndices_[low].global() != global)
      DUNE_THROW(RangeError, "Could not find entry of "<<global);
    else
      return localIndices_[low];
  }

  template<class TG, class TL, int N>
  inline bool ParallelIndexSet<TG,TL,N>::exists (const TG& global) const
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;

    while(low<high)
    {
      probe = (high + low) / 2;
      if(localIndices_[probe].global() >= global)
        high = probe;
      else
        low = probe+1;
    }

    if(probe==-1)
      return false;

    if( localIndices_[low].global() != global)
      return false;
    return true;
  }

  template<class TG, class TL, int N>
  inline IndexPair<TG,TL>& ParallelIndexSet<TG,TL,N>::operator[](const TG& global)
  {
    // perform a binary search
    int low=0, high=localIndices_.size()-1, probe=-1;

    while(low<high)
    {
      probe = (high + low) / 2;
      if(localIndices_[probe].global() >= global)
        high = probe;
      else
        low = probe+1;
    }

    return localIndices_[low];
  }
  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::iterator
  ParallelIndexSet<TG,TL,N>::begin()
  {
    return iterator(*this, localIndices_.begin());
  }


  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::iterator
  ParallelIndexSet<TG,TL,N>::end()
  {
    return iterator(*this,localIndices_.end());
  }

  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::const_iterator
  ParallelIndexSet<TG,TL,N>::begin() const
  {
    return localIndices_.begin();
  }


  template<class TG, class TL, int N>
  inline typename ParallelIndexSet<TG,TL,N>::const_iterator
  ParallelIndexSet<TG,TL,N>::end() const
  {
    return localIndices_.end();
  }

  template<class TG, class TL, int N>
  void ParallelIndexSet<TG,TL,N>::renumberLocal(){
#ifndef NDEBUG
    if(state_==RESIZE)
      DUNE_THROW(InvalidIndexSetState, "IndexSet has to be in "
                 <<"GROUND state for renumberLocal()");
#endif

    const auto end_ = end();
    uint32_t index=0;

    for(auto pair=begin(); pair!=end_; index++, ++pair)
      pair->local()=index;
  }

  template<class TG, class TL, int N>
  inline int ParallelIndexSet<TG,TL,N>::seqNo() const
  {
    return seqNo_;
  }

  template<class TG, class TL, int N>
  inline size_t ParallelIndexSet<TG,TL,N>::size() const
  {
    return localIndices_.size();
  }

  template<class I>
  GlobalLookupIndexSet<I>::GlobalLookupIndexSet(const I& indexset,
                                                std::size_t size)
    : indexSet_(indexset), size_(size),
      indices_(size_, static_cast<const IndexPair*>(0))
  {
    const_iterator end_ = indexSet_.end();

    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair) {
      assert(pair->local()<size_);
      indices_[pair->local()] = &(*pair);
    }
  }

  template<class I>
  GlobalLookupIndexSet<I>::GlobalLookupIndexSet(const I& indexset)
    : indexSet_(indexset), size_(0)
  {
    const_iterator end_ = indexSet_.end();
    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair)
      size_=std::max<std::size_t>(size_,pair->local());

    indices_.resize(++size_,  0);

    for(const_iterator pair = indexSet_.begin(); pair!=end_; ++pair)
      indices_[pair->local()] = &(*pair);
  }

  template<class I>
  GlobalLookupIndexSet<I>::~GlobalLookupIndexSet()
  {}

  template<class I>
  inline const IndexPair<typename I::GlobalIndex, typename I::LocalIndex>*
  GlobalLookupIndexSet<I>::pair(const std::size_t& local) const
  {
    return indices_[local];
  }

  template<class I>
  inline const IndexPair<typename I::GlobalIndex, typename I::LocalIndex>&
  GlobalLookupIndexSet<I>::operator[](const GlobalIndex& global) const
  {
    return indexSet_[global];
  }

  template<class I>
  typename I::const_iterator GlobalLookupIndexSet<I>::begin() const
  {
    return indexSet_.begin();
  }

  template<class I>
  typename I::const_iterator GlobalLookupIndexSet<I>::end() const
  {
    return indexSet_.end();
  }

  template<class I>
  inline size_t GlobalLookupIndexSet<I>::size() const
  {
    return size_;
  }

  template<class I>
  inline int GlobalLookupIndexSet<I>::seqNo() const
  {
    return indexSet_.seqNo();
  }

  template<typename TG, typename TL, int N, typename TG1, typename TL1, int N1>
  bool operator==(const ParallelIndexSet<TG,TL,N>& idxset,
                  const ParallelIndexSet<TG1,TL1,N1>& idxset1)
  {
    if(idxset.size()!=idxset1.size())
      return false;
    typedef typename ParallelIndexSet<TG,TL,N>::const_iterator Iter;
    typedef typename ParallelIndexSet<TG1,TL1,N1>::const_iterator Iter1;
    Iter iter=idxset.begin();
    for(Iter1 iter1=idxset1.begin(); iter1 != idxset1.end(); ++iter, ++iter1) {
      if(iter1->global()!=iter->global())
        return false;
      typedef typename ParallelIndexSet<TG,TL,N>::LocalIndex PI;
      const PI& pi=iter->local(), pi1=iter1->local();

      if(pi!=pi1)
        return false;
    }
    return true;
  }

  template<typename TG, typename TL, int N, typename TG1, typename TL1, int N1>
  bool operator!=(const ParallelIndexSet<TG,TL,N>& idxset,
                  const ParallelIndexSet<TG1,TL1,N1>& idxset1)
  {
    return !(idxset==idxset1);
  }


#endif // DOXYGEN

}

#endif // DUNE_COMMON_PARALLEL_INDEXSET_HH