indicessyncer.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_INDICESSYNCER_HH
#define DUNE_INDICESSYNCER_HH
 
#include "indexset.hh"
#include "remoteindices.hh"
#include <dune/common/stdstreams.hh>
#include <dune/common/tuples.hh>
#include <dune/common/sllist.hh>
#include <dune/common/unused.hh>
#include <cassert>
#include <cmath>
#include <limits>
#include <algorithm>
#include <functional>
#include <map>
 
#if HAVE_MPI
namespace Dune
{
  template<typename T>
  class IndicesSyncer
  {
  public:
 
    typedef T ParallelIndexSet;
 
    typedef typename ParallelIndexSet::IndexPair IndexPair;
 
    typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
 
    typedef typename ParallelIndexSet::LocalIndex::Attribute Attribute;
 
    typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
 
    IndicesSyncer(ParallelIndexSet& indexSet,
                  RemoteIndices& remoteIndices);
 
    void sync();
 
    template<typename T1>
    void sync(T1& numberer);
 
  private:
 
    ParallelIndexSet& indexSet_;
 
    RemoteIndices& remoteIndices_;
 
    char** sendBuffers_;
 
    char* receiveBuffer_;
 
    std::size_t* sendBufferSizes_;
 
    int receiveBufferSize_; // int because of MPI
 
    struct MessageInformation
    {
      MessageInformation()
        : publish(), pairs()
      {}
      int publish;
      int pairs;
    };
 
    class DefaultNumberer
    {
    public:
      std::size_t operator()(const GlobalIndex& global)
      {
        DUNE_UNUSED_PARAMETER(global);
        return std::numeric_limits<size_t>::max();
      }
    };
 
    MPI_Datatype datatype_;
 
    int rank_;
 
    typedef SLList<std::pair<GlobalIndex,Attribute>, typename RemoteIndices::Allocator> GlobalIndexList;
 
    typedef typename GlobalIndexList::ModifyIterator GlobalIndexModifier;
 
    typedef typename SLList<GlobalIndex, typename RemoteIndices::Allocator>::iterator
    GlobalIndexIterator;
 
    typedef std::map<int, GlobalIndexList> GlobalIndicesMap;
 
    GlobalIndicesMap globalMap_;
 
    typedef SLList<bool, typename RemoteIndices::Allocator> BoolList;
 
    typedef typename BoolList::iterator BoolIterator;
 
    typedef typename BoolList::ModifyIterator BoolListModifier;
 
    typedef std::map<int,BoolList> BoolMap;
 
    BoolMap oldMap_;
 
    std::map<int,MessageInformation> infoSend_;
 
    typedef typename RemoteIndices::RemoteIndexList RemoteIndexList;
 
    typedef typename RemoteIndexList::ModifyIterator RemoteIndexModifier;
 
    typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
 
    typedef typename RemoteIndexList::iterator RemoteIndexIterator;
 
    typedef typename RemoteIndexList::const_iterator ConstRemoteIndexIterator;
 
    typedef Dune::tuple<RemoteIndexModifier,GlobalIndexModifier,BoolListModifier,
        const ConstRemoteIndexIterator> IteratorTuple;
 
    class Iterators
    {
      friend class IndicesSyncer<T>;
    public:
      Iterators(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
                BoolList& booleans);
 
      Iterators();
 
      Iterators& operator++();
 
      void insert(const RemoteIndex& index,
                  const std::pair<GlobalIndex,Attribute>& global);
 
      RemoteIndex& remoteIndex() const;
 
      std::pair<GlobalIndex,Attribute>& globalIndexPair() const;
 
      Attribute& attribute() const;
 
      bool isOld() const;
 
      void reset(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
                 BoolList& booleans);
 
      bool isNotAtEnd() const;
 
      bool isAtEnd() const;
 
    private:
      IteratorTuple iterators_;
    };
 
    typedef std::map<int,Iterators> IteratorsMap;
 
    IteratorsMap iteratorsMap_;
 
    void calculateMessageSizes();
 
    void packAndSend(int destination, char* buffer, std::size_t bufferSize, MPI_Request& req);
 
    template<typename T1>
    void recvAndUnpack(T1& numberer);
 
    void registerMessageDatatype();
 
    void insertIntoRemoteIndexList(int process,
                                   const std::pair<GlobalIndex,Attribute>& global,
                                   char attribute);
 
    void resetIteratorsMap();
 
    bool checkReset();
 
    bool checkReset(const Iterators& iterators, RemoteIndexList& rlist, GlobalIndexList& gList,
                    BoolList& bList);
  };
 
  template<typename TG, typename TA>
  bool operator<(const IndexPair<TG,ParallelLocalIndex<TA> >& i1,
                 const std::pair<TG,TA>& i2)
  {
    return i1.global() < i2.first ||
           (i1.global() == i2.first && i1.local().attribute()<i2.second);
  }
 
  template<typename TG, typename TA>
  bool operator<(const std::pair<TG,TA>& i1,
                 const IndexPair<TG,ParallelLocalIndex<TA> >& i2)
  {
    return i1.first < i2.global() ||
           (i1.first == i2.global() && i1.second<i2.local().attribute());
  }
 
  template<typename TG, typename TA>
  bool operator==(const IndexPair<TG,ParallelLocalIndex<TA> >& i1,
                  const std::pair<TG,TA>& i2)
  {
    return (i1.global() == i2.first && i1.local().attribute()==i2.second);
  }
 
  template<typename TG, typename TA>
  bool operator!=(const IndexPair<TG,ParallelLocalIndex<TA> >& i1,
                  const std::pair<TG,TA>& i2)
  {
    return (i1.global() != i2.first || i1.local().attribute()!=i2.second);
  }
 
  template<typename TG, typename TA>
  bool operator==(const std::pair<TG,TA>& i2,
                  const IndexPair<TG,ParallelLocalIndex<TA> >& i1)
  {
    return (i1.global() == i2.first && i1.local().attribute()==i2.second);
  }
 
  template<typename TG, typename TA>
  bool operator!=(const std::pair<TG,TA>& i2,
                  const IndexPair<TG,ParallelLocalIndex<TA> >& i1)
  {
    return (i1.global() != i2.first || i1.local().attribute()!=i2.second);
  }
 
  template<typename T, typename A, typename A1>
  void storeGlobalIndicesOfRemoteIndices(std::map<int,SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> >& globalMap,
                                         const RemoteIndices<T,A1>& remoteIndices)
  {
    typedef typename RemoteIndices<T,A1>::const_iterator RemoteIterator;
 
    for(RemoteIterator remote = remoteIndices.begin(), end =remoteIndices.end(); remote != end; ++remote) {
      typedef typename RemoteIndices<T,A1>::RemoteIndexList RemoteIndexList;
      typedef typename RemoteIndexList::const_iterator RemoteIndexIterator;
      typedef SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> GlobalIndexList;
      GlobalIndexList& global = globalMap[remote->first];
      RemoteIndexList& rList = *(remote->second.first);
 
      for(RemoteIndexIterator index = rList.begin(), riEnd = rList.end();
          index != riEnd; ++index) {
        global.push_back(std::make_pair(index->localIndexPair().global(),
                                        index->localIndexPair().local().attribute()));
      }
    }
  }
 
  template<typename T, typename A, typename A1>
  inline void repairLocalIndexPointers(std::map<int,
                                           SLList<std::pair<typename T::GlobalIndex,
                                                   typename T::LocalIndex::Attribute>,A> >& globalMap,
                                       RemoteIndices<T,A1>& remoteIndices,
                                       const T& indexSet)
  {
    typedef typename RemoteIndices<T,A1>::RemoteIndexMap::iterator RemoteIterator;
    typedef typename RemoteIndices<T,A1>::RemoteIndexList::iterator RemoteIndexIterator;
    typedef typename T::GlobalIndex GlobalIndex;
    typedef typename T::LocalIndex::Attribute Attribute;
    typedef std::pair<GlobalIndex,Attribute> GlobalIndexPair;
    typedef SLList<GlobalIndexPair,A> GlobalIndexPairList;
    typedef typename GlobalIndexPairList::iterator GlobalIndexIterator;
 
    assert(globalMap.size()==static_cast<std::size_t>(remoteIndices.neighbours()));
    // Repair pointers to index set in remote indices.
    typename std::map<int,GlobalIndexPairList>::iterator global = globalMap.begin();
    RemoteIterator end = remoteIndices.remoteIndices_.end();
 
    for(RemoteIterator remote = remoteIndices.remoteIndices_.begin(); remote != end; ++remote, ++global) {
      typedef typename T::const_iterator IndexIterator;
 
      assert(remote->first==global->first);
      assert(remote->second.first->size() == global->second.size());
 
      RemoteIndexIterator riEnd  = remote->second.first->end();
      RemoteIndexIterator rIndex = remote->second.first->begin();
      GlobalIndexIterator gIndex = global->second.begin();
      IndexIterator index  = indexSet.begin();
 
      assert(rIndex==riEnd || gIndex != global->second.end());
      while(rIndex != riEnd) {
        // Search for the index in the set.
        assert(gIndex != global->second.end());
 
        while(!(index->global() == gIndex->first
                && index->local().attribute() == gIndex->second)) {
          ++index;
          // this is only needed for ALU, where there may exist
          // more entries with the same global index in the remote index set
          // than in the index set
          if (index->global() > gIndex->first) {
            index=indexSet.begin();
          }
        }
 
        assert(index != indexSet.end() && *index == *gIndex);
 
        rIndex->localIndex_ = &(*index);
        ++index;
        ++rIndex;
        ++gIndex;
      }
    }
    remoteIndices.sourceSeqNo_ = remoteIndices.source_->seqNo();
    remoteIndices.destSeqNo_ = remoteIndices.target_->seqNo();
  }
 
  template<typename T>
  IndicesSyncer<T>::IndicesSyncer(ParallelIndexSet& indexSet,
                                  RemoteIndices& remoteIndices)
    : indexSet_(indexSet), remoteIndices_(remoteIndices)
  {
    // index sets must match.
    assert(remoteIndices.source_ == remoteIndices.target_);
    assert(remoteIndices.source_ == &indexSet);
    MPI_Comm_rank(remoteIndices_.communicator(), &rank_);
  }
 
  template<typename T>
  IndicesSyncer<T>::Iterators::Iterators(RemoteIndexList& remoteIndices,
                                         GlobalIndexList& globalIndices,
                                         BoolList& booleans)
    : iterators_(remoteIndices.beginModify(), globalIndices.beginModify(),
                 booleans.beginModify(), remoteIndices.end())
  { }
 
  template<typename T>
  IndicesSyncer<T>::Iterators::Iterators()
    : iterators_()
  {}
 
  template<typename T>
  inline typename IndicesSyncer<T>::Iterators& IndicesSyncer<T>::Iterators::operator++()
  {
    ++(get<0>(iterators_));
    ++(get<1>(iterators_));
    ++(get<2>(iterators_));
    return *this;
  }
 
  template<typename T>
  inline void IndicesSyncer<T>::Iterators::insert(const RemoteIndex & index,
                                                  const std::pair<GlobalIndex,Attribute>& global)
  {
    get<0>(iterators_).insert(index);
    get<1>(iterators_).insert(global);
    get<2>(iterators_).insert(false);
  }
 
  template<typename T>
  inline typename IndicesSyncer<T>::RemoteIndex&
  IndicesSyncer<T>::Iterators::remoteIndex() const
  {
    return *(get<0>(iterators_));
  }
 
  template<typename T>
  inline std::pair<typename IndicesSyncer<T>::GlobalIndex,typename IndicesSyncer<T>::Attribute>&
  IndicesSyncer<T>::Iterators::globalIndexPair() const
  {
    return *(get<1>(iterators_));
  }
 
  template<typename T>
  inline bool IndicesSyncer<T>::Iterators::isOld() const
  {
    return *(get<2>(iterators_));
  }
 
  template<typename T>
  inline void IndicesSyncer<T>::Iterators::reset(RemoteIndexList& remoteIndices,
                                                 GlobalIndexList& globalIndices,
                                                 BoolList& booleans)
  {
    get<0>(iterators_) = remoteIndices.beginModify();
    get<1>(iterators_) = globalIndices.beginModify();
    get<2>(iterators_) = booleans.beginModify();
  }
 
  template<typename T>
  inline bool IndicesSyncer<T>::Iterators::isNotAtEnd() const
  {
    return get<0>(iterators_)!=get<3>(iterators_);
  }
 
  template<typename T>
  inline bool IndicesSyncer<T>::Iterators::isAtEnd() const
  {
    return get<0>(iterators_)==get<3>(iterators_);
  }
 
  template<typename T>
  void IndicesSyncer<T>::registerMessageDatatype()
  {
    MPI_Datatype type[2] = {MPI_INT, MPI_INT};
    int blocklength[2] = {1,1};
    MPI_Aint displacement[2];
    MPI_Aint base;
 
    // Compute displacement
    MessageInformation message;
 
    MPI_Get_address( &(message.publish), displacement);
    MPI_Get_address( &(message.pairs), displacement+1);
 
    // Make the displacement relative
    MPI_Get_address(&message, &base);
    displacement[0] -= base;
    displacement[1] -= base;
 
    MPI_Type_create_struct( 2, blocklength, displacement, type, &datatype_);
    MPI_Type_commit(&datatype_);
  }
 
  template<typename T>
  void IndicesSyncer<T>::calculateMessageSizes()
  {
    typedef typename ParallelIndexSet::const_iterator IndexIterator;
    typedef CollectiveIterator<T,typename RemoteIndices::Allocator> CollectiveIterator;
 
    IndexIterator iEnd = indexSet_.end();
    CollectiveIterator collIter = remoteIndices_.template iterator<true>();
 
    for(IndexIterator index = indexSet_.begin(); index != iEnd; ++index) {
      collIter.advance(index->global(), index->local().attribute());
      if(collIter.empty())
        break;
      int knownRemote=0;
 
      typedef typename CollectiveIterator::iterator ValidIterator;
      ValidIterator end = collIter.end();
 
      // Count the remote indices we know.
      for(ValidIterator valid = collIter.begin(); valid != end; ++valid) {
        ++knownRemote;
      }
 
      if(knownRemote>0) {
        Dune::dverb<<rank_<<": publishing "<<knownRemote<<" for index "<<index->global()<< " for processes ";
 
        // Update MessageInformation
        for(ValidIterator valid = collIter.begin(); valid != end; ++valid) {
          ++(infoSend_[valid.process()].publish);
          (infoSend_[valid.process()].pairs) += knownRemote;
          Dune::dverb<<valid.process()<<" ";
          Dune::dverb<<"(publish="<<infoSend_[valid.process()].publish<<", pairs="<<infoSend_[valid.process()].pairs
                     <<") ";
        }
        Dune::dverb<<std::endl;
      }
    }
 
    typedef typename std::map<int,MessageInformation>::const_iterator
    MessageIterator;
 
    const MessageIterator end = infoSend_.end();
 
    // registerMessageDatatype();
 
    // Now determine the buffersizes needed for each neighbour using MPI_Pack_size
    MessageInformation dummy;
 
    MessageIterator messageIter= infoSend_.begin();
 
    typedef typename RemoteIndices::RemoteIndexMap::const_iterator RemoteIterator;
    const RemoteIterator rend = remoteIndices_.end();
    int neighbour=0;
 
    for(RemoteIterator remote = remoteIndices_.begin(); remote != rend; ++remote, ++neighbour) {
      MessageInformation* message;
      MessageInformation recv;
 
      if(messageIter != end && messageIter->first==remote->first) {
        // We want to send message information to that process
        message = const_cast<MessageInformation*>(&(messageIter->second));
        ++messageIter;
      }else
        // We do not want to send information but the other process might.
        message = &dummy;
 
      sendBufferSizes_[neighbour]=0;
      int tsize;
      // The number of indices published
      MPI_Pack_size(1, MPI_INT,remoteIndices_.communicator(), &tsize);
      sendBufferSizes_[neighbour] += tsize;
 
      for(int i=0; i < message->publish; ++i) {
        // The global index
        MPI_Pack_size(1, MPITraits<GlobalIndex>::getType(), remoteIndices_.communicator(), &tsize);
        sendBufferSizes_[neighbour] += tsize;
        // The attribute in the local index
        MPI_Pack_size(1, MPI_CHAR, remoteIndices_.communicator(), &tsize);
        sendBufferSizes_[neighbour] += tsize;
        // The number of corresponding remote indices
        MPI_Pack_size(1, MPI_INT, remoteIndices_.communicator(), &tsize);
        sendBufferSizes_[neighbour] += tsize;
      }
      for(int i=0; i < message->pairs; ++i) {
        // The process of the remote index
        MPI_Pack_size(1, MPI_INT, remoteIndices_.communicator(), &tsize);
        sendBufferSizes_[neighbour] += tsize;
        // The attribute of the remote index
        MPI_Pack_size(1, MPI_CHAR, remoteIndices_.communicator(), &tsize);
        sendBufferSizes_[neighbour] += tsize;
      }
 
      Dune::dverb<<rank_<<": Buffer (neighbour="<<remote->first<<") size is "<< sendBufferSizes_[neighbour]<<" for publish="<<message->publish<<" pairs="<<message->pairs<<std::endl;
    }
 
  }
 
  template<typename T>
  inline void IndicesSyncer<T>::sync()
  {
    DefaultNumberer numberer;
    sync(numberer);
  }
 
  template<typename T>
  template<typename T1>
  void IndicesSyncer<T>::sync(T1& numberer)
  {
 
    // The pointers to the local indices in the remote indices
    // will become invalid due to the resorting of the index set.
    // Therefore store the corresponding global indices.
    // Mark all indices as not added
 
    typedef typename RemoteIndices::RemoteIndexMap::const_iterator
    RemoteIterator;
 
    const RemoteIterator end = remoteIndices_.end();
 
    // Number of neighbours might change during the syncing.
    // save the old neighbours
    std::size_t noOldNeighbours = remoteIndices_.neighbours();
    int* oldNeighbours = new int[noOldNeighbours];
    sendBufferSizes_ = new std::size_t[noOldNeighbours];
    std::size_t neighbourI = 0;
 
    for(RemoteIterator remote = remoteIndices_.begin(); remote != end; ++remote, ++neighbourI) {
      typedef typename RemoteIndices::RemoteIndexList::const_iterator
      RemoteIndexIterator;
 
      oldNeighbours[neighbourI] = remote->first;
 
      // Make sure we only have one remote index list.
      assert(remote->second.first==remote->second.second);
 
      RemoteIndexList& rList = *(remote->second.first);
 
      // Store the corresponding global indices.
      GlobalIndexList& global = globalMap_[remote->first];
      BoolList& added = oldMap_[remote->first];
      RemoteIndexIterator riEnd = rList.end();
 
      for(RemoteIndexIterator index = rList.begin();
          index != riEnd; ++index) {
        global.push_back(std::make_pair(index->localIndexPair().global(),
                                        index->localIndexPair().local().attribute()));
        added.push_back(true);
      }
 
      Iterators iterators(rList, global, added);
      iteratorsMap_.insert(std::make_pair(remote->first, iterators));
      assert(checkReset(iteratorsMap_[remote->first], rList,global,added));
    }
 
    // Exchange indices with each neighbour
    calculateMessageSizes();
 
    // Allocate the buffers
    receiveBufferSize_=1;
    sendBuffers_ = new char*[noOldNeighbours];
 
    for(std::size_t i=0; i<noOldNeighbours; ++i) {
      sendBuffers_[i] = new char[sendBufferSizes_[i]];
      receiveBufferSize_ = std::max(receiveBufferSize_, static_cast<int>(sendBufferSizes_[i]));
    }
 
    receiveBuffer_=new char[receiveBufferSize_];
 
    indexSet_.beginResize();
 
    Dune::dverb<<rank_<<": Neighbours: ";
 
    for(std::size_t i = 0; i<noOldNeighbours; ++i)
      Dune::dverb<<oldNeighbours[i]<<" ";
 
    Dune::dverb<<std::endl;
 
    MPI_Request* requests = new MPI_Request[noOldNeighbours];
    MPI_Status* statuses = new MPI_Status[noOldNeighbours];
 
    // Pack Message data and start the sends
    for(std::size_t i = 0; i<noOldNeighbours; ++i)
      packAndSend(oldNeighbours[i], sendBuffers_[i], sendBufferSizes_[i], requests[i]);
 
    // Probe for incoming messages, receive and unpack them
    for(std::size_t i = 0; i<noOldNeighbours; ++i)
      recvAndUnpack(numberer);
    //       }else{
    //  recvAndUnpack(oldNeighbours[i], numberer);
    //  packAndSend(oldNeighbours[i]);
    //       }
    //     }
 
    delete[] receiveBuffer_;
 
    // Wait for the completion of the sends
    // Wait for completion of sends
    if(MPI_SUCCESS!=MPI_Waitall(noOldNeighbours, requests, statuses)) {
      std::cerr<<": MPI_Error occurred while sending message"<<std::endl;
      for(std::size_t i=0; i< noOldNeighbours; i++)
        if(MPI_SUCCESS!=statuses[i].MPI_ERROR)
          std::cerr<<"Destination "<<statuses[i].MPI_SOURCE<<" error code: "<<statuses[i].MPI_ERROR<<std::endl;
    }
 
    delete[] statuses;
    delete[] requests;
 
    for(std::size_t i=0; i<noOldNeighbours; ++i)
      delete[] sendBuffers_[i];
 
    delete[] sendBuffers_;
    delete[] sendBufferSizes_;
 
    // No need for the iterator tuples any more
    iteratorsMap_.clear();
 
    indexSet_.endResize();
 
    delete[] oldNeighbours;
 
    repairLocalIndexPointers(globalMap_, remoteIndices_, indexSet_);
 
    oldMap_.clear();
    globalMap_.clear();
 
    // update the sequence number
    remoteIndices_.sourceSeqNo_ = remoteIndices_.destSeqNo_ = indexSet_.seqNo();
  }
 
  template<typename T>
  void IndicesSyncer<T>::packAndSend(int destination, char* buffer, std::size_t bufferSize, MPI_Request& request)
  {
    typedef typename ParallelIndexSet::const_iterator IndexIterator;
 
    IndexIterator iEnd       = indexSet_.end();
    int bpos       = 0;
    int published  = 0;
    int pairs      = 0;
 
    assert(checkReset());
 
    // Pack the number of indices we publish
    MPI_Pack(&(infoSend_[destination].publish), 1, MPI_INT, buffer, bufferSize, &bpos,
             remoteIndices_.communicator());
 
    for(IndexIterator index = indexSet_.begin(); index != iEnd; ++index) {
      // Search for corresponding remote indices in all iterator tuples
      typedef typename IteratorsMap::iterator Iterator;
      Iterator iteratorsEnd = iteratorsMap_.end();
 
      // advance all iterators to a position with global index >= index->global()
      for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators) {
        while(iterators->second.isNotAtEnd() &&
              iterators->second.globalIndexPair().first < index->global())
          ++(iterators->second);
        assert(!iterators->second.isNotAtEnd() || iterators->second.globalIndexPair().first >= index->global());
      }
 
      // Add all remote indices positioned at global which were already present before calling sync
      // to the message.
      // Count how many remote indices we will send
      int indices = 0;
      bool knownRemote = false; // Is the remote process supposed to know this index?
 
      for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators)
      {
        std::pair<GlobalIndex,Attribute> p;
        if (iterators->second.isNotAtEnd())
        {
          p = iterators->second.globalIndexPair();
        }
 
        if(iterators->second.isNotAtEnd() && iterators->second.isOld()
           && iterators->second.globalIndexPair().first == index->global()) {
          indices++;
          if(destination == iterators->first)
            knownRemote = true;
        }
      }
 
      if(!knownRemote)
        // We do not need to send any indices
        continue;
 
      Dune::dverb<<rank_<<": sending "<<indices<<" for index "<<index->global()<<" to "<<destination<<std::endl;
 
 
      // Pack the global index, the attribute and the number
      MPI_Pack(const_cast<GlobalIndex*>(&(index->global())), 1, MPITraits<GlobalIndex>::getType(), buffer, bufferSize, &bpos,
               remoteIndices_.communicator());
 
      char attr = index->local().attribute();
      MPI_Pack(&attr, 1, MPI_CHAR, buffer, bufferSize, &bpos,
               remoteIndices_.communicator());
 
      // Pack the number of remote indices we send.
      MPI_Pack(&indices, 1, MPI_INT, buffer, bufferSize, &bpos,
               remoteIndices_.communicator());
 
      // Pack the information about the remote indices
      for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators)
        if(iterators->second.isNotAtEnd() && iterators->second.isOld()
           && iterators->second.globalIndexPair().first == index->global()) {
          int process = iterators->first;
 
          ++pairs;
          assert(pairs <= infoSend_[destination].pairs);
          MPI_Pack(&process, 1, MPI_INT, buffer, bufferSize, &bpos,
                   remoteIndices_.communicator());
          char attr2 = iterators->second.remoteIndex().attribute();
 
          MPI_Pack(&attr2, 1, MPI_CHAR, buffer, bufferSize, &bpos,
                   remoteIndices_.communicator());
          --indices;
        }
      assert(indices==0);
      ++published;
      Dune::dvverb<<" (publish="<<published<<", pairs="<<pairs<<")"<<std::endl;
      assert(published <= infoSend_[destination].publish);
    }
 
    // Make sure we send all expected entries
    assert(published == infoSend_[destination].publish);
    assert(pairs == infoSend_[destination].pairs);
    resetIteratorsMap();
 
    Dune::dverb << rank_<<": Sending message of "<<bpos<<" bytes to "<<destination<<std::endl;
 
    MPI_Issend(buffer, bpos, MPI_PACKED, destination, 345, remoteIndices_.communicator(),&request);
  }
 
  template<typename T>
  inline void IndicesSyncer<T>::insertIntoRemoteIndexList(int process,
                                                          const std::pair<GlobalIndex,Attribute>& globalPair,
                                                          char attribute)
  {
    Dune::dverb<<"Inserting from "<<process<<" "<<globalPair.first<<", "<<
    globalPair.second<<" "<<attribute<<std::endl;
 
    resetIteratorsMap();
 
    // There might be cases where there no remote indices for that process yet
    typename IteratorsMap::iterator found = iteratorsMap_.find(process);
 
    if( found == iteratorsMap_.end() ) {
      Dune::dverb<<"Discovered new neighbour "<<process<<std::endl;
      RemoteIndexList* rlist = new RemoteIndexList();
      remoteIndices_.remoteIndices_.insert(std::make_pair(process,std::make_pair(rlist,rlist)));
      Iterators iterators = Iterators(*rlist, globalMap_[process], oldMap_[process]);
      found = iteratorsMap_.insert(std::make_pair(process, iterators)).first;
    }
 
    Iterators& iterators = found->second;
 
    // Search for the remote index
    while(iterators.isNotAtEnd() && iterators.globalIndexPair() < globalPair) {
      // Increment all iterators
      ++iterators;
 
    }
 
    if(iterators.isAtEnd() || iterators.globalIndexPair() != globalPair) {
      // The entry is not yet known
      // Insert in the list and do not change the first iterator.
      iterators.insert(RemoteIndex(Attribute(attribute)),globalPair);
      return;
    }
 
    // Global indices match
    bool indexIsThere=false;
    for(Iterators tmpIterators = iterators;
        !tmpIterators.isAtEnd() && tmpIterators.globalIndexPair() == globalPair;
        ++tmpIterators)
      //entry already exists with the same attribute
      if(tmpIterators.globalIndexPair().second == attribute) {
        indexIsThere=true;
        break;
      }
 
    if(!indexIsThere)
      // The entry is not yet known
      // Insert in the list and do not change the first iterator.
      iterators.insert(RemoteIndex(Attribute(attribute)),globalPair);
  }
 
  template<typename T>
  template<typename T1>
  void IndicesSyncer<T>::recvAndUnpack(T1& numberer)
  {
    typedef typename ParallelIndexSet::const_iterator IndexIterator;
 
    IndexIterator iEnd   = indexSet_.end();
    IndexIterator index  = indexSet_.begin();
    int bpos   = 0;
    int publish;
 
    assert(checkReset());
 
    MPI_Status status;
 
    // We have to determine the message size and source before the receive
    MPI_Probe(MPI_ANY_SOURCE, 345, remoteIndices_.communicator(), &status);
 
    int source=status.MPI_SOURCE;
    int count;
    MPI_Get_count(&status, MPI_PACKED, &count);
 
    Dune::dvverb<<rank_<<": Receiving message from "<< source<<" with "<<count<<" bytes"<<std::endl;
 
    if(count>receiveBufferSize_) {
      receiveBufferSize_=count;
      delete[] receiveBuffer_;
      receiveBuffer_ = new char[receiveBufferSize_];
    }
 
    MPI_Recv(receiveBuffer_, count, MPI_PACKED, source, 345, remoteIndices_.communicator(), &status);
 
    // How many global entries were published?
    MPI_Unpack(receiveBuffer_,  count, &bpos, &publish, 1, MPI_INT, remoteIndices_.communicator());
 
    // Now unpack the remote indices and add them.
    while(publish>0) {
 
      // Unpack information about the local index on the source process
      GlobalIndex global;           // global index of the current entry
      char sourceAttribute; // Attribute on the source process
      int pairs;
 
      MPI_Unpack(receiveBuffer_,  count, &bpos, &global, 1, MPITraits<GlobalIndex>::getType(),
                 remoteIndices_.communicator());
      MPI_Unpack(receiveBuffer_,  count, &bpos, &sourceAttribute, 1, MPI_CHAR,
                 remoteIndices_.communicator());
      MPI_Unpack(receiveBuffer_,  count, &bpos, &pairs, 1, MPI_INT,
                 remoteIndices_.communicator());
 
      // Insert the entry on the remote process to our
      // remote index list
      SLList<std::pair<int,Attribute> > sourceAttributeList;
      sourceAttributeList.push_back(std::make_pair(source,Attribute(sourceAttribute)));
#ifndef NDEBUG
      bool foundSelf = false;
#endif
      Attribute myAttribute=Attribute();
 
      // Unpack the remote indices
      for(; pairs>0; --pairs) {
        // Unpack the process id that knows the index
        int process;
        char attribute;
        MPI_Unpack(receiveBuffer_,  count, &bpos, &process, 1, MPI_INT,
                   remoteIndices_.communicator());
        // Unpack the attribute
        MPI_Unpack(receiveBuffer_,  count, &bpos, &attribute, 1, MPI_CHAR,
                   remoteIndices_.communicator());
 
        if(process==rank_) {
#ifndef NDEBUG
          foundSelf=true;
#endif
          myAttribute=Attribute(attribute);
          // Now we know the local attribute of the global index
          //Only add the index if it is unknown.
          // Do we know that global index already?
          IndexIterator pos = std::lower_bound(index, iEnd, IndexPair(global));
 
          if(pos == iEnd || pos->global() != global) {
            // no entry with this global index
            indexSet_.add(global,
                          ParallelLocalIndex<Attribute>(numberer(global),
                                                        myAttribute, true));
            Dune::dvverb << "Adding "<<global<<" "<<myAttribute<<std::endl;
            continue;
          }
 
          // because of above the global indices match. Add only if the attribute is different
          bool indexIsThere = false;
          index=pos;
 
          for(; pos->global()==global; ++pos)
            if(pos->local().attribute() == myAttribute) {
              Dune::dvverb<<"found "<<global<<" "<<myAttribute<<std::endl;
              indexIsThere = true;
              break;
            }
 
          if(!indexIsThere) {
            indexSet_.add(global,
                          ParallelLocalIndex<Attribute>(numberer(global),
                                                        myAttribute, true));
            Dune::dvverb << "Adding "<<global<<" "<<myAttribute<<std::endl;
          }
 
        }else{
          sourceAttributeList.push_back(std::make_pair(process,Attribute(attribute)));
        }
      }
      assert(foundSelf);
      // Insert remote indices
      typedef typename SLList<std::pair<int,Attribute> >::const_iterator Iter;
      for(Iter i=sourceAttributeList.begin(), end=sourceAttributeList.end();
          i!=end; ++i)
        insertIntoRemoteIndexList(i->first, std::make_pair(global, myAttribute),
                                  i->second);
      --publish;
    }
 
    resetIteratorsMap();
  }
 
  template<typename T>
  void IndicesSyncer<T>::resetIteratorsMap(){
 
    // Reset iterators in all tuples.
    typedef typename IteratorsMap::iterator Iterator;
    typedef typename RemoteIndices::RemoteIndexMap::iterator
    RemoteIterator;
    typedef typename GlobalIndicesMap::iterator GlobalIterator;
    typedef typename BoolMap::iterator BoolIterator;
 
    const RemoteIterator remoteEnd = remoteIndices_.remoteIndices_.end();
    Iterator iterators = iteratorsMap_.begin();
    GlobalIterator global = globalMap_.begin();
    BoolIterator added = oldMap_.begin();
 
    for(RemoteIterator remote = remoteIndices_.remoteIndices_.begin();
        remote != remoteEnd; ++remote, ++global, ++added, ++iterators) {
      iterators->second.reset(*(remote->second.first), global->second, added->second);
    }
  }
 
  template<typename T>
  bool IndicesSyncer<T>::checkReset(const Iterators& iterators, RemoteIndexList& rList, GlobalIndexList& gList,
                                    BoolList& bList){
 
    if(get<0>(iterators.iterators_)!=rList.begin())
      return false;
    if(get<1>(iterators.iterators_)!=gList.begin())
      return false;
    if(get<2>(iterators.iterators_)!=bList.begin())
      return false;
    return true;
  }
 
 
  template<typename T>
  bool IndicesSyncer<T>::checkReset(){
 
    // Reset iterators in all tuples.
    typedef typename IteratorsMap::iterator Iterator;
    typedef typename RemoteIndices::RemoteIndexMap::iterator
    RemoteIterator;
    typedef typename GlobalIndicesMap::iterator GlobalIterator;
    typedef typename BoolMap::iterator BoolIterator;
 
    const RemoteIterator remoteEnd = remoteIndices_.remoteIndices_.end();
    Iterator iterators = iteratorsMap_.begin();
    GlobalIterator global = globalMap_.begin();
    BoolIterator added = oldMap_.begin();
    bool ret = true;
 
    for(RemoteIterator remote = remoteIndices_.remoteIndices_.begin();
        remote != remoteEnd; ++remote, ++global, ++added, ++iterators) {
      if(!checkReset(iterators->second, *(remote->second.first), global->second,
                     added->second))
        ret=false;
    }
    return ret;
  }
}
 
#endif
#endif