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