libdune-common-dev 2.2.1-2 / usr / include / dune / common / parallel / indicessyncer.hh

// $Id$
#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<cassert>
#include<cmath>
#include<limits>
#include<algorithm>
#include<functional>
#include<map>

#if HAVE_MPI
namespace Dune
{
  /** @addtogroup Common_Parallel
   *
   * @{
   */
   /**
   * @file
   * @brief Class for adding missing indices of a distributed index set in a local 
   * communication.
   * @author Markus Blatt
   */

  /**
   * @brief Class for recomputing missing indices of a distributed index set.
   *
   * Missing local and remote indices will be added.
   */
  template<typename T>
  class IndicesSyncer
  {
  public:

    /** @brief The type of the index set. */
    typedef T ParallelIndexSet;

    /** @brief The type of the index pair */
    typedef typename ParallelIndexSet::IndexPair IndexPair;
    
    /** @brief Type of the global index used in the index set. */
    typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
    
    /** @brief Type of the attribute used in the index set. */
    typedef typename ParallelIndexSet::LocalIndex::Attribute Attribute;

    /**
     * @brief Type of the remote indices. 
     */
    typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
    
    /**
     * @brief Constructor.
     *
     * The source as well as the target index set of the remote
     * indices have to be the same as the provided index set.
     * @param indexSet The index set with the information
     * of the locally present indices.
     * @param remoteIndices The remoteIndices.
     */
    IndicesSyncer(ParallelIndexSet& indexSet, 
		  RemoteIndices& remoteIndices);
    
    /**
     * @brief Sync the index set.
     *
     * Computes the missing indices in the local and the remote index list and adds them.
     * No global communication is necessary!
     * All indices added to the index will become the local index 
     * std::numeric_limits<size_t>::max()
     * 
     */
    void sync();

    /**
     * @brief Synce the index set and assign local numbers to new indices
     *
     * Computes the missing indices in the local and the remote index list and adds them.
     * No global communication is necessary!
     * @param numberer Functor providing the local indices for the added global indices.
     * has to provide a function size_t operator()(const TG& global) that provides the 
     * local index to a global one. It will be called for ascending global indices.
     * 
     */
    template<typename T1>
    void sync(T1& numberer);
    
  private:    
    
    /** @brief The set of locally present indices.*/
    ParallelIndexSet& indexSet_;

    /** @brief The remote indices. */
    RemoteIndices& remoteIndices_;
    
    /** @brief The send buffers for the neighbour processes. */
    char** sendBuffers_;

    /** @brief The receive buffer. */
    char* receiveBuffer_;
    
    /** @brief The size of the send buffers. */
    std::size_t* sendBufferSizes_;
    
    /** @brief The size of the receive buffer in bytes. */
    int receiveBufferSize_; // int because of MPI
    
    /**
     * @brief Information about the messages to send to a neighbouring process.
     */
    struct MessageInformation
    {
      MessageInformation()
	: publish(), pairs()
      {}
      /** @brief The number of indices we publish for the other process. */
      int publish;
      /**
       * @brief The number of pairs (attribute and process number)
       * we publish to the neighbour process.
       */
      int pairs;
    };

    /**
     * @brief Default numberer for sync().
     */
    class DefaultNumberer
    {
    public:
      /**
       * @brief Provide the lcoal index, always
       * std::numeric_limits<size_t>::max()
       * @param global The global index (ignored).
       */
      std::size_t operator()(const GlobalIndex& global)
      {
	return std::numeric_limits<size_t>::max();
      }
    };
    
    /** @brief The mpi datatype for the MessageInformation */
    MPI_Datatype datatype_;

    /** @brief Our rank. */
    int rank_;
    
    /** 
     * @brief List type for temporarily storing the global indices of the
     * remote indices.
     */
    typedef SLList<std::pair<GlobalIndex,Attribute>, typename RemoteIndices::Allocator> GlobalIndexList;
    
    /** @brief The modifying iterator for the global index list. */
    typedef typename GlobalIndexList::ModifyIterator GlobalIndexModifier;
    
    /** 
     * @brief The type of the iterator of GlobalIndexList
     */
    typedef typename SLList<GlobalIndex, typename RemoteIndices::Allocator>::iterator 
    GlobalIndexIterator;

    /** @brief Type of the map of ranks onto GlobalIndexLists. */
    typedef std::map<int, GlobalIndexList> GlobalIndicesMap;
    
    /** 
     * @brief Map of global index lists onto process ranks.
     *
     * As the pointers in the remote index lists become invalid due to
     * resorting the index set entries one has store the corresponding 
     * global index for each remote index. Thus the pointers can be adjusted
     * properly as a last step.
     */
    GlobalIndicesMap globalMap_;

    /**
     * @brief The type of the single linked list of bools.
     */
    typedef SLList<bool, typename RemoteIndices::Allocator> BoolList;

    /**
     * @brief The mutable iterator of the single linked bool list.
     */
    typedef typename BoolList::iterator BoolIterator;

    /** @brief The type of the modifying iterator for the list of bools. */
    typedef typename BoolList::ModifyIterator BoolListModifier;
    
    /** @brief The type of the map of bool lists. */
    typedef std::map<int,BoolList> BoolMap;

    /**
     * @brief Map of lists of bool indicating whether the remote index was present before
     * call of sync.
     */
    BoolMap oldMap_;
    
    /** @brief Information about the messages we send. */
    std::map<int,MessageInformation> infoSend_;
    
    /** @brief The type of the remote index list. */
    typedef typename RemoteIndices::RemoteIndexList RemoteIndexList;

    /** @brief The tyoe of the modifying iterator of the remote index list. */
    typedef typename RemoteIndexList::ModifyIterator RemoteIndexModifier;
    
    /** @brief The type of the remote inde. */
    typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
    
    /** @brief The iterator of the remote index list. */
    typedef typename RemoteIndexList::iterator RemoteIndexIterator;

    /** @brief The const iterator of the remote index list. */
    typedef typename RemoteIndexList::const_iterator ConstRemoteIndexIterator;

    /** @brief Type of the tuple of iterators needed for the adding of indices. */
    typedef Dune::tuple<RemoteIndexModifier,GlobalIndexModifier,BoolListModifier,
		  const ConstRemoteIndexIterator> IteratorTuple;

    /**
     * @brief A tuple of iterators.
     *
     * Insertion into a single linked list is only possible at the position after the one of the iterator.
     * Therefore for each linked list two iterators are needed: One position before the actual entry 
     * (for insertion) and one positioned at the actual position (for searching).
     */
    class Iterators
    {
      friend class IndicesSyncer<T>;
    public:
      /**
       * @brief Constructor.
       * 
       * Initializes all iterator to first entry and the one before the first entry, respectively.
       * @param remoteIndices The list of the remote indices.
       * @param globalIndices The list of the coresponding global indices. This is needed because the
       * the pointers to the local index will become invalid due to the merging of the index sets.
       * @param booleans Whether the remote index was there before the sync process started.
       */
      Iterators(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
		BoolList& booleans);
      
      /**
       * @brief Default constructor.
       */
      Iterators();
      
       /**
       * @brief Increment all iteraors.
       */
      Iterators& operator++();
      
      /**
       * @brief Insert a new remote index to the underlying remote index list.
       * @param index The remote index.
       * @param global The global index corresponding to the remote index.
       */
      void insert(const RemoteIndex& index, 
                  const std::pair<GlobalIndex,Attribute>& global);
      
      /**
       * @brief Get the remote index at current position.
       * @return The current remote index.
       */
      RemoteIndex& remoteIndex() const;
      
      /**
       * @brief Get the global index of the remote index at current position.
       * @return The current global index.
       */
      std::pair<GlobalIndex,Attribute>& globalIndexPair() const;

      Attribute& attribute() const;
      
      /**
       * @brief Was this entry already in the remote index list before the sync process?
       * @return True if the current index wasalready in the remote index list 
       * before the sync process.
       */
      bool isOld() const;
      
      /**
       * @brief Reset all the underlying iterators.
       *
       * Position them to first list entry and the entry before the first entry respectively.
       * @param remoteIndices The list of the remote indices.
       * @param globalIndices The list of the coresponding global indices. This is needed because the
       * the pointers to the local index will become invalid due to the merging of the index sets.
       * @param booleans Whether the remote index was there before the sync process started.
       */
      void reset(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
		 BoolList& booleans);

      /**
       * @brief Are we not at the end of the list?
       * @return True if the iterators are not positioned at the end of the list
       * and the tail of the list respectively.
       */
      bool isNotAtEnd() const;
      
      /**
       * @brief Are we at the end of the list?
       * @return True if the iterators are positioned at the end of the list
       * and the tail of the list respectively.
       */
      bool isAtEnd() const;
      
    private:
      /**
       * @brief The iterator tuple.
       *
       * The tuple consists of one iterator over a single linked list of remote indices
       * initially positioned before the first entry, one over a sll of global indices
       * , one over a all of bool values both postioned at the same entry. The another three
       * iterators of the same type positioned at the first entry. Last an iterator over the 
       * sll of remote indices positioned at the end.
       */
      IteratorTuple iterators_;
    };
    
    /** @brief Type of the map from ranks to iterator tuples. */
    typedef std::map<int,Iterators> IteratorsMap;
    
    /**
     * @brief The iterator tuples mapped on the neighbours.
     *
     * The key of the map is the rank of the neighbour.
     * The first entry in the tuple is an iterator over the remote indices 
     * initially positioned before the first entry. The second entry is an
     * iterator over the corresponding global indices also initially positioned 
     * before the first entry. The third entry an iterator over remote indices
     * initially positioned at the beginning. The last entry is the iterator over
     * the remote indices positioned at the end.
     */
    IteratorsMap iteratorsMap_;
        
    /** @brief Calculates the message sizes to send. */
    void calculateMessageSizes();

    /**
     * @brief Pack and send the message for another process. 
     * @param destination The rank of the process we send to.
     * @param buffer The allocated buffer to use.
     * @param bufferSize The size of the buffer.
     * @param req The MPI_Request to setup the nonblocking send.
     */
    void packAndSend(int destination, char* buffer, std::size_t bufferSize, MPI_Request& req);
    
    /** 
     * @brief Recv and unpack the message from another process and add the indices.
     * @param numberer Functor providing local indices for added global indices.
     */
    template<typename T1>
    void recvAndUnpack(T1& numberer);

    /**
     * @brief Register the MPI datatype for the MessageInformation.
     */
    void registerMessageDatatype();
    
    /**
     * @brief Insert an entry into the  remote index list if not yet present.
     */
    void insertIntoRemoteIndexList(int process, 
                                   const std::pair<GlobalIndex,Attribute>& global, 
                                   char attribute);

    /**
     * @brief Reset the iterator tuples of all neighbouring processes.
     */
    void resetIteratorsMap();
    
    /**
     * @brief Check whether the iterator tuples of all neighbouring processes
     * are reset.
     */
    bool checkReset();
    
    /**
     * @brief Check whether the iterator tuple is reset.
     *
     * @param iterators The iterator tuple to check.
     * @param rlist The SLList of the remote indices.
     * @param gList The SLList of the global indices.
     * @param bList The SLList of the bool values.
     */
    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);
  }

  /**
   * @brief Stores the corresponding global indices of the remote index information.
   *
   * Whenever a ParallelIndexSet is resized all RemoteIndices that use it will be invalided
   * as the pointers to the index set are invalid after calling ParallelIndexSet::Resize()
   * One can rebuild them by storing the global indices in a map with this function and later 
   * repairing the pointers by calling repairLocalIndexPointers.
   * 
   * @warning The RemoteIndices class has to be build with the same index set for both the 
   * sending and receiving side 
   * @param globalMap Map to store the corresponding global indices in.
   * @param remoteIndices The remote index information we need to store the corresponding global
   * indices of.
   * @param indexSet The index set that is for both the sending and receiving side of the remote
   * index information.
   */
  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,
        			       const T& indexSet)
  {
    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()));
      }
    }
  }
  
  /**
   * @brief Repair the pointers to the local indices in the remote indices.
   *
   * @param globalMap The map of the process number to the list of global indices
   * corresponding to the remote index list of the process.
   * @param remoteIndices The known remote indices.
   * @param indexSet The set of local indices of the current process.
   */
  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_Address( &(message.publish), displacement);
    MPI_Address( &(message.pairs), displacement+1);
    
    // Make the displacement relative
    MPI_Address(&message, &base);
    displacement[0] -= base;
    displacement[1] -= base;
    
    MPI_Type_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 i=0;
    
    for(RemoteIterator remote = remoteIndices_.begin(); remote != end; ++remote, ++i){
      typedef typename RemoteIndices::RemoteIndexList::const_iterator
	RemoteIndexIterator;

      oldNeighbours[i]=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(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(i = 0; i<noOldNeighbours; ++i)
	packAndSend(oldNeighbours[i], sendBuffers_[i], sendBufferSizes_[i], requests[i]);

    // Probe for incoming messages, receive and unpack them
    for(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(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 attr = iterators->second.remoteIndex().attribute();
	  
	  MPI_Pack(&attr, 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 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 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