dune-common  2.2.0
remoteindices.hh
Go to the documentation of this file.
00001 // $Id$
00002 #ifndef DUNE_REMOTEINDICES_HH
00003 #define DUNE_REMOTEINDICES_HH
00004 
00005 #include"indexset.hh"
00006 #include <dune/common/exceptions.hh>
00007 #include"plocalindex.hh"
00008 #include<dune/common/poolallocator.hh>
00009 #include<dune/common/sllist.hh>
00010 #include<dune/common/static_assert.hh>
00011 #include <dune/common/stdstreams.hh>
00012 #include<map>
00013 #include<set>
00014 #include<utility>
00015 #include<iostream>
00016 #include<algorithm>
00017 #include<iterator>
00018 #if HAVE_MPI
00019 #include <dune/common/mpitraits.hh>
00020 #include"mpi.h"
00021 
00022 namespace Dune{
00033 
00034   template<typename TG, typename TA>
00035   class MPITraits<IndexPair<TG,ParallelLocalIndex<TA> > >
00036   {
00037   public:
00038     inline static MPI_Datatype getType();
00039   private:
00040     static MPI_Datatype type;
00041   };
00042   
00043   
00044   template<typename T, typename A>
00045   class RemoteIndices;
00046   
00047   template<typename T1, typename T2>
00048   class RemoteIndex;
00049   
00050   template<typename T>
00051   class IndicesSyncer;
00052   
00053   template<typename T1, typename T2>
00054   std::ostream& operator<<(std::ostream& os, const RemoteIndex<T1,T2>& index);
00055   
00056 
00057   template<typename T, typename A, bool mode>
00058   class RemoteIndexListModifier;
00059 
00060  
00064   template<typename T1, typename T2>
00065   class RemoteIndex
00066   {
00067     template<typename T>
00068     friend class IndicesSyncer;
00069 
00070     template<typename T, typename A, typename A1>
00071     friend void repairLocalIndexPointers(std::map<int,SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> >&, 
00072                                          RemoteIndices<T,A1>&,
00073                                          const T&);
00074 
00075     template<typename T, typename A, bool mode>
00076     friend class RemoteIndexListModifier;
00077 
00078   public:    
00083     typedef T1 GlobalIndex;
00092     typedef T2 Attribute;
00093 
00097     typedef IndexPair<GlobalIndex,ParallelLocalIndex<Attribute> >
00098     PairType;
00099     
00104     const Attribute attribute() const;
00105 
00111     const PairType& localIndexPair() const;
00112 
00116     RemoteIndex();
00117     
00118     
00124     RemoteIndex(const T2& attribute, 
00125                 const PairType* local);
00126 
00127     
00133     RemoteIndex(const T2& attribute);
00134 
00135     bool operator==(const RemoteIndex& ri) const;
00136     
00137     bool operator!=(const RemoteIndex& ri) const;
00138   private:
00140     const PairType* localIndex_;
00141 
00143     char attribute_;
00144   };
00145   
00146   template<class T, class A>
00147   std::ostream& operator<<(std::ostream& os, const RemoteIndices<T,A>& indices);
00148   
00149   class InterfaceBuilder;
00150 
00151   template<class T, class A>
00152   class CollectiveIterator;
00153 
00154   template<class T>
00155     class IndicesSyncer;
00156   
00157   // forward declaration needed for friend declaration.
00158   template<typename T1, typename T2>
00159   class OwnerOverlapCopyCommunication;
00160   
00161 
00178   template<class T, class A=std::allocator<RemoteIndex<typename T::GlobalIndex,
00179                                                        typename T::LocalIndex::Attribute> > >
00180   class RemoteIndices
00181   {
00182     friend class InterfaceBuilder;
00183     friend class IndicesSyncer<T>;
00184     template<typename T1, typename A2, typename A1>
00185     friend void repairLocalIndexPointers(std::map<int,SLList<std::pair<typename T1::GlobalIndex, typename T1::LocalIndex::Attribute>,A2> >&, 
00186                                          RemoteIndices<T1,A1>&,
00187                                          const T1&);
00188     
00189     template<class G, class T1, class T2>
00190     friend void fillIndexSetHoles(const G& graph, Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm);
00191     friend std::ostream& operator<<<>(std::ostream&, const RemoteIndices<T>&);
00192     
00193   public:
00194    
00198     typedef T ParallelIndexSet;
00199  
00202     typedef CollectiveIterator<T,A> CollectiveIteratorT;
00203     
00207     typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00208     
00209     
00213     typedef typename ParallelIndexSet::LocalIndex LocalIndex;
00214 
00218     typedef typename LocalIndex::Attribute Attribute;
00219  
00223     typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
00224 
00225     
00229     typedef typename A::template rebind<RemoteIndex>::other Allocator;
00230 
00232     typedef Dune::SLList<RemoteIndex,Allocator>
00233     RemoteIndexList;
00234     
00236     typedef std::map<int, std::pair<RemoteIndexList*,RemoteIndexList*> >
00237     RemoteIndexMap;    
00238     
00239     typedef typename RemoteIndexMap::const_iterator const_iterator;
00240     
00258     inline RemoteIndices(const ParallelIndexSet& source, const ParallelIndexSet& destination, 
00259                          const MPI_Comm& comm, const std::vector<int>& neighbours=std::vector<int>(), bool includeSelf=false);
00260 
00261     RemoteIndices();
00262     
00270     void setIncludeSelf(bool includeSelf);
00271     
00288     void setIndexSets(const ParallelIndexSet& source, const ParallelIndexSet& destination, 
00289                       const MPI_Comm& comm, const std::vector<int>& neighbours=std::vector<int>());
00290 
00291     template<typename C>
00292     void setNeighbours(const C& neighbours)
00293     {
00294       typedef typename C::const_iterator Iter;
00295       neighbourIds.clear();
00296       neighbourIds.insert(neighbours.begin(), neighbours.end());
00297         
00298     }
00299     
00300     const std::set<int>& getNeighbours() const
00301     {
00302       return neighbourIds;
00303     }
00304     
00308     ~RemoteIndices();
00309     
00319     template<bool ignorePublic>
00320     void rebuild();
00321 
00322     bool operator==(const RemoteIndices& ri);
00323     
00331     inline bool isSynced() const;
00332 
00336     inline MPI_Comm communicator() const;
00337 
00352     template<bool mode, bool send>
00353     inline RemoteIndexListModifier<T,A,mode> getModifier(int process);
00354 
00361     inline const_iterator find(int proc) const;
00362 
00367     inline const_iterator begin() const;
00368     
00373     inline const_iterator end() const;
00374     
00378     template<bool send>
00379     inline CollectiveIteratorT iterator() const;
00380 
00384     inline void free();
00385     
00390     inline int neighbours() const;
00391 
00393     inline const ParallelIndexSet& sourceIndexSet() const;
00394     
00396     inline const ParallelIndexSet& destinationIndexSet() const;
00397 
00398   private:
00400     RemoteIndices(const RemoteIndices&)
00401     {} 
00402 
00404     const ParallelIndexSet* source_;
00405 
00407     const ParallelIndexSet* target_;
00408 
00410     MPI_Comm comm_;
00411 
00414     std::set<int> neighbourIds;
00415     
00417     const static int commTag_=333;
00418     
00423     int sourceSeqNo_;
00424     
00429     int destSeqNo_;
00430 
00434     bool publicIgnored;
00435 
00439     bool firstBuild;
00440     
00441     /*
00442      * @brief If true, sending from indices of the processor to other 
00443      * indices on the same processor is enabled even if the same indexset is used 
00444      * on both the
00445      * sending and receiving side.
00446      */
00447     bool includeSelf;
00448     
00450     typedef IndexPair<GlobalIndex, LocalIndex> 
00451     PairType;
00452     
00459     RemoteIndexMap remoteIndices_;
00460     
00471     template<bool ignorePublic>
00472     inline void buildRemote(bool includeSelf);
00473 
00479     inline int noPublic(const ParallelIndexSet& indexSet);
00480     
00492     template<bool ignorePublic>
00493     inline void packEntries(PairType** myPairs, const ParallelIndexSet& indexSet,
00494                             char* p_out, MPI_Datatype type, int bufferSize, 
00495                             int* position, int n);
00496     
00510     inline void unpackIndices(RemoteIndexList& remote, int remoteEntries,
00511                               PairType** local, int localEntries, char* p_in, 
00512                               MPI_Datatype type, int* positon, int bufferSize,
00513                               bool fromOurself);
00514 
00515     inline void unpackIndices(RemoteIndexList& send, RemoteIndexList& receive,
00516                               int remoteEntries, PairType** localSource,
00517                               int localSourceEntries, PairType** localDest,
00518                               int localDestEntries, char* p_in,
00519                               MPI_Datatype type, int* position, int bufferSize);
00520 
00521     void unpackCreateRemote(char* p_in, PairType** sourcePairs, PairType** DestPairs, 
00522                             int remoteProc,  int sourcePublish, int destPublish, 
00523                             int bufferSize, bool sendTwo, bool fromOurSelf=false);
00524   };
00525   
00543   template<class T, class A, bool mode>
00544   class RemoteIndexListModifier
00545   {
00546 
00547     template<typename T1, typename A1>
00548     friend class RemoteIndices;
00549   
00550   public:
00551     class InvalidPosition : public RangeError
00552     {};
00553 
00554     enum {
00563       MODIFYINDEXSET=mode
00564     };
00565   
00569     typedef T ParallelIndexSet;
00570 
00574     typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00575 
00579     typedef typename ParallelIndexSet::LocalIndex LocalIndex;
00580     
00584     typedef typename LocalIndex::Attribute Attribute;
00585     
00589     typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
00590     
00594     typedef A Allocator;
00595 
00597     typedef Dune::SLList<RemoteIndex,Allocator>
00598     RemoteIndexList;
00599 
00603     typedef SLListModifyIterator<RemoteIndex,Allocator> ModifyIterator;
00604     
00608     typedef typename RemoteIndexList::const_iterator ConstIterator;
00609     
00623     void insert(const RemoteIndex& index) throw(InvalidPosition);
00624     
00625     
00640     void insert(const RemoteIndex& index, const GlobalIndex& global) throw(InvalidPosition);
00641 
00649     bool remove(const GlobalIndex& global) throw(InvalidPosition);
00650 
00663     void repairLocalIndexPointers() throw(InvalidIndexSetState);
00664 
00665 
00666     RemoteIndexListModifier(const RemoteIndexListModifier&);
00667     
00672     RemoteIndexListModifier()
00673       : glist_()
00674     {}
00675     
00677     ~RemoteIndexListModifier()
00678     {
00679       if(glist_)
00680         delete glist_;
00681     }
00682     
00683   private:
00684     
00690     RemoteIndexListModifier(const ParallelIndexSet& indexSet,
00691                             RemoteIndexList& rList);
00692 
00693     typedef SLList<GlobalIndex,Allocator> GlobalList;
00694     typedef typename GlobalList::ModifyIterator GlobalModifyIterator;
00695     RemoteIndices<ParallelIndexSet>* remoteIndices_;
00696     RemoteIndexList* rList_;
00697     const ParallelIndexSet* indexSet_;
00698     GlobalList* glist_;
00699     ModifyIterator iter_;
00700     GlobalModifyIterator giter_;
00701     ConstIterator end_;
00702     bool first_;
00703     GlobalIndex last_;
00704   };
00705   
00710   template<class T, class A>
00711   class CollectiveIterator
00712   {
00713     
00717     typedef T ParallelIndexSet;
00718 
00722     typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00723 
00727     typedef typename ParallelIndexSet::LocalIndex LocalIndex;
00728     
00732     typedef typename LocalIndex::Attribute Attribute;
00733 
00735     typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
00736     
00738     typedef typename A::template rebind<RemoteIndex>::other Allocator;
00739 
00741     typedef Dune::SLList<RemoteIndex,Allocator> RemoteIndexList;    
00742 
00744     typedef std::map<int,std::pair<typename RemoteIndexList::const_iterator,
00745                                    const typename RemoteIndexList::const_iterator> >
00746     Map;
00747 
00748   public:
00749     
00751     typedef std::map<int, std::pair<RemoteIndexList*,RemoteIndexList*> >
00752     RemoteIndexMap;
00753     
00759     inline CollectiveIterator(const RemoteIndexMap& map_, bool send);
00760     
00769     inline void advance(const GlobalIndex& global);
00770 
00780     inline void advance(const GlobalIndex& global, const Attribute& attribute);
00781 
00782     CollectiveIterator& operator++();
00783     
00787     inline bool empty();
00788     
00795     class iterator
00796     {
00797     public:
00798       typedef typename Map::iterator RealIterator;
00799       typedef typename Map::iterator ConstRealIterator;
00800       
00801       
00803       iterator(const RealIterator& iter, const ConstRealIterator& end, GlobalIndex& index)
00804         : iter_(iter), end_(end), index_(index), hasAttribute(false)
00805       {
00806         // Move to the first valid entry
00807         while(iter_!=end_ && iter_->second.first->localIndexPair().global()!=index_)
00808           ++iter_;
00809       }
00810       
00811       iterator(const RealIterator& iter, const ConstRealIterator& end, GlobalIndex index,
00812                Attribute attribute)
00813         : iter_(iter), end_(end), index_(index), attribute_(attribute), hasAttribute(true)
00814       {
00815         // Move to the first valid entry or the end
00816         while(iter_!=end_ && (iter_->second.first->localIndexPair().global()!=index_
00817                               || iter_->second.first->localIndexPair().local().attribute()!=attribute))
00818           ++iter_;
00819       }
00821       iterator(const iterator& other)
00822         : iter_(other.iter_), end_(other.end_), index_(other.index_)
00823       { }
00824          
00826       iterator& operator++()
00827       { 
00828         ++iter_;
00829         // If entry is not valid move on
00830         while(iter_!=end_ && (iter_->second.first->localIndexPair().global()!=index_ ||
00831               (hasAttribute &&
00832                iter_->second.first->localIndexPair().local().attribute()!=attribute_)))
00833           ++iter_;
00834         assert(iter_==end_ || 
00835                (iter_->second.first->localIndexPair().global()==index_));
00836         assert(iter_==end_ || !hasAttribute ||
00837                (iter_->second.first->localIndexPair().local().attribute()==attribute_));
00838         return *this;
00839       }
00840       
00842       const RemoteIndex& operator*()const
00843       {
00844         return *(iter_->second.first);
00845       }
00846       
00848       int process() const
00849       {
00850         return iter_->first;
00851       }
00852       
00854       const RemoteIndex* operator->()const
00855       {
00856         return iter_->second.first.operator->();
00857       }
00858       
00860       bool operator==(const iterator& other)
00861       {
00862         return other.iter_==iter_;
00863       }
00864 
00866       bool operator!=(const iterator& other)
00867       {
00868         return other.iter_!=iter_;
00869       }
00870 
00871     private:
00872       iterator();
00873       
00874       RealIterator iter_;
00875       RealIterator end_;
00876       GlobalIndex index_;
00877       Attribute attribute_;
00878       bool hasAttribute;
00879     };
00880     
00881     iterator begin();
00882     
00883     iterator end();
00884     
00885   private:
00886     
00887     Map map_;
00888     GlobalIndex index_;
00889     Attribute attribute_;
00890     bool noattribute;
00891   };
00892     
00893   template<typename TG, typename TA>
00894   MPI_Datatype MPITraits<IndexPair<TG,ParallelLocalIndex<TA> > >::getType()
00895   {
00896     if(type==MPI_DATATYPE_NULL){
00897       int length[4];
00898       MPI_Aint disp[4];
00899       MPI_Datatype types[4] = {MPI_LB, MPITraits<TG>::getType(), 
00900                                MPITraits<ParallelLocalIndex<TA> >::getType(), MPI_UB};
00901       IndexPair<TG,ParallelLocalIndex<TA> > rep[2];
00902       length[0]=length[1]=length[2]=length[3]=1;
00903       MPI_Address(rep, disp); // lower bound of the datatype
00904       MPI_Address(&(rep[0].global_), disp+1);
00905       MPI_Address(&(rep[0].local_), disp+2);
00906       MPI_Address(rep+1, disp+3); // upper bound of the datatype
00907       for(int i=3; i >= 0; --i)
00908         disp[i] -= disp[0];
00909       MPI_Type_struct(4, length, disp, types, &type);
00910       MPI_Type_commit(&type);
00911     }
00912     return type;
00913   }
00914 
00915   template<typename TG, typename TA>
00916   MPI_Datatype MPITraits<IndexPair<TG,ParallelLocalIndex<TA> > >::type=MPI_DATATYPE_NULL;
00917 
00918   template<typename T1, typename T2>
00919   RemoteIndex<T1,T2>::RemoteIndex(const T2& attribute, const PairType* local)
00920     : localIndex_(local), attribute_(attribute)
00921   {}
00922 
00923   template<typename T1, typename T2>
00924   RemoteIndex<T1,T2>::RemoteIndex(const T2& attribute)
00925     : localIndex_(0), attribute_(attribute)
00926   {}
00927 
00928   template<typename T1, typename T2>
00929   RemoteIndex<T1,T2>::RemoteIndex()
00930     : localIndex_(0), attribute_()
00931   {}
00932   template<typename T1, typename T2>
00933   inline bool RemoteIndex<T1,T2>::operator==(const RemoteIndex& ri) const
00934   {
00935     return localIndex_==ri.localIndex_ && attribute_==ri.attribute;
00936   }
00937   
00938   template<typename T1, typename T2>
00939   inline bool RemoteIndex<T1,T2>::operator!=(const RemoteIndex& ri) const
00940   {
00941     return localIndex_!=ri.localIndex_ || attribute_!=ri.attribute_;
00942   }
00943   
00944   template<typename T1, typename T2>
00945   inline const T2 RemoteIndex<T1,T2>::attribute() const
00946   {
00947     return T2(attribute_);
00948   }
00949 
00950   template<typename T1, typename T2>
00951   inline const IndexPair<T1,ParallelLocalIndex<T2> >& RemoteIndex<T1,T2>::localIndexPair() const
00952   {
00953     return *localIndex_;
00954   }
00955   
00956   template<typename T, typename A>
00957   inline RemoteIndices<T,A>::RemoteIndices(const ParallelIndexSet& source,
00958                                          const ParallelIndexSet& destination,
00959                                          const MPI_Comm& comm,
00960                                            const std::vector<int>& neighbours,
00961                                            bool includeSelf_)
00962     : source_(&source), target_(&destination), comm_(comm),
00963       sourceSeqNo_(-1), destSeqNo_(-1), publicIgnored(false), firstBuild(true),
00964       includeSelf(includeSelf_)
00965   {
00966     setNeighbours(neighbours);
00967   }
00968 
00969   template<typename T, typename A>
00970   void RemoteIndices<T,A>::setIncludeSelf(bool b)
00971   {
00972     includeSelf=b;
00973   }
00974   
00975   template<typename T, typename A>
00976   RemoteIndices<T,A>::RemoteIndices()
00977     :source_(0), target_(0), sourceSeqNo_(-1), 
00978      destSeqNo_(-1), publicIgnored(false), firstBuild(true)
00979   {}
00980   
00981   template<class T, typename A>
00982   void RemoteIndices<T,A>::setIndexSets(const ParallelIndexSet& source,
00983                                       const ParallelIndexSet& destination,
00984                                       const MPI_Comm& comm,
00985                                       const std::vector<int>& neighbours)
00986   {
00987     free();
00988     source_ = &source;
00989     target_ = &destination;
00990     comm_ = comm;
00991     firstBuild = true;
00992     setNeighbours(neighbours);
00993   }
00994 
00995   template<typename T, typename A>
00996   const typename RemoteIndices<T,A>::ParallelIndexSet& 
00997   RemoteIndices<T,A>::sourceIndexSet() const
00998   {
00999     return *source_;
01000   }
01001   
01002     
01003   template<typename T, typename A>
01004   const typename RemoteIndices<T,A>::ParallelIndexSet& 
01005   RemoteIndices<T,A>::destinationIndexSet() const
01006   {
01007     return *target_;
01008   }
01009   
01010   
01011   template<typename T, typename A>
01012   RemoteIndices<T,A>::~RemoteIndices()
01013   {
01014     free();
01015   }
01016 
01017   template<typename T, typename A>
01018   template<bool ignorePublic>
01019   inline void RemoteIndices<T,A>::packEntries(IndexPair<GlobalIndex,LocalIndex>** pairs,
01020                                                 const ParallelIndexSet& indexSet,
01021                                                 char* p_out, MPI_Datatype type, 
01022                                                 int bufferSize,
01023                                                 int *position, int n)
01024   {
01025     // fill with own indices
01026     typedef typename ParallelIndexSet::const_iterator const_iterator;
01027     typedef IndexPair<GlobalIndex,LocalIndex> PairType;
01028     const const_iterator end = indexSet.end();
01029 
01030     //Now pack the source indices
01031     int i=0;
01032     for(const_iterator index = indexSet.begin(); index != end; ++index)
01033       if(ignorePublic || index->local().isPublic()){
01034         
01035         MPI_Pack(const_cast<PairType*>(&(*index)), 1, 
01036                  type,
01037                  p_out, bufferSize, position, comm_);
01038         pairs[i++] = const_cast<PairType*>(&(*index));
01039       
01040       }
01041     assert(i==n);
01042   }
01043   
01044   template<typename T, typename A>
01045   inline int RemoteIndices<T,A>::noPublic(const ParallelIndexSet& indexSet)
01046   {
01047     typedef typename ParallelIndexSet::const_iterator const_iterator;
01048     
01049     int noPublic=0;
01050     
01051     const const_iterator end=indexSet.end();
01052     for(const_iterator index=indexSet.begin(); index!=end; ++index)
01053       if(index->local().isPublic())
01054         noPublic++;
01055     
01056     return noPublic;
01057     
01058   }
01059     
01060 
01061   template<typename T, typename A>
01062   inline void RemoteIndices<T,A>::unpackCreateRemote(char* p_in, PairType** sourcePairs,
01063                                                      PairType** destPairs, int remoteProc, 
01064                                                      int sourcePublish, int destPublish,
01065                                                      int bufferSize, bool sendTwo,
01066                                                      bool fromOurSelf)
01067   {
01068     
01069       // unpack the number of indices we received
01070       int noRemoteSource=-1, noRemoteDest=-1;
01071       char twoIndexSets=0;
01072       int position=0;
01073       // Did we receive two index sets?
01074       MPI_Unpack(p_in, bufferSize, &position, &twoIndexSets, 1, MPI_CHAR, comm_);
01075       // The number of source indices received
01076       MPI_Unpack(p_in, bufferSize, &position, &noRemoteSource, 1, MPI_INT, comm_);
01077       // The number of destination indices received
01078       MPI_Unpack(p_in, bufferSize, &position, &noRemoteDest, 1, MPI_INT, comm_);      
01079 
01080 
01081       // Indices for which we receive
01082       RemoteIndexList* receive= new RemoteIndexList();
01083       // Indices for which we send
01084       RemoteIndexList* send=0;
01085 
01086       MPI_Datatype type= MPITraits<PairType>::getType();
01087       
01088       if(!twoIndexSets){
01089         if(sendTwo){
01090           send = new RemoteIndexList();
01091           // Create both remote index sets simultaneously
01092           unpackIndices(*send, *receive, noRemoteSource, sourcePairs, sourcePublish,
01093                         destPairs, destPublish, p_in, type, &position, bufferSize);
01094         }else{
01095           // we only need one list
01096           unpackIndices(*receive, noRemoteSource, sourcePairs, sourcePublish,
01097                         p_in, type, &position, bufferSize, fromOurSelf);
01098           send=receive;
01099         }
01100       }else{
01101         
01102         int oldPos=position;
01103         // Two index sets received
01104         unpackIndices(*receive, noRemoteSource, destPairs, destPublish,
01105                       p_in, type, &position, bufferSize, fromOurSelf);
01106         if(!sendTwo)
01107           //unpack source entries again as destination entries
01108           position=oldPos;
01109         
01110         send = new RemoteIndexList();
01111         unpackIndices(*send, noRemoteDest, sourcePairs, sourcePublish, 
01112                       p_in, type, &position, bufferSize, fromOurSelf);
01113       }
01114       
01115       if(receive->empty() && send->empty()){
01116         if(send==receive){
01117           delete send;
01118         }else{
01119           delete send;
01120           delete receive;
01121         }
01122       }else{
01123         remoteIndices_.insert(std::make_pair(remoteProc, 
01124                                              std::make_pair(send,receive)));
01125       }
01126   }
01127   
01128 
01129   template<typename T, typename A>
01130   template<bool ignorePublic>
01131   inline void RemoteIndices<T,A>::buildRemote(bool includeSelf)
01132   {
01133     // Processor configuration
01134     int rank, procs;
01135     MPI_Comm_rank(comm_, &rank);
01136     MPI_Comm_size(comm_, &procs);
01137     
01138     // number of local indices to publish
01139     // The indices of the destination will be send.
01140     int sourcePublish, destPublish;
01141     
01142     // Do we need to send two index sets?
01143     char sendTwo = (source_ != target_);
01144     
01145     if(procs==1 && !(sendTwo || includeSelf))
01146       // Nothing to communicate
01147       return;   
01148 
01149     sourcePublish = (ignorePublic)? source_->size() : noPublic(*source_);
01150     
01151     if(sendTwo)
01152       destPublish = (ignorePublic)? target_->size() : noPublic(*target_);
01153     else  
01154       // we only need to send one set of indices
01155       destPublish = 0;
01156         
01157     int maxPublish, publish=sourcePublish+destPublish;
01158 
01159     // Calucate maximum number of indices send
01160     MPI_Allreduce(&publish, &maxPublish, 1, MPI_INT, MPI_MAX, comm_);
01161 
01162     // allocate buffers
01163     typedef IndexPair<GlobalIndex,LocalIndex> PairType;
01164     
01165     PairType** destPairs;
01166     PairType** sourcePairs = new PairType*[sourcePublish>0?sourcePublish:1];
01167     
01168     if(sendTwo)
01169       destPairs = new PairType*[destPublish>0?destPublish:1];
01170     else
01171       destPairs=sourcePairs;
01172     
01173     char** buffer = new char*[2];
01174     int bufferSize;    
01175     int position=0;
01176     int intSize;
01177     int charSize;
01178 
01179     // calculate buffer size
01180     MPI_Datatype type = MPITraits<PairType>::getType();
01181     
01182     MPI_Pack_size(maxPublish, type, comm_,
01183                   &bufferSize);
01184     MPI_Pack_size(1, MPI_INT, comm_,
01185                   &intSize);
01186     MPI_Pack_size(1, MPI_CHAR, comm_,
01187                   &charSize);
01188     // Our message will contain the following:
01189     // a bool wether two index sets where sent
01190     // the size of the source and the dest indexset, 
01191     // then the source and destination indices
01192     bufferSize += 2 * intSize + charSize;
01193     
01194     if(bufferSize<=0) bufferSize=1;
01195     
01196     buffer[0] = new char[bufferSize];
01197     buffer[1] = new char[bufferSize];
01198     
01199     
01200     // pack entries into buffer[0], p_out below!
01201     MPI_Pack(&sendTwo, 1, MPI_CHAR, buffer[0], bufferSize, &position,
01202                  comm_);
01203         
01204     // The number of indices we send for each index set
01205     MPI_Pack(&sourcePublish, 1, MPI_INT, buffer[0], bufferSize, &position, 
01206              comm_);
01207     MPI_Pack(&destPublish, 1, MPI_INT, buffer[0], bufferSize, &position,
01208              comm_);
01209         
01210     // Now pack the source indices and setup the destination pairs
01211     packEntries<ignorePublic>(sourcePairs, *source_, buffer[0], type, 
01212                               bufferSize, &position, sourcePublish);
01213     // If necessary send the dest indices and setup the source pairs
01214     if(sendTwo)
01215       packEntries<ignorePublic>(destPairs, *target_, buffer[0], type,
01216                                 bufferSize, &position, destPublish);
01217 
01218 
01219     // Update remote indices for ourself
01220     if(sendTwo|| includeSelf)
01221       unpackCreateRemote(buffer[0], sourcePairs, destPairs, rank, sourcePublish, 
01222                          destPublish, bufferSize, sendTwo, includeSelf);
01223 
01224     neighbourIds.erase(rank);
01225 
01226     if(neighbourIds.size()==0)
01227       {
01228         Dune::dvverb<<rank<<": Sending messages in a ring"<<std::endl;
01229         // send messages in ring
01230         for(int proc=1; proc<procs; proc++){
01231           // pointers to the current input and output buffers
01232           char* p_out = buffer[1-(proc%2)];
01233           char* p_in = buffer[proc%2];
01234           
01235           MPI_Status status;
01236           if(rank%2==0){
01237             MPI_Ssend(p_out, bufferSize, MPI_PACKED, (rank+1)%procs,
01238                       commTag_, comm_);
01239             MPI_Recv(p_in, bufferSize, MPI_PACKED, (rank+procs-1)%procs,
01240                      commTag_, comm_, &status);
01241           }else{
01242             MPI_Recv(p_in, bufferSize, MPI_PACKED, (rank+procs-1)%procs,
01243                      commTag_, comm_, &status);
01244             MPI_Ssend(p_out, bufferSize, MPI_PACKED, (rank+1)%procs,
01245                       commTag_, comm_);
01246           }
01247 
01248 
01249           // The process these indices are from
01250           int remoteProc = (rank+procs-proc)%procs;
01251           
01252           unpackCreateRemote(p_in, sourcePairs, destPairs, remoteProc, sourcePublish, 
01253                              destPublish, bufferSize, sendTwo);
01254           
01255         }
01256         
01257       }
01258     else
01259       {
01260         MPI_Request* requests=new MPI_Request[neighbourIds.size()];
01261         MPI_Request* req=requests;
01262         
01263         typedef typename std::set<int>::size_type size_type;
01264         size_type noNeighbours=neighbourIds.size();
01265 
01266         // setup sends  
01267         for(std::set<int>::iterator neighbour=neighbourIds.begin();
01268             neighbour!= neighbourIds.end(); ++neighbour){
01269             // Only send the information to the neighbouring processors
01270             MPI_Issend(buffer[0], position , MPI_PACKED, *neighbour, commTag_, comm_, req++);
01271         }
01272         
01273         //Test for received messages
01274         
01275         for(size_type received=0; received <noNeighbours; ++received)
01276           {
01277             MPI_Status status;
01278             // probe for next message
01279             MPI_Probe(MPI_ANY_SOURCE, commTag_, comm_, &status);
01280             int remoteProc=status.MPI_SOURCE;
01281             int size;
01282             MPI_Get_count(&status, MPI_PACKED, &size);
01283             // receive message
01284             MPI_Recv(buffer[1], size, MPI_PACKED, remoteProc,
01285                      commTag_, comm_, &status);
01286               
01287             unpackCreateRemote(buffer[1], sourcePairs, destPairs, remoteProc, sourcePublish, 
01288                                destPublish, bufferSize, sendTwo);
01289           }
01290         // wait for completion of pending requests
01291         MPI_Status* statuses = new MPI_Status[neighbourIds.size()];
01292         
01293         if(MPI_ERR_IN_STATUS==MPI_Waitall(neighbourIds.size(), requests, statuses)){
01294           for(size_type i=0; i < neighbourIds.size(); ++i)
01295             if(statuses[i].MPI_ERROR!=MPI_SUCCESS){
01296               std::cerr<<rank<<": MPI_Error occurred while receiving message."<<std::endl;
01297               MPI_Abort(comm_, 999);
01298           }
01299         }
01300         delete[] requests;
01301         delete[] statuses;
01302       }
01303     
01304     
01305     // delete allocated memory
01306     if(destPairs!=sourcePairs)
01307       delete[] destPairs;
01308     
01309     delete[] sourcePairs;
01310     delete[] buffer[0];
01311     delete[] buffer[1];
01312     delete[] buffer;
01313   }
01314   
01315   template<typename T, typename A>
01316   inline void RemoteIndices<T,A>::unpackIndices(RemoteIndexList& remote,
01317                                                 int remoteEntries,
01318                                                 PairType** local,
01319                                                 int localEntries,
01320                                                 char* p_in,
01321                                                 MPI_Datatype type,
01322                                                 int* position,
01323                                                 int bufferSize,
01324                                                 bool fromOurSelf)
01325   {
01326     if(remoteEntries==0)
01327       return;
01328   
01329     PairType index(-1);
01330     MPI_Unpack(p_in, bufferSize, position, &index, 1, 
01331                type, comm_);
01332     GlobalIndex oldGlobal=index.global();
01333     int n_in=0, localIndex=0;
01334         
01335     //Check if we know the global index
01336     while(localIndex<localEntries){
01337       if(local[localIndex]->global()==index.global()){
01338           int oldLocalIndex=localIndex;
01339           
01340           while(localIndex<localEntries && 
01341                 local[localIndex]->global()==index.global()){
01342             if(!fromOurSelf || index.local().attribute() != 
01343                local[localIndex]->local().attribute())
01344               // if index is from us it has to have a different attribute
01345               remote.push_back(RemoteIndex(index.local().attribute(), 
01346                                            local[localIndex]));
01347             localIndex++;
01348           }
01349           
01350         // unpack next remote index
01351         if((++n_in) < remoteEntries){
01352           MPI_Unpack(p_in, bufferSize, position, &index, 1, 
01353                      type, comm_);
01354           if(index.global()==oldGlobal)
01355             // Restart comparison for the same global indices
01356             localIndex=oldLocalIndex;
01357           else
01358             oldGlobal=index.global();
01359         }else{
01360           // No more received indices
01361           break;
01362         }
01363         continue;
01364       }
01365       
01366       if (local[localIndex]->global()<index.global()){
01367         // compare with next entry in our list
01368         ++localIndex;
01369       }else{
01370         // We do not know the index, unpack next
01371         if((++n_in) < remoteEntries){
01372           MPI_Unpack(p_in, bufferSize, position, &index, 1, 
01373                      type, comm_);
01374           oldGlobal=index.global();
01375         }else
01376           // No more received indices
01377           break;
01378       }
01379     }
01380     
01381     // Unpack the other received indices without doing anything
01382     while(++n_in < remoteEntries)
01383       MPI_Unpack(p_in, bufferSize, position, &index, 1, 
01384                type, comm_);
01385   }
01386   
01387     
01388   template<typename T, typename A>
01389   inline void RemoteIndices<T,A>::unpackIndices(RemoteIndexList& send,
01390                                                   RemoteIndexList& receive,
01391                                                   int remoteEntries,
01392                                                   PairType** localSource,
01393                                                   int localSourceEntries,
01394                                                   PairType** localDest,
01395                                                   int localDestEntries,
01396                                                   char* p_in,
01397                                                   MPI_Datatype type,
01398                                                   int* position,
01399                                                   int bufferSize)
01400   {             
01401     int n_in=0, sourceIndex=0, destIndex=0;
01402         
01403     //Check if we know the global index
01404     while(n_in<remoteEntries && (sourceIndex<localSourceEntries || destIndex<localDestEntries)){
01405       // Unpack next index
01406       PairType index;
01407       MPI_Unpack(p_in, bufferSize, position, &index, 1, 
01408                  type, comm_);
01409       n_in++;
01410       
01411       // Advance until global index in localSource and localDest are >= than the one in the unpacked index
01412       while(sourceIndex<localSourceEntries && localSource[sourceIndex]->global()<index.global())
01413         sourceIndex++;
01414       
01415       while(destIndex<localDestEntries && localDest[destIndex]->global()<index.global())
01416         destIndex++;
01417 
01418       // Add a remote index if we found the global index.
01419       if(sourceIndex<localSourceEntries && localSource[sourceIndex]->global()==index.global())
01420           send.push_back(RemoteIndex(index.local().attribute(), 
01421                                               localSource[sourceIndex]));
01422 
01423       if(destIndex < localDestEntries && localDest[destIndex]->global() == index.global())
01424           receive.push_back(RemoteIndex(index.local().attribute(), 
01425                                         localDest[sourceIndex]));
01426     }
01427     
01428   }
01429     
01430   template<typename T, typename A>
01431   inline void RemoteIndices<T,A>::free()
01432   {
01433     typedef typename RemoteIndexMap::iterator Iterator;
01434     Iterator lend = remoteIndices_.end();
01435     for(Iterator lists=remoteIndices_.begin(); lists != lend; ++lists){
01436       if(lists->second.first==lists->second.second){
01437         // there is only one remote index list.
01438         delete lists->second.first;
01439       }else{
01440         delete lists->second.first;
01441         delete lists->second.second;
01442       }
01443     }
01444     remoteIndices_.clear();
01445     firstBuild=true;
01446   }
01447 
01448   template<typename T, typename A>
01449   inline int RemoteIndices<T,A>::neighbours() const
01450   {
01451     return remoteIndices_.size();
01452   }
01453   
01454   template<typename T, typename A>
01455   template<bool ignorePublic>
01456   inline void RemoteIndices<T,A>::rebuild()
01457   {
01458     // Test wether a rebuild is Needed.
01459     if(firstBuild || 
01460        ignorePublic!=publicIgnored || !
01461        isSynced()){
01462       free();
01463       
01464       buildRemote<ignorePublic>(includeSelf);
01465 
01466       sourceSeqNo_ = source_->seqNo();
01467       destSeqNo_ = target_->seqNo();
01468       firstBuild=false;
01469       publicIgnored=ignorePublic;
01470     }
01471     
01472     
01473   }
01474   
01475   template<typename T, typename A>
01476   inline bool RemoteIndices<T,A>::isSynced() const
01477   {
01478     return sourceSeqNo_==source_->seqNo() && destSeqNo_ ==target_->seqNo();
01479   }
01480 
01481   template<typename T, typename A>
01482   template<bool mode, bool send>
01483   RemoteIndexListModifier<T,A,mode> RemoteIndices<T,A>::getModifier(int process)
01484   {
01485 
01486     // The user are on their own now!
01487     // We assume they know what they are doing and just set the
01488     // remote indices to synced status.
01489     sourceSeqNo_ = source_->seqNo();
01490     destSeqNo_ = target_->seqNo();
01491 
01492     typename RemoteIndexMap::iterator found = remoteIndices_.find(process);
01493     
01494     if(found == remoteIndices_.end())
01495     {
01496       if(source_ != target_)
01497         remoteIndices_.insert(std::make_pair(process,
01498                                              std::make_pair(new RemoteIndexList(), 
01499                                                             new RemoteIndexList())));
01500       else{
01501         RemoteIndexList* rlist = new RemoteIndexList();
01502         remoteIndices_.insert(std::make_pair(process,
01503                                              std::make_pair(rlist, rlist)));
01504         
01505         found = remoteIndices_.find(process);
01506       }
01507     }
01508     
01509     firstBuild = false;
01510     
01511     if(send)
01512       return RemoteIndexListModifier<T,A,mode>(*source_, *(found->second.first));
01513     else
01514       return RemoteIndexListModifier<T,A,mode>(*target_, *(found->second.second));
01515   }
01516 
01517   template<typename T, typename A>
01518   inline typename RemoteIndices<T,A>::const_iterator
01519   RemoteIndices<T,A>::find(int proc) const
01520   {
01521     return remoteIndices_.find(proc);
01522   }
01523 
01524   template<typename T, typename A>
01525   inline typename RemoteIndices<T,A>::const_iterator
01526   RemoteIndices<T,A>::begin() const
01527   {
01528     return remoteIndices_.begin();
01529   }
01530   
01531   template<typename T, typename A>
01532   inline typename RemoteIndices<T,A>::const_iterator
01533   RemoteIndices<T,A>::end() const
01534   {
01535     return remoteIndices_.end();
01536   }
01537 
01538 
01539   template<typename T, typename A>
01540   bool RemoteIndices<T,A>::operator==(const RemoteIndices& ri)
01541   {
01542     if(neighbours()!=ri.neighbours())
01543       return false;
01544     
01545     typedef RemoteIndexList RList;
01546     typedef typename std::map<int,std::pair<RList*,RList*> >::const_iterator const_iterator;
01547     
01548     const const_iterator rend = remoteIndices_.end();
01549 
01550     for(const_iterator rindex = remoteIndices_.begin(), rindex1=ri.remoteIndices_.begin(); rindex!=rend; ++rindex, ++rindex1){
01551       if(rindex->first != rindex1->first)
01552         return false;
01553       if(*(rindex->second.first) != *(rindex1->second.first))
01554         return false;
01555       if(*(rindex->second.second) != *(rindex1->second.second))
01556         return false;
01557     }
01558     return true;
01559   }
01560   
01561   template<class T, class A, bool mode>
01562   RemoteIndexListModifier<T,A,mode>::RemoteIndexListModifier(const ParallelIndexSet& indexSet,
01563                                                            RemoteIndexList& rList)
01564     : rList_(&rList), indexSet_(&indexSet), glist_(new GlobalList()), iter_(rList.beginModify()), end_(rList.end()), first_(true)
01565   {
01566     if(MODIFYINDEXSET){
01567       assert(indexSet_);
01568       for(ConstIterator iter=iter_; iter != end_; ++iter)
01569         glist_->push_back(iter->localIndexPair().global());
01570       giter_ = glist_->beginModify();
01571     }
01572   }
01573 
01574   template<typename T, typename A, bool mode>
01575   RemoteIndexListModifier<T,A,mode>::RemoteIndexListModifier(const RemoteIndexListModifier<T,A,mode>& other)
01576     : rList_(other.rList_), indexSet_(other.indexSet_), 
01577       glist_(other.glist_), iter_(other.iter_), giter_(other.giter_), end_(other.end_), 
01578       first_(other.first_), last_(other.last_)
01579   {}
01580     
01581   template<typename T, typename A, bool mode>
01582   inline void RemoteIndexListModifier<T,A,mode>::repairLocalIndexPointers() throw(InvalidIndexSetState)
01583   { 
01584     if(MODIFYINDEXSET){
01585       // repair pointers to local index set.
01586 #ifdef DUNE_ISTL_WITH_CHECKING
01587       if(indexSet_->state()!=GROUND)
01588         DUNE_THROW(InvalidIndexSetState, "Index has to be in ground mode for repairing pointers to indices");
01589 #endif
01590       typedef typename ParallelIndexSet::const_iterator IndexIterator;
01591       typedef typename GlobalList::const_iterator GlobalIterator;
01592       typedef typename RemoteIndexList::iterator Iterator;
01593       GlobalIterator giter = glist_->begin();
01594       IndexIterator index = indexSet_->begin();
01595       
01596       for(Iterator iter=rList_->begin(); iter != end_; ++iter){
01597         while(index->global()<*giter){
01598           ++index;
01599 #ifdef DUNE_ISTL_WITH_CHECKING
01600           if(index == indexSet_.end())
01601             DUNE_THROW(InvalidPosition, "No such global index in set!");
01602 #endif
01603         }
01604 
01605 #ifdef DUNE_ISTL_WITH_CHECKING
01606           if(index->global() != *giter)
01607             DUNE_THROW(InvalidPosition, "No such global index in set!");
01608 #endif
01609           iter->localIndex_ = &(*index);
01610       }
01611     }
01612   }
01613   
01614   template<typename T, typename A, bool mode>
01615   inline void RemoteIndexListModifier<T,A,mode>::insert(const RemoteIndex& index) throw(InvalidPosition)
01616   {
01617     dune_static_assert(!mode,"Not allowed if the mode indicates that new indices"
01618                        "might be added to the underlying index set. Use "
01619                        "insert(const RemoteIndex&, const GlobalIndex&) instead");
01620     
01621 #ifdef DUNE_ISTL_WITH_CHECKING
01622     if(!first_ && index.localIndexPair().global()<last_)
01623       DUNE_THROW(InvalidPosition, "Modifcation of remote indices have to occur with ascending global index.");
01624 #endif
01625     // Move to the correct position
01626     while(iter_ != end_ && iter_->localIndexPair().global() < index.localIndexPair().global()){
01627       ++iter_;
01628     }
01629     
01630     // No duplicate entries allowed
01631     assert(iter_==end_ || iter_->localIndexPair().global() != index.localIndexPair().global());
01632     iter_.insert(index);
01633     last_ = index.localIndexPair().global();
01634     first_ = false;
01635   }
01636   
01637   template<typename T, typename A, bool mode>
01638   inline void RemoteIndexListModifier<T,A,mode>::insert(const RemoteIndex& index, const GlobalIndex& global) throw(InvalidPosition)
01639   {
01640      dune_static_assert(mode,"Not allowed if the mode indicates that no new indices"
01641                        "might be added to the underlying index set. Use "
01642                        "insert(const RemoteIndex&) instead");
01643 #ifdef DUNE_ISTL_WITH_CHECKING
01644     if(!first_ && global<last_)
01645       DUNE_THROW(InvalidPosition, "Modification of remote indices have to occur with ascending global index.");
01646 #endif
01647     // Move to the correct position
01648     while(iter_ != end_ && *giter_ < global){
01649       ++giter_;
01650       ++iter_;
01651     }
01652     
01653     // No duplicate entries allowed
01654     assert(iter_->localIndexPair().global() != global);
01655     iter_.insert(index);
01656     giter_.insert(global);
01657     
01658     last_ = global;
01659     first_ = false;
01660   }
01661   
01662   template<typename T, typename A, bool mode>
01663   bool RemoteIndexListModifier<T,A,mode>::remove(const GlobalIndex& global) throw(InvalidPosition)
01664   {
01665 #ifdef DUNE_ISTL_WITH_CHECKING
01666     if(!first_ && global<last_)
01667       DUNE_THROW(InvalidPosition, "Modifcation of remote indices have to occur with ascending global index.");
01668 #endif
01669 
01670     bool found= false;
01671     
01672     if(MODIFYINDEXSET){
01673     // Move to the correct position
01674       while(iter_!=end_ && *giter_< global){
01675         ++giter_;
01676         ++iter_;
01677       }
01678       if(*giter_ == global){
01679         giter_.remove();
01680         iter_.remove();
01681         found=true;
01682       }
01683     }else{
01684       while(iter_!=end_ && iter_->localIndexPair().global() < global)
01685         ++iter_;
01686         
01687       if(iter_->localIndexPair().global()==global){
01688         iter_.remove();
01689         found = true;
01690       }
01691     }
01692     
01693     last_ = global;
01694     first_ = false;
01695     return found;
01696   }
01697 
01698   template<typename T, typename A>
01699   template<bool send>
01700   inline typename RemoteIndices<T,A>::CollectiveIteratorT RemoteIndices<T,A>::iterator() const
01701   {
01702     return CollectiveIterator<T,A>(remoteIndices_, send);
01703   }
01704 
01705   template<typename T, typename A>
01706   inline MPI_Comm RemoteIndices<T,A>::communicator() const
01707   {
01708     return comm_;
01709     
01710   }
01711   
01712   template<typename T, typename A>
01713   CollectiveIterator<T,A>::CollectiveIterator(const RemoteIndexMap& pmap, bool send)
01714   {
01715     typedef typename RemoteIndexMap::const_iterator const_iterator;
01716     typedef typename RemoteIndexMap::iterator iterator;
01717     
01718     const const_iterator end=pmap.end();
01719     for(const_iterator process=pmap.begin(); process != end; ++process){
01720       const RemoteIndexList* list = send? process->second.first : process->second.second;
01721       typedef typename RemoteIndexList::const_iterator iterator;
01722       map_.insert(std::make_pair(process->first, 
01723                                  std::pair<iterator, const iterator>(list->begin(), list->end())));
01724     }
01725   }
01726 
01727   template<typename T, typename A>
01728   inline void CollectiveIterator<T,A>::advance(const GlobalIndex& index)
01729   {
01730     typedef typename Map::iterator iterator;
01731     typedef typename Map::const_iterator const_iterator;
01732     const const_iterator end = map_.end();
01733     
01734     for(iterator iter = map_.begin(); iter != end;){
01735       // Step the iterator until we are >= index
01736       typename RemoteIndexList::const_iterator current = iter->second.first;
01737       typename RemoteIndexList::const_iterator rend = iter->second.second;
01738       RemoteIndex remoteIndex;
01739       if(current != rend)
01740         remoteIndex = *current;
01741       
01742       while(iter->second.first!=iter->second.second && iter->second.first->localIndexPair().global()<index)
01743         ++(iter->second.first);
01744       
01745       // erase from the map if there are no more entries.
01746       if(iter->second.first == iter->second.second)
01747         map_.erase(iter++);
01748       else{
01749         ++iter;
01750       }
01751     }
01752     index_=index;
01753     noattribute=true;
01754   }
01755 
01756   template<typename T, typename A>
01757   inline void CollectiveIterator<T,A>::advance(const GlobalIndex& index,
01758                                                const Attribute& attribute)
01759   {
01760     typedef typename Map::iterator iterator;
01761     typedef typename Map::const_iterator const_iterator;
01762     const const_iterator end = map_.end();
01763     
01764     for(iterator iter = map_.begin(); iter != end;){
01765       // Step the iterator until we are >= index
01766       typename RemoteIndexList::const_iterator current = iter->second.first;
01767       typename RemoteIndexList::const_iterator rend = iter->second.second;
01768       RemoteIndex remoteIndex;
01769       if(current != rend)
01770         remoteIndex = *current;
01771 
01772       // Move to global index or bigger
01773       while(iter->second.first!=iter->second.second && iter->second.first->localIndexPair().global()<index)
01774         ++(iter->second.first);
01775       
01776       // move to attribute or bigger
01777       while(iter->second.first!=iter->second.second 
01778             && iter->second.first->localIndexPair().global()==index
01779             && iter->second.first->localIndexPair().local().attribute()<attribute)
01780         ++(iter->second.first);
01781 
01782       // erase from the map if there are no more entries.
01783       if(iter->second.first == iter->second.second)
01784         map_.erase(iter++);
01785       else{
01786         ++iter;
01787       }
01788     }
01789     index_=index;
01790     attribute_=attribute;
01791     noattribute=false;
01792   }
01793 
01794   template<typename T, typename A>
01795   inline  CollectiveIterator<T,A>& CollectiveIterator<T,A>::operator++()
01796   {
01797     typedef typename Map::iterator iterator;
01798     typedef typename Map::const_iterator const_iterator;
01799     const const_iterator end = map_.end();
01800     
01801     for(iterator iter = map_.begin(); iter != end;){
01802       // Step the iterator until we are >= index
01803       typename RemoteIndexList::const_iterator current = iter->second.first;
01804       typename RemoteIndexList::const_iterator rend = iter->second.second;
01805 
01806       // move all iterators pointing to the current global index to next value
01807       if(iter->second.first->localIndexPair().global()==index_ &&
01808          (noattribute || iter->second.first->localIndexPair().local().attribute() == attribute_))
01809         ++(iter->second.first);
01810       
01811       // erase from the map if there are no more entries.
01812       if(iter->second.first == iter->second.second)
01813         map_.erase(iter++);
01814       else{
01815         ++iter;
01816       }
01817     }
01818     return *this;
01819   }
01820 
01821   template<typename T, typename A>
01822   inline bool CollectiveIterator<T,A>::empty()
01823   {
01824     return map_.empty();
01825   }
01826     
01827   template<typename T, typename A>
01828   inline typename CollectiveIterator<T,A>::iterator
01829   CollectiveIterator<T,A>::begin()
01830   {
01831     if(noattribute)
01832       return iterator(map_.begin(), map_.end(), index_);
01833     else
01834       return iterator(map_.begin(), map_.end(), index_,
01835                       attribute_);
01836   }
01837    
01838   template<typename T, typename A>
01839   inline typename CollectiveIterator<T,A>::iterator
01840   CollectiveIterator<T,A>::end()
01841   {
01842     return iterator(map_.end(), map_.end(), index_);
01843   }
01844   
01845   template<typename TG, typename TA>
01846   inline std::ostream& operator<<(std::ostream& os, const RemoteIndex<TG,TA>& index)
01847   {
01848     os<<"[global="<<index.localIndexPair().global()<<", remote attribute="<<index.attribute()<<" local attribute="<<index.localIndexPair().local().attribute()<<"]";
01849     return os;
01850   }
01851   
01852   template<typename T, typename A>
01853   inline std::ostream& operator<<(std::ostream& os, const RemoteIndices<T,A>& indices)
01854   {
01855     int rank;
01856     MPI_Comm_rank(indices.comm_, &rank);
01857 
01858     typedef typename RemoteIndices<T,A>::RemoteIndexList RList;
01859     typedef typename std::map<int,std::pair<RList*,RList*> >::const_iterator const_iterator;
01860 
01861     const const_iterator rend = indices.remoteIndices_.end();
01862 
01863     for(const_iterator rindex = indices.remoteIndices_.begin(); rindex!=rend; ++rindex){
01864       os<<rank<<": Prozess "<<rindex->first<<":";
01865       
01866       if(!rindex->second.first->empty()){
01867         os<<" send:";
01868 
01869         const typename RList::const_iterator send= rindex->second.first->end();
01870       
01871         for(typename RList::const_iterator index = rindex->second.first->begin(); 
01872             index != send; ++index)
01873           os<<*index<<" ";
01874         os<<std::endl;
01875       }
01876       if(!rindex->second.second->empty()){
01877         os<<rank<<": Prozess "<<rindex->first<<": "<<"receive: ";
01878         
01879         const typename RList::const_iterator rend= rindex->second.second->end();
01880         
01881         for(typename RList::const_iterator index = rindex->second.second->begin(); 
01882             index != rend; ++index)
01883           os<<*index<<" ";
01884       }
01885       os<<std::endl<<std::flush;
01886     }
01887     return os;
01888   }
01890 }
01891 
01892 #endif
01893 #endif