dune-common  2.2.0
indicessyncer.hh
Go to the documentation of this file.
00001 // $Id$
00002 #ifndef DUNE_INDICESSYNCER_HH
00003 #define DUNE_INDICESSYNCER_HH
00004 
00005 #include"indexset.hh"
00006 #include"remoteindices.hh"
00007 #include<dune/common/stdstreams.hh>
00008 #include<dune/common/tuples.hh>
00009 #include<dune/common/sllist.hh>
00010 #include<cassert>
00011 #include<cmath>
00012 #include<limits>
00013 #include<algorithm>
00014 #include<functional>
00015 #include<map>
00016 
00017 #if HAVE_MPI
00018 namespace Dune
00019 {
00036   template<typename T>
00037   class IndicesSyncer
00038   {
00039   public:
00040 
00042     typedef T ParallelIndexSet;
00043 
00045     typedef typename ParallelIndexSet::IndexPair IndexPair;
00046     
00048     typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00049     
00051     typedef typename ParallelIndexSet::LocalIndex::Attribute Attribute;
00052 
00056     typedef Dune::RemoteIndices<ParallelIndexSet> RemoteIndices;
00057     
00067     IndicesSyncer(ParallelIndexSet& indexSet, 
00068                   RemoteIndices& remoteIndices);
00069     
00079     void sync();
00080 
00091     template<typename T1>
00092     void sync(T1& numberer);
00093     
00094   private:    
00095     
00097     ParallelIndexSet& indexSet_;
00098 
00100     RemoteIndices& remoteIndices_;
00101     
00103     char** sendBuffers_;
00104 
00106     char* receiveBuffer_;
00107     
00109     std::size_t* sendBufferSizes_;
00110     
00112     int receiveBufferSize_; // int because of MPI
00113     
00117     struct MessageInformation
00118     {
00119       MessageInformation()
00120         : publish(), pairs()
00121       {}
00123       int publish;
00128       int pairs;
00129     };
00130 
00134     class DefaultNumberer
00135     {
00136     public:
00142       std::size_t operator()(const GlobalIndex& global)
00143       {
00144         return std::numeric_limits<size_t>::max();
00145       }
00146     };
00147     
00149     MPI_Datatype datatype_;
00150 
00152     int rank_;
00153     
00158     typedef SLList<std::pair<GlobalIndex,Attribute>, typename RemoteIndices::Allocator> GlobalIndexList;
00159     
00161     typedef typename GlobalIndexList::ModifyIterator GlobalIndexModifier;
00162     
00166     typedef typename SLList<GlobalIndex, typename RemoteIndices::Allocator>::iterator 
00167     GlobalIndexIterator;
00168 
00170     typedef std::map<int, GlobalIndexList> GlobalIndicesMap;
00171     
00180     GlobalIndicesMap globalMap_;
00181 
00185     typedef SLList<bool, typename RemoteIndices::Allocator> BoolList;
00186 
00190     typedef typename BoolList::iterator BoolIterator;
00191 
00193     typedef typename BoolList::ModifyIterator BoolListModifier;
00194     
00196     typedef std::map<int,BoolList> BoolMap;
00197 
00202     BoolMap oldMap_;
00203     
00205     std::map<int,MessageInformation> infoSend_;
00206     
00208     typedef typename RemoteIndices::RemoteIndexList RemoteIndexList;
00209 
00211     typedef typename RemoteIndexList::ModifyIterator RemoteIndexModifier;
00212     
00214     typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
00215     
00217     typedef typename RemoteIndexList::iterator RemoteIndexIterator;
00218 
00220     typedef typename RemoteIndexList::const_iterator ConstRemoteIndexIterator;
00221 
00223     typedef Dune::tuple<RemoteIndexModifier,GlobalIndexModifier,BoolListModifier,
00224                   const ConstRemoteIndexIterator> IteratorTuple;
00225 
00233     class Iterators
00234     {
00235       friend class IndicesSyncer<T>;
00236     public:
00246       Iterators(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
00247                 BoolList& booleans);
00248       
00252       Iterators();
00253       
00257       Iterators& operator++();
00258       
00264       void insert(const RemoteIndex& index, 
00265                   const std::pair<GlobalIndex,Attribute>& global);
00266       
00271       RemoteIndex& remoteIndex() const;
00272       
00277       std::pair<GlobalIndex,Attribute>& globalIndexPair() const;
00278 
00279       Attribute& attribute() const;
00280       
00286       bool isOld() const;
00287       
00297       void reset(RemoteIndexList& remoteIndices, GlobalIndexList& globalIndices,
00298                  BoolList& booleans);
00299 
00305       bool isNotAtEnd() const;
00306       
00312       bool isAtEnd() const;
00313       
00314     private:
00324       IteratorTuple iterators_;
00325     };
00326     
00328     typedef std::map<int,Iterators> IteratorsMap;
00329     
00341     IteratorsMap iteratorsMap_;
00342         
00344     void calculateMessageSizes();
00345 
00353     void packAndSend(int destination, char* buffer, std::size_t bufferSize, MPI_Request& req);
00354     
00359     template<typename T1>
00360     void recvAndUnpack(T1& numberer);
00361 
00365     void registerMessageDatatype();
00366     
00370     void insertIntoRemoteIndexList(int process, 
00371                                    const std::pair<GlobalIndex,Attribute>& global, 
00372                                    char attribute);
00373 
00377     void resetIteratorsMap();
00378     
00383     bool checkReset();
00384     
00393     bool checkReset(const Iterators& iterators, RemoteIndexList& rlist, GlobalIndexList& gList,
00394                     BoolList& bList);
00395   };
00396 
00397   template<typename TG, typename TA>
00398   bool operator<(const IndexPair<TG,ParallelLocalIndex<TA> >& i1,
00399                  const std::pair<TG,TA>& i2)
00400   {
00401     return i1.global() < i2.first || 
00402       (i1.global() == i2.first && i1.local().attribute()<i2.second);
00403   }
00404 
00405   template<typename TG, typename TA>
00406   bool operator<(const std::pair<TG,TA>& i1,
00407                   const IndexPair<TG,ParallelLocalIndex<TA> >& i2)
00408   {
00409     return i1.first < i2.global() || 
00410       (i1.first == i2.global() && i1.second<i2.local().attribute());
00411   }
00412 
00413   template<typename TG, typename TA>
00414   bool operator==(const IndexPair<TG,ParallelLocalIndex<TA> >& i1, 
00415                   const std::pair<TG,TA>& i2)
00416   {
00417     return (i1.global() == i2.first && i1.local().attribute()==i2.second);
00418   }
00419   
00420   template<typename TG, typename TA>
00421   bool operator!=(const IndexPair<TG,ParallelLocalIndex<TA> >& i1, 
00422                   const std::pair<TG,TA>& i2)
00423   {
00424     return (i1.global() != i2.first || i1.local().attribute()!=i2.second);
00425   }
00426 
00427   template<typename TG, typename TA>
00428   bool operator==(const std::pair<TG,TA>& i2,
00429                   const IndexPair<TG,ParallelLocalIndex<TA> >& i1)
00430   {
00431     return (i1.global() == i2.first && i1.local().attribute()==i2.second);
00432   }
00433   
00434   template<typename TG, typename TA>
00435   bool operator!=(const std::pair<TG,TA>& i2,
00436                   const IndexPair<TG,ParallelLocalIndex<TA> >& i1)
00437   {
00438     return (i1.global() != i2.first || i1.local().attribute()!=i2.second);
00439   }
00440 
00457   template<typename T, typename A, typename A1>
00458   void storeGlobalIndicesOfRemoteIndices(std::map<int,SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> >& globalMap,
00459                                        const RemoteIndices<T,A1>& remoteIndices,
00460                                        const T& indexSet)
00461   {
00462     typedef typename RemoteIndices<T,A1>::const_iterator RemoteIterator;
00463     
00464     for(RemoteIterator remote = remoteIndices.begin(), end =remoteIndices.end(); remote != end; ++remote){
00465       typedef typename RemoteIndices<T,A1>::RemoteIndexList RemoteIndexList;
00466       typedef typename RemoteIndexList::const_iterator RemoteIndexIterator;
00467       typedef SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> GlobalIndexList;
00468       GlobalIndexList& global = globalMap[remote->first];
00469       RemoteIndexList& rList = *(remote->second.first);
00470       
00471       for(RemoteIndexIterator index = rList.begin(), riEnd = rList.end();
00472           index != riEnd; ++index){
00473         global.push_back(std::make_pair(index->localIndexPair().global(),
00474                                         index->localIndexPair().local().attribute()));
00475       }
00476     }
00477   }
00478   
00487   template<typename T, typename A, typename A1>
00488   inline void repairLocalIndexPointers(std::map<int,
00489                                        SLList<std::pair<typename T::GlobalIndex,
00490                                        typename T::LocalIndex::Attribute>,A> >& globalMap,
00491                                        RemoteIndices<T,A1>& remoteIndices,
00492                                        const T& indexSet)
00493   {
00494     typedef typename RemoteIndices<T,A1>::RemoteIndexMap::iterator RemoteIterator;
00495     typedef typename RemoteIndices<T,A1>::RemoteIndexList::iterator RemoteIndexIterator;
00496     typedef typename T::GlobalIndex GlobalIndex;
00497     typedef typename T::LocalIndex::Attribute Attribute;
00498     typedef std::pair<GlobalIndex,Attribute> GlobalIndexPair;
00499     typedef SLList<GlobalIndexPair,A> GlobalIndexPairList;
00500     typedef typename GlobalIndexPairList::iterator GlobalIndexIterator;
00501 
00502     assert(globalMap.size()==static_cast<std::size_t>(remoteIndices.neighbours()));
00503     // Repair pointers to index set in remote indices.
00504     typename std::map<int,GlobalIndexPairList>::iterator global = globalMap.begin();
00505     RemoteIterator end = remoteIndices.remoteIndices_.end();
00506     
00507     for(RemoteIterator remote = remoteIndices.remoteIndices_.begin(); remote != end; ++remote, ++global){
00508       typedef typename T::const_iterator IndexIterator;
00509 
00510       assert(remote->first==global->first);
00511       assert(remote->second.first->size() == global->second.size());
00512       
00513       RemoteIndexIterator riEnd  = remote->second.first->end();
00514       RemoteIndexIterator rIndex = remote->second.first->begin();
00515       GlobalIndexIterator gIndex = global->second.begin();
00516       IndexIterator       index  = indexSet.begin();
00517 
00518       assert(rIndex==riEnd || gIndex != global->second.end());
00519       while(rIndex != riEnd){
00520         // Search for the index in the set.
00521         assert(gIndex != global->second.end());
00522         
00523         while(!(index->global() == gIndex->first
00524                 && index->local().attribute() == gIndex->second)){
00525           ++index;
00526           // this is only needed for ALU, where there may exist 
00527           // more entries with the same global index in the remote index set
00528           // than in the index set
00529           if (index->global() > gIndex->first){
00530             index=indexSet.begin();
00531           }
00532         }
00533         
00534         assert(index != indexSet.end() && *index == *gIndex);
00535         
00536         rIndex->localIndex_ = &(*index);
00537         ++index;
00538         ++rIndex;
00539         ++gIndex;
00540       }
00541     }
00542     remoteIndices.sourceSeqNo_ = remoteIndices.source_->seqNo();
00543     remoteIndices.destSeqNo_ = remoteIndices.target_->seqNo();
00544   }
00545   
00546   template<typename T>
00547   IndicesSyncer<T>::IndicesSyncer(ParallelIndexSet& indexSet, 
00548                                       RemoteIndices& remoteIndices)
00549     : indexSet_(indexSet), remoteIndices_(remoteIndices)
00550   {
00551     // index sets must match.
00552     assert(remoteIndices.source_ == remoteIndices.target_);
00553     assert(remoteIndices.source_ == &indexSet);
00554     MPI_Comm_rank(remoteIndices_.communicator(), &rank_);
00555   }
00556     
00557   template<typename T>
00558   IndicesSyncer<T>::Iterators::Iterators(RemoteIndexList& remoteIndices, 
00559                                                GlobalIndexList& globalIndices,
00560                                                BoolList& booleans)
00561     : iterators_(remoteIndices.beginModify(), globalIndices.beginModify(),
00562                     booleans.beginModify(), remoteIndices.end())
00563   { }
00564 
00565   template<typename T>
00566   IndicesSyncer<T>::Iterators::Iterators()
00567     : iterators_()
00568   {}
00569   
00570   template<typename T>
00571   inline typename IndicesSyncer<T>::Iterators& IndicesSyncer<T>::Iterators::operator++()
00572   {
00573       ++(get<0>(iterators_));
00574       ++(get<1>(iterators_));
00575       ++(get<2>(iterators_));
00576       return *this;
00577   }
00578 
00579   template<typename T>
00580   inline void IndicesSyncer<T>::Iterators::insert(const RemoteIndex& index, 
00581                                                  const std::pair<GlobalIndex,Attribute>& global)
00582   {
00583     get<0>(iterators_).insert(index);
00584     get<1>(iterators_).insert(global);
00585     get<2>(iterators_).insert(false);
00586   }
00587   
00588   template<typename T>
00589   inline typename IndicesSyncer<T>::RemoteIndex& 
00590   IndicesSyncer<T>::Iterators::remoteIndex() const
00591   {
00592     return *(get<0>(iterators_));
00593   }
00594   
00595   template<typename T>
00596     inline std::pair<typename IndicesSyncer<T>::GlobalIndex,typename IndicesSyncer<T>::Attribute>&  
00597     IndicesSyncer<T>::Iterators::globalIndexPair() const
00598   {
00599     return *(get<1>(iterators_));
00600   }
00601   
00602   template<typename T>
00603   inline bool IndicesSyncer<T>::Iterators::isOld() const
00604   {
00605     return *(get<2>(iterators_));
00606   }
00607   
00608   template<typename T>
00609   inline void IndicesSyncer<T>::Iterators::reset(RemoteIndexList& remoteIndices, 
00610                                                 GlobalIndexList& globalIndices,
00611                                                 BoolList& booleans)
00612   {
00613     get<0>(iterators_) = remoteIndices.beginModify();
00614     get<1>(iterators_) = globalIndices.beginModify();
00615     get<2>(iterators_) = booleans.beginModify();
00616   }
00617   
00618   template<typename T>
00619   inline bool IndicesSyncer<T>::Iterators::isNotAtEnd() const
00620   {
00621     return get<0>(iterators_)!=get<3>(iterators_);
00622   }
00623   
00624   template<typename T>
00625   inline bool IndicesSyncer<T>::Iterators::isAtEnd() const
00626   {
00627     return get<0>(iterators_)==get<3>(iterators_);
00628   }
00629 
00630   template<typename T>
00631   void IndicesSyncer<T>::registerMessageDatatype()
00632   {
00633     MPI_Datatype type[2] = {MPI_INT, MPI_INT};
00634     int blocklength[2] = {1,1};
00635     MPI_Aint displacement[2];
00636     MPI_Aint base;
00637     
00638     // Compute displacement
00639     MessageInformation message;
00640     
00641     MPI_Address( &(message.publish), displacement);
00642     MPI_Address( &(message.pairs), displacement+1);
00643     
00644     // Make the displacement relative
00645     MPI_Address(&message, &base);
00646     displacement[0] -= base;
00647     displacement[1] -= base;
00648     
00649     MPI_Type_struct( 2, blocklength, displacement, type, &datatype_); 
00650     MPI_Type_commit(&datatype_);
00651   }
00652   
00653   template<typename T>
00654   void IndicesSyncer<T>::calculateMessageSizes()
00655   {    
00656     typedef typename ParallelIndexSet::const_iterator IndexIterator;
00657     typedef CollectiveIterator<T,typename RemoteIndices::Allocator> CollectiveIterator;
00658     
00659     IndexIterator iEnd = indexSet_.end();
00660     CollectiveIterator collIter = remoteIndices_.template iterator<true>();
00661     
00662     for(IndexIterator index = indexSet_.begin(); index != iEnd; ++index){
00663       collIter.advance(index->global(), index->local().attribute());
00664       if(collIter.empty())
00665         break;
00666       int knownRemote=0;
00667 
00668       typedef typename CollectiveIterator::iterator ValidIterator;
00669       ValidIterator end = collIter.end();
00670 
00671       // Count the remote indices we know.
00672       for(ValidIterator valid = collIter.begin(); valid != end; ++valid){
00673         ++knownRemote;
00674       }
00675       
00676       if(knownRemote>0){
00677         Dune::dverb<<rank_<<": publishing "<<knownRemote<<" for index "<<index->global()<< " for processes ";
00678         
00679         // Update MessageInformation
00680         for(ValidIterator valid = collIter.begin(); valid != end; ++valid){
00681           ++(infoSend_[valid.process()].publish);
00682           (infoSend_[valid.process()].pairs) += knownRemote;
00683           Dune::dverb<<valid.process()<<" ";
00684           Dune::dverb<<"(publish="<<infoSend_[valid.process()].publish<<", pairs="<<infoSend_[valid.process()].pairs
00685                <<") ";
00686         }
00687         Dune::dverb<<std::endl;
00688       }
00689     }
00690     
00691     typedef typename std::map<int,MessageInformation>::const_iterator 
00692       MessageIterator;
00693     
00694     const MessageIterator end = infoSend_.end();
00695 
00696     // registerMessageDatatype();
00697     
00698     // Now determine the buffersizes needed for each neighbour using MPI_Pack_size
00699     MessageInformation dummy;
00700     
00701     MessageIterator messageIter= infoSend_.begin();
00702 
00703     typedef typename RemoteIndices::RemoteIndexMap::const_iterator RemoteIterator;
00704     const RemoteIterator rend = remoteIndices_.end();
00705     int neighbour=0;
00706     
00707     for(RemoteIterator remote = remoteIndices_.begin(); remote != rend; ++remote, ++neighbour){
00708       MessageInformation* message;
00709       MessageInformation recv;
00710       
00711       if(messageIter != end && messageIter->first==remote->first){
00712         // We want to send message information to that process
00713         message = const_cast<MessageInformation*>(&(messageIter->second));
00714         ++messageIter;
00715       }else
00716         // We do not want to send information but the other process might.
00717         message = &dummy;
00718 
00719       sendBufferSizes_[neighbour]=0;
00720       int tsize;
00721       // The number of indices published
00722       MPI_Pack_size(1, MPI_INT,remoteIndices_.communicator(), &tsize);
00723       sendBufferSizes_[neighbour] += tsize;
00724 
00725       for(int i=0; i < message->publish; ++i){
00726         // The global index
00727         MPI_Pack_size(1, MPITraits<GlobalIndex>::getType(), remoteIndices_.communicator(), &tsize);
00728         sendBufferSizes_[neighbour] += tsize;
00729         // The attribute in the local index
00730         MPI_Pack_size(1, MPI_CHAR, remoteIndices_.communicator(), &tsize);
00731         sendBufferSizes_[neighbour] += tsize;
00732         // The number of corresponding remote indices
00733         MPI_Pack_size(1, MPI_INT, remoteIndices_.communicator(), &tsize);
00734         sendBufferSizes_[neighbour] += tsize;
00735       }
00736       for(int i=0; i < message->pairs; ++i){
00737         // The process of the remote index
00738         MPI_Pack_size(1, MPI_INT, remoteIndices_.communicator(), &tsize);
00739         sendBufferSizes_[neighbour] += tsize;
00740         // The attribute of the remote index
00741         MPI_Pack_size(1, MPI_CHAR, remoteIndices_.communicator(), &tsize);
00742         sendBufferSizes_[neighbour] += tsize;
00743       }
00744           
00745       Dune::dverb<<rank_<<": Buffer (neighbour="<<remote->first<<") size is "<< sendBufferSizes_[neighbour]<<" for publish="<<message->publish<<" pairs="<<message->pairs<<std::endl;
00746     }
00747 
00748   }
00749   
00750   template<typename T>
00751   inline void IndicesSyncer<T>::sync()
00752   {
00753     DefaultNumberer numberer;
00754     sync(numberer);
00755   }
00756   
00757   template<typename T>
00758   template<typename T1>
00759   void IndicesSyncer<T>::sync(T1& numberer)
00760   {
00761     
00762     // The pointers to the local indices in the remote indices
00763     // will become invalid due to the resorting of the index set.
00764     // Therefore store the corresponding global indices.
00765     // Mark all indices as not added
00766     
00767     typedef typename RemoteIndices::RemoteIndexMap::const_iterator 
00768       RemoteIterator;
00769 
00770     const RemoteIterator end = remoteIndices_.end();
00771     
00772     // Number of neighbours might change during the syncing.
00773     // save the old neighbours
00774     std::size_t noOldNeighbours = remoteIndices_.neighbours();
00775     int* oldNeighbours = new int[noOldNeighbours];
00776     sendBufferSizes_ = new std::size_t[noOldNeighbours];
00777     std::size_t i=0;
00778     
00779     for(RemoteIterator remote = remoteIndices_.begin(); remote != end; ++remote, ++i){
00780       typedef typename RemoteIndices::RemoteIndexList::const_iterator
00781         RemoteIndexIterator;
00782 
00783       oldNeighbours[i]=remote->first;
00784 
00785       // Make sure we only have one remote index list.
00786       assert(remote->second.first==remote->second.second);
00787 
00788       RemoteIndexList& rList = *(remote->second.first);
00789             
00790       // Store the corresponding global indices.
00791       GlobalIndexList& global = globalMap_[remote->first];
00792       BoolList& added = oldMap_[remote->first];
00793       RemoteIndexIterator riEnd = rList.end();
00794       
00795       for(RemoteIndexIterator index = rList.begin();
00796           index != riEnd; ++index){
00797         global.push_back(std::make_pair(index->localIndexPair().global(),
00798                                         index->localIndexPair().local().attribute()));
00799         added.push_back(true);
00800       }
00801 
00802       Iterators iterators(rList, global, added);
00803       iteratorsMap_.insert(std::make_pair(remote->first, iterators));
00804       assert(checkReset(iteratorsMap_[remote->first], rList,global,added));
00805     }
00806     
00807     // Exchange indices with each neighbour
00808     calculateMessageSizes();
00809     
00810     // Allocate the buffers
00811     receiveBufferSize_=1;
00812     sendBuffers_ = new char*[noOldNeighbours];
00813     
00814     for(std::size_t i=0; i<noOldNeighbours; ++i){
00815       sendBuffers_[i] = new char[sendBufferSizes_[i]];
00816       receiveBufferSize_ = std::max(receiveBufferSize_, static_cast<int>(sendBufferSizes_[i]));
00817     }
00818     
00819     receiveBuffer_=new char[receiveBufferSize_];
00820     
00821     indexSet_.beginResize();
00822     
00823     Dune::dverb<<rank_<<": Neighbours: ";
00824     
00825     for(i = 0; i<noOldNeighbours; ++i)
00826       Dune::dverb<<oldNeighbours[i]<<" ";
00827     
00828     Dune::dverb<<std::endl;
00829 
00830     MPI_Request* requests = new MPI_Request[noOldNeighbours];
00831     MPI_Status* statuses = new MPI_Status[noOldNeighbours];
00832 
00833     // Pack Message data and start the sends
00834     for(i = 0; i<noOldNeighbours; ++i)
00835         packAndSend(oldNeighbours[i], sendBuffers_[i], sendBufferSizes_[i], requests[i]);
00836 
00837     // Probe for incoming messages, receive and unpack them
00838     for(i = 0; i<noOldNeighbours; ++i)
00839       recvAndUnpack(numberer);
00840 //       }else{
00841 //      recvAndUnpack(oldNeighbours[i], numberer);
00842 //      packAndSend(oldNeighbours[i]);
00843 //       }
00844 //     }
00845 
00846     delete[] receiveBuffer_;
00847     
00848     // Wait for the completion of the sends
00849     // Wait for completion of sends
00850     if(MPI_SUCCESS!=MPI_Waitall(noOldNeighbours, requests, statuses)){
00851       std::cerr<<": MPI_Error occurred while sending message"<<std::endl;
00852       for(i=0;i< noOldNeighbours;i++)
00853         if(MPI_SUCCESS!=statuses[i].MPI_ERROR)
00854           std::cerr<<"Destination "<<statuses[i].MPI_SOURCE<<" error code: "<<statuses[i].MPI_ERROR<<std::endl;
00855     }
00856   
00857     delete[] statuses;
00858     delete[] requests;
00859     
00860     for(std::size_t i=0; i<noOldNeighbours; ++i)
00861       delete[] sendBuffers_[i];
00862 
00863     delete[] sendBuffers_;
00864     delete[] sendBufferSizes_;
00865     
00866     // No need for the iterator tuples any more
00867     iteratorsMap_.clear();
00868     
00869     indexSet_.endResize();
00870     
00871     delete[] oldNeighbours;
00872     
00873     repairLocalIndexPointers(globalMap_, remoteIndices_, indexSet_);
00874     
00875     oldMap_.clear();    
00876     globalMap_.clear();
00877         
00878     // update the sequence number
00879     remoteIndices_.sourceSeqNo_ = remoteIndices_.destSeqNo_ = indexSet_.seqNo();    
00880   }
00881   
00882   template<typename T>
00883   void IndicesSyncer<T>::packAndSend(int destination, char* buffer, std::size_t bufferSize, MPI_Request& request)
00884   {
00885     typedef typename ParallelIndexSet::const_iterator IndexIterator;
00886     typedef typename RemoteIndexList::const_iterator RemoteIndexIterator;
00887     typedef typename GlobalIndexList::const_iterator GlobalIterator;
00888     typedef typename BoolList::const_iterator BoolIterator;
00889         
00890     IndexIterator      iEnd       = indexSet_.end();
00891     int                bpos       = 0;
00892     int                published  = 0;
00893     int                pairs      = 0;
00894     
00895     assert(checkReset());
00896     
00897     // Pack the number of indices we publish
00898     MPI_Pack(&(infoSend_[destination].publish), 1, MPI_INT, buffer, bufferSize, &bpos, 
00899              remoteIndices_.communicator());
00900 
00901     for(IndexIterator index = indexSet_.begin(); index != iEnd; ++index){
00902       // Search for corresponding remote indices in all iterator tuples
00903       typedef typename IteratorsMap::iterator Iterator;
00904       Iterator iteratorsEnd = iteratorsMap_.end();
00905 
00906       // advance all iterators to a position with global index >= index->global()
00907       for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators){
00908         while(iterators->second.isNotAtEnd() && 
00909               iterators->second.globalIndexPair().first < index->global())
00910           ++(iterators->second);
00911         assert(!iterators->second.isNotAtEnd() || iterators->second.globalIndexPair().first >= index->global());
00912       }
00913       
00914       // Add all remote indices positioned at global which were already present before calling sync
00915       // to the message.
00916       // Count how many remote indices we will send
00917       int indices = 0;
00918       bool knownRemote = false; // Is the remote process supposed to know this index?
00919             
00920       for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators)
00921         {
00922           std::pair<GlobalIndex,Attribute> p;
00923           if (iterators->second.isNotAtEnd())
00924             {
00925               p = iterators->second.globalIndexPair();
00926             }
00927           
00928         if(iterators->second.isNotAtEnd() && iterators->second.isOld() 
00929            && iterators->second.globalIndexPair().first == index->global()){
00930           indices++;
00931           if(destination == iterators->first)
00932             knownRemote = true;
00933         }
00934         }
00935       
00936       if(!knownRemote)
00937         // We do not need to send any indices
00938         continue;
00939       
00940       Dune::dverb<<rank_<<": sending "<<indices<<" for index "<<index->global()<<" to "<<destination<<std::endl;
00941 
00942 
00943       // Pack the global index, the attribute and the number
00944       MPI_Pack(const_cast<GlobalIndex*>(&(index->global())), 1, MPITraits<GlobalIndex>::getType(), buffer, bufferSize, &bpos, 
00945                remoteIndices_.communicator());
00946       
00947       char attr = index->local().attribute();
00948       MPI_Pack(&attr, 1, MPI_CHAR, buffer, bufferSize, &bpos, 
00949                remoteIndices_.communicator());
00950 
00951       // Pack the number of remote indices we send.
00952       MPI_Pack(&indices, 1, MPI_INT, buffer, bufferSize, &bpos, 
00953                remoteIndices_.communicator());
00954         
00955       // Pack the information about the remote indices
00956       for(Iterator iterators = iteratorsMap_.begin(); iteratorsEnd != iterators; ++iterators)
00957         if(iterators->second.isNotAtEnd() && iterators->second.isOld() 
00958            && iterators->second.globalIndexPair().first == index->global()){
00959           int process = iterators->first;
00960 
00961           ++pairs;
00962           assert(pairs <= infoSend_[destination].pairs);
00963           MPI_Pack(&process, 1, MPI_INT, buffer, bufferSize, &bpos, 
00964                    remoteIndices_.communicator());
00965           char attr = iterators->second.remoteIndex().attribute();
00966           
00967           MPI_Pack(&attr, 1, MPI_CHAR, buffer, bufferSize, &bpos, 
00968                    remoteIndices_.communicator());
00969           --indices;
00970         }
00971       assert(indices==0);
00972       ++published;
00973       Dune::dvverb<<" (publish="<<published<<", pairs="<<pairs<<")"<<std::endl;
00974       assert(published <= infoSend_[destination].publish);
00975     }
00976 
00977     // Make sure we send all expected entries
00978     assert(published == infoSend_[destination].publish);
00979     assert(pairs == infoSend_[destination].pairs);
00980     resetIteratorsMap();
00981 
00982     Dune::dverb << rank_<<": Sending message of "<<bpos<<" bytes to "<<destination<<std::endl;
00983     
00984     MPI_Issend(buffer, bpos, MPI_PACKED, destination, 345, remoteIndices_.communicator(),&request);
00985   }
00986 
00987   template<typename T>
00988     inline void IndicesSyncer<T>::insertIntoRemoteIndexList(int process, 
00989                                                             const std::pair<GlobalIndex,Attribute>& globalPair,
00990                                                             char attribute)
00991   {
00992     Dune::dverb<<"Inserting from "<<process<<" "<<globalPair.first<<", "<<
00993                 globalPair.second<<" "<<attribute<<std::endl;
00994 
00995     resetIteratorsMap();
00996 
00997     // There might be cases where there no remote indices for that process yet
00998     typename IteratorsMap::iterator found = iteratorsMap_.find(process);
00999     
01000     if( found == iteratorsMap_.end() ){
01001       Dune::dverb<<"Discovered new neighbour "<<process<<std::endl;
01002       RemoteIndexList* rlist = new RemoteIndexList();
01003       remoteIndices_.remoteIndices_.insert(std::make_pair(process,std::make_pair(rlist,rlist)));
01004       Iterators iterators = Iterators(*rlist, globalMap_[process], oldMap_[process]);
01005       found = iteratorsMap_.insert(std::make_pair(process, iterators)).first;
01006     }
01007     
01008     Iterators& iterators = found->second;
01009     
01010     // Search for the remote index
01011     while(iterators.isNotAtEnd() && iterators.globalIndexPair() < globalPair){
01012       // Increment all iterators
01013       ++iterators;
01014       
01015     }
01016     
01017     if(iterators.isAtEnd() || iterators.globalIndexPair() != globalPair){
01018       // The entry is not yet known
01019       // Insert in the the list and do not change the first iterator.
01020       iterators.insert(RemoteIndex(Attribute(attribute)),globalPair);
01021       return;
01022     }
01023     
01024     // Global indices match
01025     bool indexIsThere=false;
01026     for(Iterators tmpIterators = iterators; 
01027         !tmpIterators.isAtEnd() && tmpIterators.globalIndexPair() == globalPair;
01028         ++tmpIterators)
01029       //entry already exists with the same attribute
01030       if(tmpIterators.globalIndexPair().second == attribute){
01031         indexIsThere=true;
01032         break;
01033       }
01034     
01035     if(!indexIsThere)
01036       // The entry is not yet known
01037       // Insert in the the list and do not change the first iterator.
01038       iterators.insert(RemoteIndex(Attribute(attribute)),globalPair);
01039   }
01040   
01041   template<typename T>
01042   template<typename T1>
01043   void IndicesSyncer<T>::recvAndUnpack(T1& numberer)
01044   {
01045     typedef typename ParallelIndexSet::const_iterator IndexIterator;
01046     typedef typename RemoteIndexList::iterator RemoteIndexIterator;
01047     typedef typename GlobalIndexList::iterator GlobalIndexIterator;
01048 
01049     IndexIterator    iEnd   = indexSet_.end();
01050     IndexIterator    index  = indexSet_.begin();
01051     int              bpos   = 0;
01052     int              publish;
01053     
01054     assert(checkReset());
01055 
01056     MPI_Status status;
01057 
01058     // We have to determine the message size and source before the receive
01059     MPI_Probe(MPI_ANY_SOURCE, 345, remoteIndices_.communicator(), &status);
01060     
01061     int source=status.MPI_SOURCE;
01062     int count;
01063     MPI_Get_count(&status, MPI_PACKED, &count);
01064 
01065     Dune::dvverb<<rank_<<": Receiving message from "<< source<<" with "<<count<<" bytes"<<std::endl;
01066 
01067     if(count>receiveBufferSize_){
01068       receiveBufferSize_=count;
01069       delete[] receiveBuffer_;
01070       receiveBuffer_ = new char[receiveBufferSize_];
01071     }
01072 
01073     MPI_Recv(receiveBuffer_, count, MPI_PACKED, source, 345, remoteIndices_.communicator(), &status);
01074     
01075     // How many global entries were published?
01076     MPI_Unpack(receiveBuffer_,  count, &bpos, &publish, 1, MPI_INT, remoteIndices_.communicator());
01077     
01078     // Now unpack the remote indices and add them.
01079     while(publish>0){
01080       
01081       // Unpack information about the local index on the source process
01082       GlobalIndex  global;          // global index of the current entry
01083       char sourceAttribute; // Attribute on the source process
01084       int pairs;
01085       
01086       MPI_Unpack(receiveBuffer_,  count, &bpos, &global, 1, MPITraits<GlobalIndex>::getType(), 
01087                  remoteIndices_.communicator());
01088       MPI_Unpack(receiveBuffer_,  count, &bpos, &sourceAttribute, 1, MPI_CHAR, 
01089                  remoteIndices_.communicator());
01090       MPI_Unpack(receiveBuffer_,  count, &bpos, &pairs, 1, MPI_INT, 
01091                  remoteIndices_.communicator());
01092     
01093       // Insert the entry on the remote process to our
01094       // remote index list
01095       SLList<std::pair<int,Attribute> > sourceAttributeList;
01096       sourceAttributeList.push_back(std::make_pair(source,Attribute(sourceAttribute)));
01097 #ifndef NDEBUG
01098       bool foundSelf = false;
01099 #endif
01100       Attribute myAttribute=Attribute();
01101       
01102       // Unpack the remote indices
01103       for(;pairs>0;--pairs){
01104         // Unpack the process id that knows the index
01105         int process;
01106         char attribute;
01107         MPI_Unpack(receiveBuffer_,  count, &bpos, &process, 1, MPI_INT, 
01108                  remoteIndices_.communicator());
01109         // Unpack the attribute
01110         MPI_Unpack(receiveBuffer_,  count, &bpos, &attribute, 1, MPI_CHAR, 
01111                    remoteIndices_.communicator());
01112 
01113         if(process==rank_){
01114 #ifndef NDEBUG
01115           foundSelf=true;
01116 #endif
01117           myAttribute=Attribute(attribute);
01118           // Now we know the local attribute of the global index
01119           //Only add the index if it is unknown.
01120           // Do we know that global index already?
01121           IndexIterator pos = std::lower_bound(index, iEnd, IndexPair(global));
01122 
01123           if(pos == iEnd || pos->global() != global){
01124             // no entry with this global index
01125             indexSet_.add(global,
01126                           ParallelLocalIndex<Attribute>(numberer(global),
01127                                                         myAttribute, true));
01128             Dune::dvverb << "Adding "<<global<<" "<<myAttribute<<std::endl;
01129             continue;
01130           }
01131           
01132           // because of above the global indices match. Add only if the attribute is different
01133           bool indexIsThere = false;
01134           index=pos;
01135 
01136           for(;pos->global()==global;++pos)
01137             if(pos->local().attribute() == myAttribute){
01138               Dune::dvverb<<"found "<<global<<" "<<myAttribute<<std::endl;
01139               indexIsThere = true;
01140               break;
01141             }
01142           
01143           if(!indexIsThere){
01144             indexSet_.add(global,
01145                           ParallelLocalIndex<Attribute>(numberer(global),
01146                                                         myAttribute, true));  
01147             Dune::dvverb << "Adding "<<global<<" "<<myAttribute<<std::endl;
01148           }
01149           
01150         }else{
01151           sourceAttributeList.push_back(std::make_pair(process,Attribute(attribute)));
01152         }
01153       }
01154       assert(foundSelf);
01155       // Insert remote indices
01156       typedef typename SLList<std::pair<int,Attribute> >::const_iterator Iter;
01157       for(Iter i=sourceAttributeList.begin(), end=sourceAttributeList.end();
01158           i!=end; ++i)
01159         insertIntoRemoteIndexList(i->first, std::make_pair(global, myAttribute),
01160                                   i->second);
01161       --publish;
01162     }
01163     
01164     resetIteratorsMap();
01165   }
01166   
01167   template<typename T>
01168   void IndicesSyncer<T>::resetIteratorsMap(){
01169     
01170     // Reset iterators in all tuples.
01171     typedef typename IteratorsMap::iterator Iterator;
01172     typedef typename RemoteIndices::RemoteIndexMap::iterator 
01173       RemoteIterator;
01174     typedef typename GlobalIndicesMap::iterator GlobalIterator;
01175     typedef typename BoolMap::iterator BoolIterator;
01176     
01177     const RemoteIterator remoteEnd = remoteIndices_.remoteIndices_.end();
01178     Iterator iterators = iteratorsMap_.begin();
01179     GlobalIterator global = globalMap_.begin();
01180     BoolIterator added = oldMap_.begin();
01181     
01182     for(RemoteIterator remote = remoteIndices_.remoteIndices_.begin();
01183         remote != remoteEnd; ++remote, ++global, ++added, ++iterators){
01184       iterators->second.reset(*(remote->second.first), global->second, added->second);
01185     }
01186   }
01187  
01188   template<typename T>
01189   bool IndicesSyncer<T>::checkReset(const Iterators& iterators, RemoteIndexList& rList, GlobalIndexList& gList,
01190                                           BoolList& bList){
01191     
01192     if(get<0>(iterators.iterators_)!=rList.begin())
01193       return false;
01194     if(get<1>(iterators.iterators_)!=gList.begin())
01195       return false;
01196     if(get<2>(iterators.iterators_)!=bList.begin())
01197       return false;
01198     return true;
01199   }
01200   
01201 
01202   template<typename T>
01203   bool IndicesSyncer<T>::checkReset(){
01204     
01205     // Reset iterators in all tuples.
01206     typedef typename IteratorsMap::iterator Iterator;
01207     typedef typename RemoteIndices::RemoteIndexMap::iterator 
01208       RemoteIterator;
01209     typedef typename GlobalIndicesMap::iterator GlobalIterator;
01210     typedef typename BoolMap::iterator BoolIterator;
01211     
01212     const RemoteIterator remoteEnd = remoteIndices_.remoteIndices_.end();
01213     Iterator iterators = iteratorsMap_.begin();
01214     GlobalIterator global = globalMap_.begin();
01215     BoolIterator added = oldMap_.begin();
01216     bool ret = true;
01217     
01218     for(RemoteIterator remote = remoteIndices_.remoteIndices_.begin();
01219         remote != remoteEnd; ++remote, ++global, ++added, ++iterators){
01220       if(!checkReset(iterators->second, *(remote->second.first), global->second,
01221                      added->second))
01222         ret=false;
01223     }
01224     return ret;
01225   } 
01226 }
01227 
01228 #endif
01229 #endif