dune-common  2.2.1
remoteindices.hh
Go to the documentation of this file.
1 // $Id$
2 #ifndef DUNE_REMOTEINDICES_HH
3 #define DUNE_REMOTEINDICES_HH
4 
5 #include"indexset.hh"
7 #include"plocalindex.hh"
9 #include<dune/common/sllist.hh>
12 #include<map>
13 #include<set>
14 #include<utility>
15 #include<iostream>
16 #include<algorithm>
17 #include<iterator>
18 #if HAVE_MPI
19 #include <dune/common/mpitraits.hh>
20 #include"mpi.h"
21 
22 namespace Dune{
33 
34  template<typename TG, typename TA>
36  {
37  public:
38  inline static MPI_Datatype getType();
39  private:
40  static MPI_Datatype type;
41  };
42 
43 
44  template<typename T, typename A>
46 
47  template<typename T1, typename T2>
48  class RemoteIndex;
49 
50  template<typename T>
51  class IndicesSyncer;
52 
53  template<typename T1, typename T2>
54  std::ostream& operator<<(std::ostream& os, const RemoteIndex<T1,T2>& index);
55 
56 
57  template<typename T, typename A, bool mode>
59 
60 
64  template<typename T1, typename T2>
65  class RemoteIndex
66  {
67  template<typename T>
68  friend class IndicesSyncer;
69 
70  template<typename T, typename A, typename A1>
71  friend void repairLocalIndexPointers(std::map<int,SLList<std::pair<typename T::GlobalIndex, typename T::LocalIndex::Attribute>,A> >&,
73  const T&);
74 
75  template<typename T, typename A, bool mode>
77 
78  public:
83  typedef T1 GlobalIndex;
92  typedef T2 Attribute;
93 
99 
104  const Attribute attribute() const;
105 
111  const PairType& localIndexPair() const;
112 
116  RemoteIndex();
117 
118 
124  RemoteIndex(const T2& attribute,
125  const PairType* local);
126 
127 
133  RemoteIndex(const T2& attribute);
134 
135  bool operator==(const RemoteIndex& ri) const;
136 
137  bool operator!=(const RemoteIndex& ri) const;
138  private:
140  const PairType* localIndex_;
141 
143  char attribute_;
144  };
145 
146  template<class T, class A>
147  std::ostream& operator<<(std::ostream& os, const RemoteIndices<T,A>& indices);
148 
149  class InterfaceBuilder;
150 
151  template<class T, class A>
153 
154  template<class T>
155  class IndicesSyncer;
156 
157  // forward declaration needed for friend declaration.
158  template<typename T1, typename T2>
160 
161 
178  template<class T, class A=std::allocator<RemoteIndex<typename T::GlobalIndex,
179  typename T::LocalIndex::Attribute> > >
180  class RemoteIndices
181  {
182  friend class InterfaceBuilder;
183  friend class IndicesSyncer<T>;
184  template<typename T1, typename A2, typename A1>
185  friend void repairLocalIndexPointers(std::map<int,SLList<std::pair<typename T1::GlobalIndex, typename T1::LocalIndex::Attribute>,A2> >&,
187  const T1&);
188 
189  template<class G, class T1, class T2>
190  friend void fillIndexSetHoles(const G& graph, Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm);
191  friend std::ostream& operator<<<>(std::ostream&, const RemoteIndices<T>&);
192 
193  public:
194 
198  typedef T ParallelIndexSet;
199 
203 
208 
209 
214 
218  typedef typename LocalIndex::Attribute Attribute;
219 
224 
225 
229  typedef typename A::template rebind<RemoteIndex>::other Allocator;
230 
234 
236  typedef std::map<int, std::pair<RemoteIndexList*,RemoteIndexList*> >
238 
239  typedef typename RemoteIndexMap::const_iterator const_iterator;
240 
258  inline RemoteIndices(const ParallelIndexSet& source, const ParallelIndexSet& destination,
259  const MPI_Comm& comm, const std::vector<int>& neighbours=std::vector<int>(), bool includeSelf=false);
260 
261  RemoteIndices();
262 
270  void setIncludeSelf(bool includeSelf);
271 
288  void setIndexSets(const ParallelIndexSet& source, const ParallelIndexSet& destination,
289  const MPI_Comm& comm, const std::vector<int>& neighbours=std::vector<int>());
290 
291  template<typename C>
292  void setNeighbours(const C& neighbours)
293  {
294  neighbourIds.clear();
295  neighbourIds.insert(neighbours.begin(), neighbours.end());
296 
297  }
298 
299  const std::set<int>& getNeighbours() const
300  {
301  return neighbourIds;
302  }
303 
307  ~RemoteIndices();
308 
318  template<bool ignorePublic>
319  void rebuild();
320 
321  bool operator==(const RemoteIndices& ri);
322 
330  inline bool isSynced() const;
331 
335  inline MPI_Comm communicator() const;
336 
351  template<bool mode, bool send>
353 
360  inline const_iterator find(int proc) const;
361 
366  inline const_iterator begin() const;
367 
372  inline const_iterator end() const;
373 
377  template<bool send>
378  inline CollectiveIteratorT iterator() const;
379 
383  inline void free();
384 
389  inline int neighbours() const;
390 
392  inline const ParallelIndexSet& sourceIndexSet() const;
393 
395  inline const ParallelIndexSet& destinationIndexSet() const;
396 
397  private:
400  {}
401 
403  const ParallelIndexSet* source_;
404 
406  const ParallelIndexSet* target_;
407 
409  MPI_Comm comm_;
410 
413  std::set<int> neighbourIds;
414 
416  const static int commTag_=333;
417 
422  int sourceSeqNo_;
423 
428  int destSeqNo_;
429 
433  bool publicIgnored;
434 
438  bool firstBuild;
439 
440  /*
441  * @brief If true, sending from indices of the processor to other
442  * indices on the same processor is enabled even if the same indexset is used
443  * on both the
444  * sending and receiving side.
445  */
446  bool includeSelf;
447 
449  typedef IndexPair<GlobalIndex, LocalIndex>
450  PairType;
451 
458  RemoteIndexMap remoteIndices_;
459 
470  template<bool ignorePublic>
471  inline void buildRemote(bool includeSelf);
472 
478  inline int noPublic(const ParallelIndexSet& indexSet);
479 
491  template<bool ignorePublic>
492  inline void packEntries(PairType** myPairs, const ParallelIndexSet& indexSet,
493  char* p_out, MPI_Datatype type, int bufferSize,
494  int* position, int n);
495 
509  inline void unpackIndices(RemoteIndexList& remote, int remoteEntries,
510  PairType** local, int localEntries, char* p_in,
511  MPI_Datatype type, int* positon, int bufferSize,
512  bool fromOurself);
513 
514  inline void unpackIndices(RemoteIndexList& send, RemoteIndexList& receive,
515  int remoteEntries, PairType** localSource,
516  int localSourceEntries, PairType** localDest,
517  int localDestEntries, char* p_in,
518  MPI_Datatype type, int* position, int bufferSize);
519 
520  void unpackCreateRemote(char* p_in, PairType** sourcePairs, PairType** DestPairs,
521  int remoteProc, int sourcePublish, int destPublish,
522  int bufferSize, bool sendTwo, bool fromOurSelf=false);
523  };
524 
542  template<class T, class A, bool mode>
543  class RemoteIndexListModifier
544  {
545 
546  template<typename T1, typename A1>
547  friend class RemoteIndices;
548 
549  public:
551  {};
552 
553  enum {
563  };
564 
568  typedef T ParallelIndexSet;
569 
574 
579 
583  typedef typename LocalIndex::Attribute Attribute;
584 
589 
593  typedef A Allocator;
594 
598 
603 
608 
622  void insert(const RemoteIndex& index) throw(InvalidPosition);
623 
624 
639  void insert(const RemoteIndex& index, const GlobalIndex& global) throw(InvalidPosition);
640 
648  bool remove(const GlobalIndex& global) throw(InvalidPosition);
649 
663 
664 
666 
671  RemoteIndexListModifier()
672  : glist_()
673  {}
674 
677  {
678  if(glist_)
679  delete glist_;
680  }
681 
682  private:
683 
690  RemoteIndexList& rList);
691 
692  typedef SLList<GlobalIndex,Allocator> GlobalList;
693  typedef typename GlobalList::ModifyIterator GlobalModifyIterator;
694  RemoteIndices<ParallelIndexSet>* remoteIndices_;
695  RemoteIndexList* rList_;
696  const ParallelIndexSet* indexSet_;
697  GlobalList* glist_;
698  ModifyIterator iter_;
699  GlobalModifyIterator giter_;
700  ConstIterator end_;
701  bool first_;
702  GlobalIndex last_;
703  };
704 
709  template<class T, class A>
710  class CollectiveIterator
711  {
712 
716  typedef T ParallelIndexSet;
717 
721  typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
722 
726  typedef typename ParallelIndexSet::LocalIndex LocalIndex;
727 
731  typedef typename LocalIndex::Attribute Attribute;
732 
734  typedef Dune::RemoteIndex<GlobalIndex,Attribute> RemoteIndex;
735 
737  typedef typename A::template rebind<RemoteIndex>::other Allocator;
738 
740  typedef Dune::SLList<RemoteIndex,Allocator> RemoteIndexList;
741 
743  typedef std::map<int,std::pair<typename RemoteIndexList::const_iterator,
744  const typename RemoteIndexList::const_iterator> >
745  Map;
746 
747  public:
748 
750  typedef std::map<int, std::pair<RemoteIndexList*,RemoteIndexList*> >
752 
758  inline CollectiveIterator(const RemoteIndexMap& map_, bool send);
759 
768  inline void advance(const GlobalIndex& global);
769 
779  inline void advance(const GlobalIndex& global, const Attribute& attribute);
780 
782 
786  inline bool empty();
787 
794  class iterator
795  {
796  public:
797  typedef typename Map::iterator RealIterator;
798  typedef typename Map::iterator ConstRealIterator;
799 
800 
802  iterator(const RealIterator& iter, const ConstRealIterator& end, GlobalIndex& index)
803  : iter_(iter), end_(end), index_(index), hasAttribute(false)
804  {
805  // Move to the first valid entry
806  while(iter_!=end_ && iter_->second.first->localIndexPair().global()!=index_)
807  ++iter_;
808  }
809 
810  iterator(const RealIterator& iter, const ConstRealIterator& end, GlobalIndex index,
811  Attribute attribute)
812  : iter_(iter), end_(end), index_(index), attribute_(attribute), hasAttribute(true)
813  {
814  // Move to the first valid entry or the end
815  while(iter_!=end_ && (iter_->second.first->localIndexPair().global()!=index_
816  || iter_->second.first->localIndexPair().local().attribute()!=attribute))
817  ++iter_;
818  }
820  iterator(const iterator& other)
821  : iter_(other.iter_), end_(other.end_), index_(other.index_)
822  { }
823 
826  {
827  ++iter_;
828  // If entry is not valid move on
829  while(iter_!=end_ && (iter_->second.first->localIndexPair().global()!=index_ ||
830  (hasAttribute &&
831  iter_->second.first->localIndexPair().local().attribute()!=attribute_)))
832  ++iter_;
833  assert(iter_==end_ ||
834  (iter_->second.first->localIndexPair().global()==index_));
835  assert(iter_==end_ || !hasAttribute ||
836  (iter_->second.first->localIndexPair().local().attribute()==attribute_));
837  return *this;
838  }
839 
841  const RemoteIndex& operator*()const
842  {
843  return *(iter_->second.first);
844  }
845 
847  int process() const
848  {
849  return iter_->first;
850  }
851 
853  const RemoteIndex* operator->()const
854  {
855  return iter_->second.first.operator->();
856  }
857 
859  bool operator==(const iterator& other)
860  {
861  return other.iter_==iter_;
862  }
863 
865  bool operator!=(const iterator& other)
866  {
867  return other.iter_!=iter_;
868  }
869 
870  private:
871  iterator();
872 
873  RealIterator iter_;
874  RealIterator end_;
875  GlobalIndex index_;
876  Attribute attribute_;
877  bool hasAttribute;
878  };
879 
880  iterator begin();
881 
882  iterator end();
883 
884  private:
885 
886  Map map_;
887  GlobalIndex index_;
888  Attribute attribute_;
889  bool noattribute;
890  };
891 
892  template<typename TG, typename TA>
894  {
895  if(type==MPI_DATATYPE_NULL){
896  int length[4];
897  MPI_Aint disp[4];
898  MPI_Datatype types[4] = {MPI_LB, MPITraits<TG>::getType(),
899  MPITraits<ParallelLocalIndex<TA> >::getType(), MPI_UB};
901  length[0]=length[1]=length[2]=length[3]=1;
902  MPI_Address(rep, disp); // lower bound of the datatype
903  MPI_Address(&(rep[0].global_), disp+1);
904  MPI_Address(&(rep[0].local_), disp+2);
905  MPI_Address(rep+1, disp+3); // upper bound of the datatype
906  for(int i=3; i >= 0; --i)
907  disp[i] -= disp[0];
908  MPI_Type_struct(4, length, disp, types, &type);
909  MPI_Type_commit(&type);
910  }
911  return type;
912  }
913 
914  template<typename TG, typename TA>
915  MPI_Datatype MPITraits<IndexPair<TG,ParallelLocalIndex<TA> > >::type=MPI_DATATYPE_NULL;
916 
917  template<typename T1, typename T2>
918  RemoteIndex<T1,T2>::RemoteIndex(const T2& attribute, const PairType* local)
919  : localIndex_(local), attribute_(attribute)
920  {}
921 
922  template<typename T1, typename T2>
923  RemoteIndex<T1,T2>::RemoteIndex(const T2& attribute)
924  : localIndex_(0), attribute_(attribute)
925  {}
926 
927  template<typename T1, typename T2>
929  : localIndex_(0), attribute_()
930  {}
931  template<typename T1, typename T2>
932  inline bool RemoteIndex<T1,T2>::operator==(const RemoteIndex& ri) const
933  {
934  return localIndex_==ri.localIndex_ && attribute_==ri.attribute;
935  }
936 
937  template<typename T1, typename T2>
938  inline bool RemoteIndex<T1,T2>::operator!=(const RemoteIndex& ri) const
939  {
940  return localIndex_!=ri.localIndex_ || attribute_!=ri.attribute_;
941  }
942 
943  template<typename T1, typename T2>
944  inline const T2 RemoteIndex<T1,T2>::attribute() const
945  {
946  return T2(attribute_);
947  }
948 
949  template<typename T1, typename T2>
951  {
952  return *localIndex_;
953  }
954 
955  template<typename T, typename A>
957  const ParallelIndexSet& destination,
958  const MPI_Comm& comm,
959  const std::vector<int>& neighbours,
960  bool includeSelf_)
961  : source_(&source), target_(&destination), comm_(comm),
962  sourceSeqNo_(-1), destSeqNo_(-1), publicIgnored(false), firstBuild(true),
963  includeSelf(includeSelf_)
964  {
965  setNeighbours(neighbours);
966  }
967 
968  template<typename T, typename A>
970  {
971  includeSelf=b;
972  }
973 
974  template<typename T, typename A>
976  :source_(0), target_(0), sourceSeqNo_(-1),
977  destSeqNo_(-1), publicIgnored(false), firstBuild(true)
978  {}
979 
980  template<class T, typename A>
982  const ParallelIndexSet& destination,
983  const MPI_Comm& comm,
984  const std::vector<int>& neighbours)
985  {
986  free();
987  source_ = &source;
988  target_ = &destination;
989  comm_ = comm;
990  firstBuild = true;
991  setNeighbours(neighbours);
992  }
993 
994  template<typename T, typename A>
995  const typename RemoteIndices<T,A>::ParallelIndexSet&
997  {
998  return *source_;
999  }
1000 
1001 
1002  template<typename T, typename A>
1003  const typename RemoteIndices<T,A>::ParallelIndexSet&
1005  {
1006  return *target_;
1007  }
1008 
1009 
1010  template<typename T, typename A>
1012  {
1013  free();
1014  }
1015 
1016  template<typename T, typename A>
1017  template<bool ignorePublic>
1019  const ParallelIndexSet& indexSet,
1020  char* p_out, MPI_Datatype type,
1021  int bufferSize,
1022  int *position, int n)
1023  {
1024  // fill with own indices
1027  const const_iterator end = indexSet.end();
1028 
1029  //Now pack the source indices
1030  int i=0;
1031  for(const_iterator index = indexSet.begin(); index != end; ++index)
1032  if(ignorePublic || index->local().isPublic()){
1033 
1034  MPI_Pack(const_cast<PairType*>(&(*index)), 1,
1035  type,
1036  p_out, bufferSize, position, comm_);
1037  pairs[i++] = const_cast<PairType*>(&(*index));
1038 
1039  }
1040  assert(i==n);
1041  }
1042 
1043  template<typename T, typename A>
1044  inline int RemoteIndices<T,A>::noPublic(const ParallelIndexSet& indexSet)
1045  {
1046  typedef typename ParallelIndexSet::const_iterator const_iterator;
1047 
1048  int noPublic=0;
1049 
1050  const const_iterator end=indexSet.end();
1051  for(const_iterator index=indexSet.begin(); index!=end; ++index)
1052  if(index->local().isPublic())
1053  noPublic++;
1054 
1055  return noPublic;
1056 
1057  }
1058 
1059 
1060  template<typename T, typename A>
1061  inline void RemoteIndices<T,A>::unpackCreateRemote(char* p_in, PairType** sourcePairs,
1062  PairType** destPairs, int remoteProc,
1063  int sourcePublish, int destPublish,
1064  int bufferSize, bool sendTwo,
1065  bool fromOurSelf)
1066  {
1067 
1068  // unpack the number of indices we received
1069  int noRemoteSource=-1, noRemoteDest=-1;
1070  char twoIndexSets=0;
1071  int position=0;
1072  // Did we receive two index sets?
1073  MPI_Unpack(p_in, bufferSize, &position, &twoIndexSets, 1, MPI_CHAR, comm_);
1074  // The number of source indices received
1075  MPI_Unpack(p_in, bufferSize, &position, &noRemoteSource, 1, MPI_INT, comm_);
1076  // The number of destination indices received
1077  MPI_Unpack(p_in, bufferSize, &position, &noRemoteDest, 1, MPI_INT, comm_);
1078 
1079 
1080  // Indices for which we receive
1081  RemoteIndexList* receive= new RemoteIndexList();
1082  // Indices for which we send
1083  RemoteIndexList* send=0;
1084 
1085  MPI_Datatype type= MPITraits<PairType>::getType();
1086 
1087  if(!twoIndexSets){
1088  if(sendTwo){
1089  send = new RemoteIndexList();
1090  // Create both remote index sets simultaneously
1091  unpackIndices(*send, *receive, noRemoteSource, sourcePairs, sourcePublish,
1092  destPairs, destPublish, p_in, type, &position, bufferSize);
1093  }else{
1094  // we only need one list
1095  unpackIndices(*receive, noRemoteSource, sourcePairs, sourcePublish,
1096  p_in, type, &position, bufferSize, fromOurSelf);
1097  send=receive;
1098  }
1099  }else{
1100 
1101  int oldPos=position;
1102  // Two index sets received
1103  unpackIndices(*receive, noRemoteSource, destPairs, destPublish,
1104  p_in, type, &position, bufferSize, fromOurSelf);
1105  if(!sendTwo)
1106  //unpack source entries again as destination entries
1107  position=oldPos;
1108 
1109  send = new RemoteIndexList();
1110  unpackIndices(*send, noRemoteDest, sourcePairs, sourcePublish,
1111  p_in, type, &position, bufferSize, fromOurSelf);
1112  }
1113 
1114  if(receive->empty() && send->empty()){
1115  if(send==receive){
1116  delete send;
1117  }else{
1118  delete send;
1119  delete receive;
1120  }
1121  }else{
1122  remoteIndices_.insert(std::make_pair(remoteProc,
1123  std::make_pair(send,receive)));
1124  }
1125  }
1126 
1127 
1128  template<typename T, typename A>
1129  template<bool ignorePublic>
1130  inline void RemoteIndices<T,A>::buildRemote(bool includeSelf)
1131  {
1132  // Processor configuration
1133  int rank, procs;
1134  MPI_Comm_rank(comm_, &rank);
1135  MPI_Comm_size(comm_, &procs);
1136 
1137  // number of local indices to publish
1138  // The indices of the destination will be send.
1139  int sourcePublish, destPublish;
1140 
1141  // Do we need to send two index sets?
1142  char sendTwo = (source_ != target_);
1143 
1144  if(procs==1 && !(sendTwo || includeSelf))
1145  // Nothing to communicate
1146  return;
1147 
1148  sourcePublish = (ignorePublic)? source_->size() : noPublic(*source_);
1149 
1150  if(sendTwo)
1151  destPublish = (ignorePublic)? target_->size() : noPublic(*target_);
1152  else
1153  // we only need to send one set of indices
1154  destPublish = 0;
1155 
1156  int maxPublish, publish=sourcePublish+destPublish;
1157 
1158  // Calucate maximum number of indices send
1159  MPI_Allreduce(&publish, &maxPublish, 1, MPI_INT, MPI_MAX, comm_);
1160 
1161  // allocate buffers
1162  typedef IndexPair<GlobalIndex,LocalIndex> PairType;
1163 
1164  PairType** destPairs;
1165  PairType** sourcePairs = new PairType*[sourcePublish>0?sourcePublish:1];
1166 
1167  if(sendTwo)
1168  destPairs = new PairType*[destPublish>0?destPublish:1];
1169  else
1170  destPairs=sourcePairs;
1171 
1172  char** buffer = new char*[2];
1173  int bufferSize;
1174  int position=0;
1175  int intSize;
1176  int charSize;
1177 
1178  // calculate buffer size
1179  MPI_Datatype type = MPITraits<PairType>::getType();
1180 
1181  MPI_Pack_size(maxPublish, type, comm_,
1182  &bufferSize);
1183  MPI_Pack_size(1, MPI_INT, comm_,
1184  &intSize);
1185  MPI_Pack_size(1, MPI_CHAR, comm_,
1186  &charSize);
1187  // Our message will contain the following:
1188  // a bool wether two index sets where sent
1189  // the size of the source and the dest indexset,
1190  // then the source and destination indices
1191  bufferSize += 2 * intSize + charSize;
1192 
1193  if(bufferSize<=0) bufferSize=1;
1194 
1195  buffer[0] = new char[bufferSize];
1196  buffer[1] = new char[bufferSize];
1197 
1198 
1199  // pack entries into buffer[0], p_out below!
1200  MPI_Pack(&sendTwo, 1, MPI_CHAR, buffer[0], bufferSize, &position,
1201  comm_);
1202 
1203  // The number of indices we send for each index set
1204  MPI_Pack(&sourcePublish, 1, MPI_INT, buffer[0], bufferSize, &position,
1205  comm_);
1206  MPI_Pack(&destPublish, 1, MPI_INT, buffer[0], bufferSize, &position,
1207  comm_);
1208 
1209  // Now pack the source indices and setup the destination pairs
1210  packEntries<ignorePublic>(sourcePairs, *source_, buffer[0], type,
1211  bufferSize, &position, sourcePublish);
1212  // If necessary send the dest indices and setup the source pairs
1213  if(sendTwo)
1214  packEntries<ignorePublic>(destPairs, *target_, buffer[0], type,
1215  bufferSize, &position, destPublish);
1216 
1217 
1218  // Update remote indices for ourself
1219  if(sendTwo|| includeSelf)
1220  unpackCreateRemote(buffer[0], sourcePairs, destPairs, rank, sourcePublish,
1221  destPublish, bufferSize, sendTwo, includeSelf);
1222 
1223  neighbourIds.erase(rank);
1224 
1225  if(neighbourIds.size()==0)
1226  {
1227  Dune::dvverb<<rank<<": Sending messages in a ring"<<std::endl;
1228  // send messages in ring
1229  for(int proc=1; proc<procs; proc++){
1230  // pointers to the current input and output buffers
1231  char* p_out = buffer[1-(proc%2)];
1232  char* p_in = buffer[proc%2];
1233 
1234  MPI_Status status;
1235  if(rank%2==0){
1236  MPI_Ssend(p_out, bufferSize, MPI_PACKED, (rank+1)%procs,
1237  commTag_, comm_);
1238  MPI_Recv(p_in, bufferSize, MPI_PACKED, (rank+procs-1)%procs,
1239  commTag_, comm_, &status);
1240  }else{
1241  MPI_Recv(p_in, bufferSize, MPI_PACKED, (rank+procs-1)%procs,
1242  commTag_, comm_, &status);
1243  MPI_Ssend(p_out, bufferSize, MPI_PACKED, (rank+1)%procs,
1244  commTag_, comm_);
1245  }
1246 
1247 
1248  // The process these indices are from
1249  int remoteProc = (rank+procs-proc)%procs;
1250 
1251  unpackCreateRemote(p_in, sourcePairs, destPairs, remoteProc, sourcePublish,
1252  destPublish, bufferSize, sendTwo);
1253 
1254  }
1255 
1256  }
1257  else
1258  {
1259  MPI_Request* requests=new MPI_Request[neighbourIds.size()];
1260  MPI_Request* req=requests;
1261 
1262  typedef typename std::set<int>::size_type size_type;
1263  size_type noNeighbours=neighbourIds.size();
1264 
1265  // setup sends
1266  for(std::set<int>::iterator neighbour=neighbourIds.begin();
1267  neighbour!= neighbourIds.end(); ++neighbour){
1268  // Only send the information to the neighbouring processors
1269  MPI_Issend(buffer[0], position , MPI_PACKED, *neighbour, commTag_, comm_, req++);
1270  }
1271 
1272  //Test for received messages
1273 
1274  for(size_type received=0; received <noNeighbours; ++received)
1275  {
1276  MPI_Status status;
1277  // probe for next message
1278  MPI_Probe(MPI_ANY_SOURCE, commTag_, comm_, &status);
1279  int remoteProc=status.MPI_SOURCE;
1280  int size;
1281  MPI_Get_count(&status, MPI_PACKED, &size);
1282  // receive message
1283  MPI_Recv(buffer[1], size, MPI_PACKED, remoteProc,
1284  commTag_, comm_, &status);
1285 
1286  unpackCreateRemote(buffer[1], sourcePairs, destPairs, remoteProc, sourcePublish,
1287  destPublish, bufferSize, sendTwo);
1288  }
1289  // wait for completion of pending requests
1290  MPI_Status* statuses = new MPI_Status[neighbourIds.size()];
1291 
1292  if(MPI_ERR_IN_STATUS==MPI_Waitall(neighbourIds.size(), requests, statuses)){
1293  for(size_type i=0; i < neighbourIds.size(); ++i)
1294  if(statuses[i].MPI_ERROR!=MPI_SUCCESS){
1295  std::cerr<<rank<<": MPI_Error occurred while receiving message."<<std::endl;
1296  MPI_Abort(comm_, 999);
1297  }
1298  }
1299  delete[] requests;
1300  delete[] statuses;
1301  }
1302 
1303 
1304  // delete allocated memory
1305  if(destPairs!=sourcePairs)
1306  delete[] destPairs;
1307 
1308  delete[] sourcePairs;
1309  delete[] buffer[0];
1310  delete[] buffer[1];
1311  delete[] buffer;
1312  }
1313 
1314  template<typename T, typename A>
1315  inline void RemoteIndices<T,A>::unpackIndices(RemoteIndexList& remote,
1316  int remoteEntries,
1317  PairType** local,
1318  int localEntries,
1319  char* p_in,
1320  MPI_Datatype type,
1321  int* position,
1322  int bufferSize,
1323  bool fromOurSelf)
1324  {
1325  if(remoteEntries==0)
1326  return;
1327 
1328  PairType index(-1);
1329  MPI_Unpack(p_in, bufferSize, position, &index, 1,
1330  type, comm_);
1331  GlobalIndex oldGlobal=index.global();
1332  int n_in=0, localIndex=0;
1333 
1334  //Check if we know the global index
1335  while(localIndex<localEntries){
1336  if(local[localIndex]->global()==index.global()){
1337  int oldLocalIndex=localIndex;
1338 
1339  while(localIndex<localEntries &&
1340  local[localIndex]->global()==index.global()){
1341  if(!fromOurSelf || index.local().attribute() !=
1342  local[localIndex]->local().attribute())
1343  // if index is from us it has to have a different attribute
1344  remote.push_back(RemoteIndex(index.local().attribute(),
1345  local[localIndex]));
1346  localIndex++;
1347  }
1348 
1349  // unpack next remote index
1350  if((++n_in) < remoteEntries){
1351  MPI_Unpack(p_in, bufferSize, position, &index, 1,
1352  type, comm_);
1353  if(index.global()==oldGlobal)
1354  // Restart comparison for the same global indices
1355  localIndex=oldLocalIndex;
1356  else
1357  oldGlobal=index.global();
1358  }else{
1359  // No more received indices
1360  break;
1361  }
1362  continue;
1363  }
1364 
1365  if (local[localIndex]->global()<index.global()){
1366  // compare with next entry in our list
1367  ++localIndex;
1368  }else{
1369  // We do not know the index, unpack next
1370  if((++n_in) < remoteEntries){
1371  MPI_Unpack(p_in, bufferSize, position, &index, 1,
1372  type, comm_);
1373  oldGlobal=index.global();
1374  }else
1375  // No more received indices
1376  break;
1377  }
1378  }
1379 
1380  // Unpack the other received indices without doing anything
1381  while(++n_in < remoteEntries)
1382  MPI_Unpack(p_in, bufferSize, position, &index, 1,
1383  type, comm_);
1384  }
1385 
1386 
1387  template<typename T, typename A>
1388  inline void RemoteIndices<T,A>::unpackIndices(RemoteIndexList& send,
1389  RemoteIndexList& receive,
1390  int remoteEntries,
1391  PairType** localSource,
1392  int localSourceEntries,
1393  PairType** localDest,
1394  int localDestEntries,
1395  char* p_in,
1396  MPI_Datatype type,
1397  int* position,
1398  int bufferSize)
1399  {
1400  int n_in=0, sourceIndex=0, destIndex=0;
1401 
1402  //Check if we know the global index
1403  while(n_in<remoteEntries && (sourceIndex<localSourceEntries || destIndex<localDestEntries)){
1404  // Unpack next index
1405  PairType index;
1406  MPI_Unpack(p_in, bufferSize, position, &index, 1,
1407  type, comm_);
1408  n_in++;
1409 
1410  // Advance until global index in localSource and localDest are >= than the one in the unpacked index
1411  while(sourceIndex<localSourceEntries && localSource[sourceIndex]->global()<index.global())
1412  sourceIndex++;
1413 
1414  while(destIndex<localDestEntries && localDest[destIndex]->global()<index.global())
1415  destIndex++;
1416 
1417  // Add a remote index if we found the global index.
1418  if(sourceIndex<localSourceEntries && localSource[sourceIndex]->global()==index.global())
1419  send.push_back(RemoteIndex(index.local().attribute(),
1420  localSource[sourceIndex]));
1421 
1422  if(destIndex < localDestEntries && localDest[destIndex]->global() == index.global())
1423  receive.push_back(RemoteIndex(index.local().attribute(),
1424  localDest[sourceIndex]));
1425  }
1426 
1427  }
1428 
1429  template<typename T, typename A>
1431  {
1432  typedef typename RemoteIndexMap::iterator Iterator;
1433  Iterator lend = remoteIndices_.end();
1434  for(Iterator lists=remoteIndices_.begin(); lists != lend; ++lists){
1435  if(lists->second.first==lists->second.second){
1436  // there is only one remote index list.
1437  delete lists->second.first;
1438  }else{
1439  delete lists->second.first;
1440  delete lists->second.second;
1441  }
1442  }
1443  remoteIndices_.clear();
1444  firstBuild=true;
1445  }
1446 
1447  template<typename T, typename A>
1449  {
1450  return remoteIndices_.size();
1451  }
1452 
1453  template<typename T, typename A>
1454  template<bool ignorePublic>
1456  {
1457  // Test wether a rebuild is Needed.
1458  if(firstBuild ||
1459  ignorePublic!=publicIgnored || !
1460  isSynced()){
1461  free();
1462 
1463  buildRemote<ignorePublic>(includeSelf);
1464 
1465  sourceSeqNo_ = source_->seqNo();
1466  destSeqNo_ = target_->seqNo();
1467  firstBuild=false;
1468  publicIgnored=ignorePublic;
1469  }
1470 
1471 
1472  }
1473 
1474  template<typename T, typename A>
1475  inline bool RemoteIndices<T,A>::isSynced() const
1476  {
1477  return sourceSeqNo_==source_->seqNo() && destSeqNo_ ==target_->seqNo();
1478  }
1479 
1480  template<typename T, typename A>
1481  template<bool mode, bool send>
1483  {
1484 
1485  // The user are on their own now!
1486  // We assume they know what they are doing and just set the
1487  // remote indices to synced status.
1488  sourceSeqNo_ = source_->seqNo();
1489  destSeqNo_ = target_->seqNo();
1490 
1491  typename RemoteIndexMap::iterator found = remoteIndices_.find(process);
1492 
1493  if(found == remoteIndices_.end())
1494  {
1495  if(source_ != target_)
1496  remoteIndices_.insert(std::make_pair(process,
1497  std::make_pair(new RemoteIndexList(),
1498  new RemoteIndexList())));
1499  else{
1500  RemoteIndexList* rlist = new RemoteIndexList();
1501  remoteIndices_.insert(std::make_pair(process,
1502  std::make_pair(rlist, rlist)));
1503 
1504  found = remoteIndices_.find(process);
1505  }
1506  }
1507 
1508  firstBuild = false;
1509 
1510  if(send)
1511  return RemoteIndexListModifier<T,A,mode>(*source_, *(found->second.first));
1512  else
1513  return RemoteIndexListModifier<T,A,mode>(*target_, *(found->second.second));
1514  }
1515 
1516  template<typename T, typename A>
1517  inline typename RemoteIndices<T,A>::const_iterator
1519  {
1520  return remoteIndices_.find(proc);
1521  }
1522 
1523  template<typename T, typename A>
1524  inline typename RemoteIndices<T,A>::const_iterator
1526  {
1527  return remoteIndices_.begin();
1528  }
1529 
1530  template<typename T, typename A>
1531  inline typename RemoteIndices<T,A>::const_iterator
1533  {
1534  return remoteIndices_.end();
1535  }
1536 
1537 
1538  template<typename T, typename A>
1540  {
1541  if(neighbours()!=ri.neighbours())
1542  return false;
1543 
1544  typedef RemoteIndexList RList;
1545  typedef typename std::map<int,std::pair<RList*,RList*> >::const_iterator const_iterator;
1546 
1547  const const_iterator rend = remoteIndices_.end();
1548 
1549  for(const_iterator rindex = remoteIndices_.begin(), rindex1=ri.remoteIndices_.begin(); rindex!=rend; ++rindex, ++rindex1){
1550  if(rindex->first != rindex1->first)
1551  return false;
1552  if(*(rindex->second.first) != *(rindex1->second.first))
1553  return false;
1554  if(*(rindex->second.second) != *(rindex1->second.second))
1555  return false;
1556  }
1557  return true;
1558  }
1559 
1560  template<class T, class A, bool mode>
1562  RemoteIndexList& rList)
1563  : rList_(&rList), indexSet_(&indexSet), glist_(new GlobalList()), iter_(rList.beginModify()), end_(rList.end()), first_(true)
1564  {
1565  if(MODIFYINDEXSET){
1566  assert(indexSet_);
1567  for(ConstIterator iter=iter_; iter != end_; ++iter)
1568  glist_->push_back(iter->localIndexPair().global());
1569  giter_ = glist_->beginModify();
1570  }
1571  }
1572 
1573  template<typename T, typename A, bool mode>
1575  : rList_(other.rList_), indexSet_(other.indexSet_),
1576  glist_(other.glist_), iter_(other.iter_), giter_(other.giter_), end_(other.end_),
1577  first_(other.first_), last_(other.last_)
1578  {}
1579 
1580  template<typename T, typename A, bool mode>
1582  {
1583  if(MODIFYINDEXSET){
1584  // repair pointers to local index set.
1585 #ifdef DUNE_ISTL_WITH_CHECKING
1586  if(indexSet_->state()!=GROUND)
1587  DUNE_THROW(InvalidIndexSetState, "Index has to be in ground mode for repairing pointers to indices");
1588 #endif
1589  typedef typename ParallelIndexSet::const_iterator IndexIterator;
1590  typedef typename GlobalList::const_iterator GlobalIterator;
1591  typedef typename RemoteIndexList::iterator Iterator;
1592  GlobalIterator giter = glist_->begin();
1593  IndexIterator index = indexSet_->begin();
1594 
1595  for(Iterator iter=rList_->begin(); iter != end_; ++iter){
1596  while(index->global()<*giter){
1597  ++index;
1598 #ifdef DUNE_ISTL_WITH_CHECKING
1599  if(index == indexSet_.end())
1600  DUNE_THROW(InvalidPosition, "No such global index in set!");
1601 #endif
1602  }
1603 
1604 #ifdef DUNE_ISTL_WITH_CHECKING
1605  if(index->global() != *giter)
1606  DUNE_THROW(InvalidPosition, "No such global index in set!");
1607 #endif
1608  iter->localIndex_ = &(*index);
1609  }
1610  }
1611  }
1612 
1613  template<typename T, typename A, bool mode>
1615  {
1616  dune_static_assert(!mode,"Not allowed if the mode indicates that new indices"
1617  "might be added to the underlying index set. Use "
1618  "insert(const RemoteIndex&, const GlobalIndex&) instead");
1619 
1620 #ifdef DUNE_ISTL_WITH_CHECKING
1621  if(!first_ && index.localIndexPair().global()<last_)
1622  DUNE_THROW(InvalidPosition, "Modifcation of remote indices have to occur with ascending global index.");
1623 #endif
1624  // Move to the correct position
1625  while(iter_ != end_ && iter_->localIndexPair().global() < index.localIndexPair().global()){
1626  ++iter_;
1627  }
1628 
1629  // No duplicate entries allowed
1630  assert(iter_==end_ || iter_->localIndexPair().global() != index.localIndexPair().global());
1631  iter_.insert(index);
1632  last_ = index.localIndexPair().global();
1633  first_ = false;
1634  }
1635 
1636  template<typename T, typename A, bool mode>
1638  {
1639  dune_static_assert(mode,"Not allowed if the mode indicates that no new indices"
1640  "might be added to the underlying index set. Use "
1641  "insert(const RemoteIndex&) instead");
1642 #ifdef DUNE_ISTL_WITH_CHECKING
1643  if(!first_ && global<last_)
1644  DUNE_THROW(InvalidPosition, "Modification of remote indices have to occur with ascending global index.");
1645 #endif
1646  // Move to the correct position
1647  while(iter_ != end_ && *giter_ < global){
1648  ++giter_;
1649  ++iter_;
1650  }
1651 
1652  // No duplicate entries allowed
1653  assert(iter_->localIndexPair().global() != global);
1654  iter_.insert(index);
1655  giter_.insert(global);
1656 
1657  last_ = global;
1658  first_ = false;
1659  }
1660 
1661  template<typename T, typename A, bool mode>
1663  {
1664 #ifdef DUNE_ISTL_WITH_CHECKING
1665  if(!first_ && global<last_)
1666  DUNE_THROW(InvalidPosition, "Modifcation of remote indices have to occur with ascending global index.");
1667 #endif
1668 
1669  bool found= false;
1670 
1671  if(MODIFYINDEXSET){
1672  // Move to the correct position
1673  while(iter_!=end_ && *giter_< global){
1674  ++giter_;
1675  ++iter_;
1676  }
1677  if(*giter_ == global){
1678  giter_.remove();
1679  iter_.remove();
1680  found=true;
1681  }
1682  }else{
1683  while(iter_!=end_ && iter_->localIndexPair().global() < global)
1684  ++iter_;
1685 
1686  if(iter_->localIndexPair().global()==global){
1687  iter_.remove();
1688  found = true;
1689  }
1690  }
1691 
1692  last_ = global;
1693  first_ = false;
1694  return found;
1695  }
1696 
1697  template<typename T, typename A>
1698  template<bool send>
1700  {
1701  return CollectiveIterator<T,A>(remoteIndices_, send);
1702  }
1703 
1704  template<typename T, typename A>
1705  inline MPI_Comm RemoteIndices<T,A>::communicator() const
1706  {
1707  return comm_;
1708 
1709  }
1710 
1711  template<typename T, typename A>
1713  {
1714  typedef typename RemoteIndexMap::const_iterator const_iterator;
1715 
1716  const const_iterator end=pmap.end();
1717  for(const_iterator process=pmap.begin(); process != end; ++process){
1718  const RemoteIndexList* list = send? process->second.first : process->second.second;
1719  typedef typename RemoteIndexList::const_iterator iterator;
1720  map_.insert(std::make_pair(process->first,
1721  std::pair<iterator, const iterator>(list->begin(), list->end())));
1722  }
1723  }
1724 
1725  template<typename T, typename A>
1726  inline void CollectiveIterator<T,A>::advance(const GlobalIndex& index)
1727  {
1728  typedef typename Map::iterator iterator;
1729  typedef typename Map::const_iterator const_iterator;
1730  const const_iterator end = map_.end();
1731 
1732  for(iterator iter = map_.begin(); iter != end;){
1733  // Step the iterator until we are >= index
1734  typename RemoteIndexList::const_iterator current = iter->second.first;
1735  typename RemoteIndexList::const_iterator rend = iter->second.second;
1736  RemoteIndex remoteIndex;
1737  if(current != rend)
1738  remoteIndex = *current;
1739 
1740  while(iter->second.first!=iter->second.second && iter->second.first->localIndexPair().global()<index)
1741  ++(iter->second.first);
1742 
1743  // erase from the map if there are no more entries.
1744  if(iter->second.first == iter->second.second)
1745  map_.erase(iter++);
1746  else{
1747  ++iter;
1748  }
1749  }
1750  index_=index;
1751  noattribute=true;
1752  }
1753 
1754  template<typename T, typename A>
1755  inline void CollectiveIterator<T,A>::advance(const GlobalIndex& index,
1756  const Attribute& attribute)
1757  {
1758  typedef typename Map::iterator iterator;
1759  typedef typename Map::const_iterator const_iterator;
1760  const const_iterator end = map_.end();
1761 
1762  for(iterator iter = map_.begin(); iter != end;){
1763  // Step the iterator until we are >= index
1764  typename RemoteIndexList::const_iterator current = iter->second.first;
1765  typename RemoteIndexList::const_iterator rend = iter->second.second;
1766  RemoteIndex remoteIndex;
1767  if(current != rend)
1768  remoteIndex = *current;
1769 
1770  // Move to global index or bigger
1771  while(iter->second.first!=iter->second.second && iter->second.first->localIndexPair().global()<index)
1772  ++(iter->second.first);
1773 
1774  // move to attribute or bigger
1775  while(iter->second.first!=iter->second.second
1776  && iter->second.first->localIndexPair().global()==index
1777  && iter->second.first->localIndexPair().local().attribute()<attribute)
1778  ++(iter->second.first);
1779 
1780  // erase from the map if there are no more entries.
1781  if(iter->second.first == iter->second.second)
1782  map_.erase(iter++);
1783  else{
1784  ++iter;
1785  }
1786  }
1787  index_=index;
1788  attribute_=attribute;
1789  noattribute=false;
1790  }
1791 
1792  template<typename T, typename A>
1794  {
1795  typedef typename Map::iterator iterator;
1796  typedef typename Map::const_iterator const_iterator;
1797  const const_iterator end = map_.end();
1798 
1799  for(iterator iter = map_.begin(); iter != end;){
1800  // Step the iterator until we are >= index
1801  typename RemoteIndexList::const_iterator current = iter->second.first;
1802  typename RemoteIndexList::const_iterator rend = iter->second.second;
1803 
1804  // move all iterators pointing to the current global index to next value
1805  if(iter->second.first->localIndexPair().global()==index_ &&
1806  (noattribute || iter->second.first->localIndexPair().local().attribute() == attribute_))
1807  ++(iter->second.first);
1808 
1809  // erase from the map if there are no more entries.
1810  if(iter->second.first == iter->second.second)
1811  map_.erase(iter++);
1812  else{
1813  ++iter;
1814  }
1815  }
1816  return *this;
1817  }
1818 
1819  template<typename T, typename A>
1821  {
1822  return map_.empty();
1823  }
1824 
1825  template<typename T, typename A>
1826  inline typename CollectiveIterator<T,A>::iterator
1828  {
1829  if(noattribute)
1830  return iterator(map_.begin(), map_.end(), index_);
1831  else
1832  return iterator(map_.begin(), map_.end(), index_,
1833  attribute_);
1834  }
1835 
1836  template<typename T, typename A>
1837  inline typename CollectiveIterator<T,A>::iterator
1839  {
1840  return iterator(map_.end(), map_.end(), index_);
1841  }
1842 
1843  template<typename TG, typename TA>
1844  inline std::ostream& operator<<(std::ostream& os, const RemoteIndex<TG,TA>& index)
1845  {
1846  os<<"[global="<<index.localIndexPair().global()<<", remote attribute="<<index.attribute()<<" local attribute="<<index.localIndexPair().local().attribute()<<"]";
1847  return os;
1848  }
1849 
1850  template<typename T, typename A>
1851  inline std::ostream& operator<<(std::ostream& os, const RemoteIndices<T,A>& indices)
1852  {
1853  int rank;
1854  MPI_Comm_rank(indices.comm_, &rank);
1855 
1856  typedef typename RemoteIndices<T,A>::RemoteIndexList RList;
1857  typedef typename std::map<int,std::pair<RList*,RList*> >::const_iterator const_iterator;
1858 
1859  const const_iterator rend = indices.remoteIndices_.end();
1860 
1861  for(const_iterator rindex = indices.remoteIndices_.begin(); rindex!=rend; ++rindex){
1862  os<<rank<<": Prozess "<<rindex->first<<":";
1863 
1864  if(!rindex->second.first->empty()){
1865  os<<" send:";
1866 
1867  const typename RList::const_iterator send= rindex->second.first->end();
1868 
1869  for(typename RList::const_iterator index = rindex->second.first->begin();
1870  index != send; ++index)
1871  os<<*index<<" ";
1872  os<<std::endl;
1873  }
1874  if(!rindex->second.second->empty()){
1875  os<<rank<<": Prozess "<<rindex->first<<": "<<"receive: ";
1876 
1877  const typename RList::const_iterator rend= rindex->second.second->end();
1878 
1879  for(typename RList::const_iterator index = rindex->second.second->begin();
1880  index != rend; ++index)
1881  os<<*index<<" ";
1882  }
1883  os<<std::endl<<std::flush;
1884  }
1885  return os;
1886  }
1888 }
1889 
1890 #endif
1891 #endif