dune-istl  2.8.0
repartition.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_REPARTITION_HH
4 #define DUNE_ISTL_REPARTITION_HH
5 
6 #include <cassert>
7 #include <map>
8 #include <utility>
9 #include <cmath>
10 
11 #if HAVE_PARMETIS
12 // Explicitly use C linkage as scotch does not extern "C" in its headers.
13 // Works because ParMETIS/METIS checks whether compiler is C++ and otherwise
14 // does not use extern "C". Therfore no nested extern "C" will be created
15 extern "C"
16 {
17 #include <parmetis.h>
18 }
19 #endif
20 
21 #include <dune/common/timer.hh>
22 #include <dune/common/enumset.hh>
23 #include <dune/common/stdstreams.hh>
24 #include <dune/common/parallel/mpitraits.hh>
25 #include <dune/common/parallel/communicator.hh>
26 #include <dune/common/parallel/indexset.hh>
27 #include <dune/common/parallel/indicessyncer.hh>
28 #include <dune/common/parallel/remoteindices.hh>
29 #include <dune/common/rangeutilities.hh>
30 
32 #include <dune/istl/paamg/graph.hh>
33 
42 namespace Dune
43 {
44  namespace Metis
45  {
46  // Explicitly specify a real_t and idx_t for older (Par)METIS versions that do not
47  // provide these typedefs
48 #if HAVE_PARMETIS && defined(REALTYPEWIDTH)
49  using real_t = ::real_t;
50 #else
51  using real_t = float;
52 #endif
53 
54 #if HAVE_PARMETIS && defined(IDXTYPEWIDTH)
55  using idx_t = ::idx_t;
56 #elif HAVE_PARMETIS && defined(HAVE_SCOTCH_NUM_TYPE)
57  using idx_t = SCOTCH_Num;
58 #elif HAVE_PARMETIS
59  using idx_t = int;
60 #else
61  using idx_t = std::size_t;
62 #endif
63  }
64 
65 
66 #if HAVE_MPI
80  template<class G, class T1, class T2>
82  {
84  typedef typename IndexSet::LocalIndex::Attribute Attribute;
85 
86  IndexSet& indexSet = oocomm.indexSet();
88 
89  std::size_t sum=0, needed = graph.noVertices()-indexSet.size();
90  std::vector<std::size_t> neededall(oocomm.communicator().size(), 0);
91 
92  MPI_Allgather(&needed, 1, MPITraits<std::size_t>::getType() , &(neededall[0]), 1, MPITraits<std::size_t>::getType(), oocomm.communicator());
93  for(int i=0; i<oocomm.communicator().size(); ++i)
94  sum=sum+neededall[i]; // MAke this for generic
95 
96  if(sum==0)
97  // Nothing to do
98  return;
99 
100  //Compute Maximum Global Index
101  T1 maxgi=0;
102  auto end = indexSet.end();
103  for(auto it = indexSet.begin(); it != end; ++it)
104  maxgi=std::max(maxgi,it->global());
105 
106  //Process p creates global indices consecutively
107  //starting atmaxgi+\sum_{i=1}^p neededall[i]
108  // All created indices are owned by the process
109  maxgi=oocomm.communicator().max(maxgi);
110  ++maxgi; //Sart with the next free index.
111 
112  for(int i=0; i<oocomm.communicator().rank(); ++i)
113  maxgi=maxgi+neededall[i]; // TODO: make this more generic
114 
115  // Store the global index information for repairing the remote index information
116  std::map<int,SLList<std::pair<T1,Attribute> > > globalIndices;
117  storeGlobalIndicesOfRemoteIndices(globalIndices, oocomm.remoteIndices());
118  indexSet.beginResize();
119 
120  for(auto vertex = graph.begin(), vend=graph.end(); vertex != vend; ++vertex) {
121  const typename IndexSet::IndexPair* pair=lookup.pair(*vertex);
122  if(pair==0) {
123  // No index yet, add new one
124  indexSet.add(maxgi, typename IndexSet::LocalIndex(*vertex, OwnerOverlapCopyAttributeSet::owner, false));
125  ++maxgi;
126  }
127  }
128 
129  indexSet.endResize();
130 
131  repairLocalIndexPointers(globalIndices, oocomm.remoteIndices(), indexSet);
132 
133  oocomm.freeGlobalLookup();
134  oocomm.buildGlobalLookup();
135 #ifdef DEBUG_REPART
136  std::cout<<"Holes are filled!"<<std::endl;
137  std::cout<<oocomm.communicator().rank()<<": "<<oocomm.indexSet()<<std::endl;
138 #endif
139  }
140 
141  namespace
142  {
143 
144  class ParmetisDuneIndexMap
145  {
146  public:
147  template<class Graph, class OOComm>
148  ParmetisDuneIndexMap(const Graph& graph, const OOComm& com);
149  int toParmetis(int i) const
150  {
151  return duneToParmetis[i];
152  }
153  int toLocalParmetis(int i) const
154  {
155  return duneToParmetis[i]-base_;
156  }
157  int operator[](int i) const
158  {
159  return duneToParmetis[i];
160  }
161  int toDune(int i) const
162  {
163  return parmetisToDune[i];
164  }
165  std::vector<int>::size_type numOfOwnVtx() const
166  {
167  return parmetisToDune.size();
168  }
169  Metis::idx_t* vtxDist()
170  {
171  return &vtxDist_[0];
172  }
174  private:
175  int base_;
176  std::vector<int> duneToParmetis;
177  std::vector<int> parmetisToDune;
178  // range of vertices for processor i: vtxdist[i] to vtxdist[i+1] (parmetis global)
179  std::vector<Metis::idx_t> vtxDist_;
180  };
181 
182  template<class G, class OOComm>
183  ParmetisDuneIndexMap::ParmetisDuneIndexMap(const G& graph, const OOComm& oocomm)
184  : duneToParmetis(graph.noVertices(), -1), vtxDist_(oocomm.communicator().size()+1)
185  {
186  int npes=oocomm.communicator().size(), mype=oocomm.communicator().rank();
187 
188  typedef typename OOComm::OwnerSet OwnerSet;
189 
190  int numOfOwnVtx=0;
191  auto end = oocomm.indexSet().end();
192  for(auto index = oocomm.indexSet().begin(); index != end; ++index) {
193  if (OwnerSet::contains(index->local().attribute())) {
194  numOfOwnVtx++;
195  }
196  }
197  parmetisToDune.resize(numOfOwnVtx);
198  std::vector<int> globalNumOfVtx(npes);
199  // make this number available to all processes
200  MPI_Allgather(&numOfOwnVtx, 1, MPI_INT, &(globalNumOfVtx[0]), 1, MPI_INT, oocomm.communicator());
201 
202  int base=0;
203  vtxDist_[0] = 0;
204  for(int i=0; i<npes; i++) {
205  if (i<mype) {
206  base += globalNumOfVtx[i];
207  }
208  vtxDist_[i+1] = vtxDist_[i] + globalNumOfVtx[i];
209  }
210  globalOwnerVertices=vtxDist_[npes];
211  base_=base;
212 
213 #ifdef DEBUG_REPART
214  std::cout << oocomm.communicator().rank()<<" vtxDist: ";
215  for(int i=0; i<= npes; ++i)
216  std::cout << vtxDist_[i]<<" ";
217  std::cout<<std::endl;
218 #endif
219 
220  // Traverse the graph and assign a new consecutive number/index
221  // starting by "base" to all owner vertices.
222  // The new index is used as the ParMETIS global index and is
223  // stored in the vector "duneToParmetis"
224  auto vend = graph.end();
225  for(auto vertex = graph.begin(); vertex != vend; ++vertex) {
226  const typename OOComm::ParallelIndexSet::IndexPair* index=oocomm.globalLookup().pair(*vertex);
227  assert(index);
228  if (OwnerSet::contains(index->local().attribute())) {
229  // assign and count the index
230  parmetisToDune[base-base_]=index->local();
231  duneToParmetis[index->local()] = base++;
232  }
233  }
234 
235  // At this point, every process knows the ParMETIS global index
236  // of it's owner vertices. The next step is to get the
237  // ParMETIS global index of the overlap vertices from the
238  // associated processes. To do this, the Dune::Interface class
239  // is used.
240 #ifdef DEBUG_REPART
241  std::cout <<oocomm.communicator().rank()<<": before ";
242  for(std::size_t i=0; i<duneToParmetis.size(); ++i)
243  std::cout<<duneToParmetis[i]<<" ";
244  std::cout<<std::endl;
245 #endif
246  oocomm.copyOwnerToAll(duneToParmetis,duneToParmetis);
247 #ifdef DEBUG_REPART
248  std::cout <<oocomm.communicator().rank()<<": after ";
249  for(std::size_t i=0; i<duneToParmetis.size(); ++i)
250  std::cout<<duneToParmetis[i]<<" ";
251  std::cout<<std::endl;
252 #endif
253  }
254  }
255 
257  : public Interface
258  {
259  void setCommunicator(MPI_Comm comm)
260  {
261  communicator_=comm;
262  }
263  template<class Flags,class IS>
264  void buildSendInterface(const std::vector<int>& toPart, const IS& idxset)
265  {
266  std::map<int,int> sizes;
267 
268  for(auto i=idxset.begin(), end=idxset.end(); i!=end; ++i)
269  if(Flags::contains(i->local().attribute()))
270  ++sizes[toPart[i->local()]];
271 
272  // Allocate the necessary space
273  for(auto i=sizes.begin(), end=sizes.end(); i!=end; ++i)
274  interfaces()[i->first].first.reserve(i->second);
275 
276  //Insert the interface information
277  for(auto i=idxset.begin(), end=idxset.end(); i!=end; ++i)
278  if(Flags::contains(i->local().attribute()))
279  interfaces()[toPart[i->local()]].first.add(i->local());
280  }
281 
282  void reserveSpaceForReceiveInterface(int proc, int size)
283  {
284  interfaces()[proc].second.reserve(size);
285  }
286  void addReceiveIndex(int proc, std::size_t idx)
287  {
288  interfaces()[proc].second.add(idx);
289  }
290  template<typename TG>
291  void buildReceiveInterface(std::vector<std::pair<TG,int> >& indices)
292  {
293  std::size_t i=0;
294  for(auto idx=indices.begin(); idx!= indices.end(); ++idx) {
295  interfaces()[idx->second].second.add(i++);
296  }
297  }
298 
300  {}
301 
302  };
303 
304  namespace
305  {
315  template<class GI>
316  void createSendBuf(std::vector<GI>& ownerVec, std::set<GI>& overlapVec, std::set<int>& neighbors, char *sendBuf, int buffersize, MPI_Comm comm) {
317  // Pack owner vertices
318  std::size_t s=ownerVec.size();
319  int pos=0;
320  if(s==0)
321  ownerVec.resize(1); // otherwise would read beyond the memory bound
322  MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
323  MPI_Pack(&(ownerVec[0]), s, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
324  s = overlapVec.size();
325  MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
326  for(auto i=overlapVec.begin(), end= overlapVec.end(); i != end; ++i)
327  MPI_Pack(const_cast<GI*>(&(*i)), 1, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
328 
329  s=neighbors.size();
330  MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
331 
332  for(auto i=neighbors.begin(), end= neighbors.end(); i != end; ++i)
333  MPI_Pack(const_cast<int*>(&(*i)), 1, MPI_INT, sendBuf, buffersize, &pos, comm);
334  }
343  template<class GI>
344  void saveRecvBuf(char *recvBuf, int bufferSize, std::vector<std::pair<GI,int> >& ownerVec,
345  std::set<GI>& overlapVec, std::set<int>& neighbors, RedistributeInterface& inf, int from, MPI_Comm comm) {
346  std::size_t size;
347  int pos=0;
348  // unpack owner vertices
349  MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
350  inf.reserveSpaceForReceiveInterface(from, size);
351  ownerVec.reserve(ownerVec.size()+size);
352  for(; size!=0; --size) {
353  GI gi;
354  MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
355  ownerVec.push_back(std::make_pair(gi,from));
356  }
357  // unpack overlap vertices
358  MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
359  typename std::set<GI>::iterator ipos = overlapVec.begin();
360  Dune::dverb << "unpacking "<<size<<" overlap"<<std::endl;
361  for(; size!=0; --size) {
362  GI gi;
363  MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
364  ipos=overlapVec.insert(ipos, gi);
365  }
366  //unpack neighbors
367  MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
368  Dune::dverb << "unpacking "<<size<<" neighbors"<<std::endl;
369  typename std::set<int>::iterator npos = neighbors.begin();
370  for(; size!=0; --size) {
371  int n;
372  MPI_Unpack(recvBuf, bufferSize, &pos, &n, 1, MPI_INT, comm);
373  npos=neighbors.insert(npos, n);
374  }
375  }
376 
390  template<typename T>
391  void getDomain(const MPI_Comm& comm, T *part, int numOfOwnVtx, int nparts, int *myDomain, std::vector<int> &domainMapping) {
392  int npes, mype;
393  MPI_Comm_size(comm, &npes);
394  MPI_Comm_rank(comm, &mype);
395  MPI_Status status;
396 
397  *myDomain = -1;
398  int i=0;
399  int j=0;
400 
401  std::vector<int> domain(nparts, 0);
402  std::vector<int> assigned(npes, 0);
403  // init domain Mapping
404  domainMapping.assign(domainMapping.size(), -1);
405 
406  // count the occurrence of domains
407  for (i=0; i<numOfOwnVtx; i++) {
408  domain[part[i]]++;
409  }
410 
411  std::vector<int> domainMatrix(npes * nparts, -1);
412 
413  // init buffer with the own domain
414  int *buf = new int[nparts];
415  for (i=0; i<nparts; i++) {
416  buf[i] = domain[i];
417  domainMatrix[mype*nparts+i] = domain[i];
418  }
419  int pe=0;
420  int src = (mype-1+npes)%npes;
421  int dest = (mype+1)%npes;
422  // ring communication, we need n-1 communications for n processors
423  for (i=0; i<npes-1; i++) {
424  MPI_Sendrecv_replace(buf, nparts, MPI_INT, dest, 0, src, 0, comm, &status);
425  // pe is the process of the actual received buffer
426  pe = ((mype-1-i)+npes)%npes;
427  for(j=0; j<nparts; j++) {
428  // save the values to the domain matrix
429  domainMatrix[pe*nparts+j] = buf[j];
430  }
431  }
432  delete[] buf;
433 
434  // Start the domain calculation.
435  // The process which contains the maximum number of vertices of a
436  // particular domain is selected to choose it's favorate domain
437  int maxOccurance = 0;
438  pe = -1;
439  std::set<std::size_t> unassigned;
440 
441  for(i=0; i<nparts; i++) {
442  for(j=0; j<npes; j++) {
443  // process has no domain assigned
444  if (assigned[j]==0) {
445  if (maxOccurance < domainMatrix[j*nparts+i]) {
446  maxOccurance = domainMatrix[j*nparts+i];
447  pe = j;
448  }
449  }
450 
451  }
452  if (pe!=-1) {
453  // process got a domain, ...
454  domainMapping[i] = pe;
455  // ...mark as assigned
456  assigned[pe] = 1;
457  if (pe==mype) {
458  *myDomain = i;
459  }
460  pe = -1;
461  }
462  else
463  {
464  unassigned.insert(i);
465  }
466  maxOccurance = 0;
467  }
468 
469  typename std::vector<int>::iterator next_free = assigned.begin();
470 
471  for(auto udomain = unassigned.begin(),
472  end = unassigned.end(); udomain != end; ++udomain)
473  {
474  next_free = std::find_if(next_free, assigned.end(), std::bind(std::less<int>(), std::placeholders::_1, 1));
475  assert(next_free != assigned.end());
476  domainMapping[*udomain] = next_free-assigned.begin();
477  *next_free = 1;
478  }
479  }
480 
481  struct SortFirst
482  {
483  template<class T>
484  bool operator()(const T& t1, const T& t2) const
485  {
486  return t1<t2;
487  }
488  };
489 
490 
501  template<class GI>
502  void mergeVec(std::vector<std::pair<GI, int> >& ownerVec, std::set<GI>& overlapSet) {
503 
504 #ifdef DEBUG_REPART
505  // Safety check for duplicates.
506  if(ownerVec.size()>0)
507  {
508  auto old=ownerVec.begin();
509  for(auto i=old+1, end=ownerVec.end(); i != end; old=i++)
510  {
511  if(i->first==old->first)
512  {
513  std::cerr<<"Value at indes"<<old-ownerVec.begin()<<" is the same as at index "
514  <<i-ownerVec.begin()<<" ["<<old->first<<","<<old->second<<"]==["
515  <<i->first<<","<<i->second<<"]"<<std::endl;
516  throw "Huch!";
517  }
518  }
519  }
520 
521 #endif
522 
523  auto v=ownerVec.begin(), vend=ownerVec.end();
524  for(auto s=overlapSet.begin(), send=overlapSet.end(); s!=send;)
525  {
526  while(v!=vend && v->first<*s) ++v;
527  if(v!=vend && v->first==*s) {
528  // Move to the next element before erasing
529  // thus s stays valid!
530  auto tmp=s;
531  ++s;
532  overlapSet.erase(tmp);
533  }else
534  ++s;
535  }
536  }
537 
538 
552  template<class OwnerSet, class Graph, class IS, class GI>
553  void getNeighbor(const Graph& g, std::vector<int>& part,
554  typename Graph::VertexDescriptor vtx, const IS& indexSet,
555  int toPe, std::set<GI>& neighbor, std::set<int>& neighborProcs) {
556  for(auto edge=g.beginEdges(vtx), end=g.endEdges(vtx); edge!=end; ++edge)
557  {
558  const typename IS::IndexPair* pindex = indexSet.pair(edge.target());
559  assert(pindex);
560  if(part[pindex->local()]!=toPe || !OwnerSet::contains(pindex->local().attribute()))
561  {
562  // is sent to another process and therefore becomes overlap
563  neighbor.insert(pindex->global());
564  neighborProcs.insert(part[pindex->local()]);
565  }
566  }
567  }
568 
569  template<class T, class I>
570  void my_push_back(std::vector<T>& ownerVec, const I& index, [[maybe_unused]] int proc)
571  {
572  ownerVec.push_back(index);
573  }
574 
575  template<class T, class I>
576  void my_push_back(std::vector<std::pair<T,int> >& ownerVec, const I& index, int proc)
577  {
578  ownerVec.push_back(std::make_pair(index,proc));
579  }
580  template<class T>
581  void reserve(std::vector<T>&, RedistributeInterface&, int)
582  {}
583  template<class T>
584  void reserve(std::vector<std::pair<T,int> >& ownerVec, RedistributeInterface& redist, int proc)
585  {
586  redist.reserveSpaceForReceiveInterface(proc, ownerVec.size());
587  }
588 
589 
607  template<class OwnerSet, class G, class IS, class T, class GI>
608  void getOwnerOverlapVec(const G& graph, std::vector<int>& part, IS& indexSet,
609  [[maybe_unused]] int myPe, int toPe, std::vector<T>& ownerVec, std::set<GI>& overlapSet,
610  RedistributeInterface& redist, std::set<int>& neighborProcs) {
611  for(auto index = indexSet.begin(); index != indexSet.end(); ++index) {
612  // Only Process owner vertices, the others are not in the parmetis graph.
613  if(OwnerSet::contains(index->local().attribute()))
614  {
615  if(part[index->local()]==toPe)
616  {
617  getNeighbor<OwnerSet>(graph, part, index->local(), indexSet,
618  toPe, overlapSet, neighborProcs);
619  my_push_back(ownerVec, index->global(), toPe);
620  }
621  }
622  }
623  reserve(ownerVec, redist, toPe);
624 
625  }
626 
627 
634  template<class F, class IS>
635  inline bool isOwner(IS& indexSet, int index) {
636 
637  const typename IS::IndexPair* pindex=indexSet.pair(index);
638 
639  assert(pindex);
640  return F::contains(pindex->local().attribute());
641  }
642 
643 
644  class BaseEdgeFunctor
645  {
646  public:
647  BaseEdgeFunctor(Metis::idx_t* adj,const ParmetisDuneIndexMap& data)
648  : i_(), adj_(adj), data_(data)
649  {}
650 
651  template<class T>
652  void operator()(const T& edge)
653  {
654  // Get the egde weight
655  // const Weight& weight=edge.weight();
656  adj_[i_] = data_.toParmetis(edge.target());
657  i_++;
658  }
659  std::size_t index()
660  {
661  return i_;
662  }
663 
664  private:
665  std::size_t i_;
666  Metis::idx_t* adj_;
667  const ParmetisDuneIndexMap& data_;
668  };
669 
670  template<typename G>
671  struct EdgeFunctor
672  : public BaseEdgeFunctor
673  {
674  EdgeFunctor(Metis::idx_t* adj, const ParmetisDuneIndexMap& data, std::size_t)
675  : BaseEdgeFunctor(adj, data)
676  {}
677 
678  Metis::idx_t* getWeights()
679  {
680  return NULL;
681  }
682  void free(){}
683  };
684 
685  template<class G, class V, class E, class VM, class EM>
686  class EdgeFunctor<Dune::Amg::PropertiesGraph<G,V,E,VM,EM> >
687  : public BaseEdgeFunctor
688  {
689  public:
690  EdgeFunctor(Metis::idx_t* adj, const ParmetisDuneIndexMap& data, std::size_t s)
691  : BaseEdgeFunctor(adj, data)
692  {
693  weight_=new Metis::idx_t[s];
694  }
695 
696  template<class T>
697  void operator()(const T& edge)
698  {
699  weight_[index()]=edge.properties().depends() ? 3 : 1;
700  BaseEdgeFunctor::operator()(edge);
701  }
702  Metis::idx_t* getWeights()
703  {
704  return weight_;
705  }
706  void free(){
707  if(weight_!=0) {
708  delete weight_;
709  weight_=0;
710  }
711  }
712  private:
713  Metis::idx_t* weight_;
714  };
715 
716 
717 
731  template<class F, class G, class IS, class EW>
732  void getAdjArrays(G& graph, IS& indexSet, Metis::idx_t *xadj,
733  EW& ew)
734  {
735  int j=0;
736  auto vend = graph.end();
737 
738  for(auto vertex = graph.begin(); vertex != vend; ++vertex) {
739  if (isOwner<F>(indexSet,*vertex)) {
740  // The type of const edge iterator.
741  auto eend = vertex.end();
742  xadj[j] = ew.index();
743  j++;
744  for(auto edge = vertex.begin(); edge != eend; ++edge) {
745  ew(edge);
746  }
747  }
748  }
749  xadj[j] = ew.index();
750  }
751  } // end anonymous namespace
752 
753  template<class G, class T1, class T2>
754  bool buildCommunication(const G& graph, std::vector<int>& realparts,
756  std::shared_ptr<Dune::OwnerOverlapCopyCommunication<T1,T2>>& outcomm,
757  RedistributeInterface& redistInf,
758  bool verbose=false);
759 #if HAVE_PARMETIS
760 #ifndef METIS_VER_MAJOR
761  extern "C"
762  {
763  // backwards compatibility to parmetis < 4.0.0
764  void METIS_PartGraphKway(int *nvtxs, Metis::idx_t *xadj, Metis::idx_t *adjncy, Metis::idx_t *vwgt,
765  Metis::idx_t *adjwgt, int *wgtflag, int *numflag, int *nparts,
766  int *options, int *edgecut, Metis::idx_t *part);
767 
768  void METIS_PartGraphRecursive(int *nvtxs, Metis::idx_t *xadj, Metis::idx_t *adjncy, Metis::idx_t *vwgt,
769  Metis::idx_t *adjwgt, int *wgtflag, int *numflag, int *nparts,
770  int *options, int *edgecut, Metis::idx_t *part);
771  }
772 #endif
773 #endif // HAVE_PARMETIS
774 
775  template<class S, class T>
776  inline void print_carray(S& os, T* array, std::size_t l)
777  {
778  for(T *cur=array, *end=array+l; cur!=end; ++cur)
779  os<<*cur<<" ";
780  }
781 
782  template<class S, class T>
783  inline bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T* xadj,
784  T* adjncy, bool checkSymmetry)
785  {
786  bool correct=true;
787 
788  using std::signbit;
789  for(Metis::idx_t vtx=0; vtx<(Metis::idx_t)noVtx; ++vtx) {
790  if(static_cast<S>(xadj[vtx])>noEdges || signbit(xadj[vtx])) {
791  std::cerr <<"Check graph: xadj["<<vtx<<"]="<<xadj[vtx]<<" (>"
792  <<noEdges<<") out of range!"<<std::endl;
793  correct=false;
794  }
795  if(static_cast<S>(xadj[vtx+1])>noEdges || signbit(xadj[vtx+1])) {
796  std::cerr <<"Check graph: xadj["<<vtx+1<<"]="<<xadj[vtx+1]<<" (>"
797  <<noEdges<<") out of range!"<<std::endl;
798  correct=false;
799  }
800  // Check numbers in adjncy
801  for(Metis::idx_t i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
802  if(signbit(adjncy[i]) || ((std::size_t)adjncy[i])>gnoVtx) {
803  std::cerr<<" Edge "<<adjncy[i]<<" out of range ["<<0<<","<<noVtx<<")"
804  <<std::endl;
805  correct=false;
806  }
807  }
808  if(checkSymmetry) {
809  for(Metis::idx_t i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
810  Metis::idx_t target=adjncy[i];
811  // search for symmetric edge
812  int found=0;
813  for(Metis::idx_t j=xadj[target]; j< xadj[target+1]; ++j)
814  if(adjncy[j]==vtx)
815  found++;
816  if(found!=1) {
817  std::cerr<<"Edge ("<<target<<","<<vtx<<") "<<i<<" time"<<std::endl;
818  correct=false;
819  }
820  }
821  }
822  }
823  return correct;
824  }
825 
826  template<class M, class T1, class T2>
828  Metis::idx_t nparts,
829  std::shared_ptr<Dune::OwnerOverlapCopyCommunication<T1,T2>>& outcomm,
830  RedistributeInterface& redistInf,
831  bool verbose=false)
832  {
833  if(verbose && oocomm.communicator().rank()==0)
834  std::cout<<"Repartitioning from "<<oocomm.communicator().size()
835  <<" to "<<nparts<<" parts"<<std::endl;
836  Timer time;
837  int rank = oocomm.communicator().rank();
838 #if !HAVE_PARMETIS
839  int* part = new int[1];
840  part[0]=0;
841 #else
842  Metis::idx_t* part = new Metis::idx_t[1]; // where all our data moves to
843 
844  if(nparts>1) {
845 
846  part[0]=rank;
847 
848  { // sublock for automatic memory deletion
849 
850  // Build the graph of the communication scheme and create an appropriate indexset.
851  // calculate the neighbour vertices
852  int noNeighbours = oocomm.remoteIndices().neighbours();
853 
854  for(auto n= oocomm.remoteIndices().begin(); n != oocomm.remoteIndices().end();
855  ++n)
856  if(n->first==rank) {
857  //do not include ourselves.
858  --noNeighbours;
859  break;
860  }
861 
862  // A parmetis graph representing the communication graph.
863  // The diagonal entries are the number of nodes on the process.
864  // The offdiagonal entries are the number of edges leading to other processes.
865 
866  Metis::idx_t *xadj=new Metis::idx_t[2];
867  Metis::idx_t *vtxdist=new Metis::idx_t[oocomm.communicator().size()+1];
868  Metis::idx_t *adjncy=new Metis::idx_t[noNeighbours];
869 #ifdef USE_WEIGHTS
870  Metis::idx_t *vwgt = 0;
871  Metis::idx_t *adjwgt = 0;
872 #endif
873 
874  // each process has exactly one vertex!
875  for(int i=0; i<oocomm.communicator().size(); ++i)
876  vtxdist[i]=i;
877  vtxdist[oocomm.communicator().size()]=oocomm.communicator().size();
878 
879  xadj[0]=0;
880  xadj[1]=noNeighbours;
881 
882  // count edges to other processor
883  // a vector mapping the index to the owner
884  // std::vector<int> owner(mat.N(), oocomm.communicator().rank());
885  // for(NeighbourIterator n= oocomm.remoteIndices().begin(); n != oocomm.remoteIndices().end();
886  // ++n)
887  // {
888  // if(n->first!=oocomm.communicator().rank()){
889  // typedef typename RemoteIndices::RemoteIndexList RIList;
890  // const RIList& rlist = *(n->second.first);
891  // typedef typename RIList::const_iterator LIter;
892  // for(LIter entry=rlist.begin(); entry!=rlist.end(); ++entry){
893  // if(entry->attribute()==OwnerOverlapCopyAttributeSet::owner)
894  // owner[entry->localIndexPair().local()] = n->first;
895  // }
896  // }
897  // }
898 
899  // std::map<int,Metis::idx_t> edgecount; // edges to other processors
900  // typedef typename M::ConstRowIterator RIter;
901  // typedef typename M::ConstColIterator CIter;
902 
903  // // calculate edge count
904  // for(RIter row=mat.begin(), endr=mat.end(); row != endr; ++row)
905  // if(owner[row.index()]==OwnerOverlapCopyAttributeSet::owner)
906  // for(CIter entry= row->begin(), ende = row->end(); entry != ende; ++entry)
907  // ++edgecount[owner[entry.index()]];
908 
909  // setup edge and weight pattern
910 
911  Metis::idx_t* adjp=adjncy;
912 
913 #ifdef USE_WEIGHTS
914  vwgt = new Metis::idx_t[1];
915  vwgt[0]= mat.N(); // weight is numer of rows TODO: Should actually be the nonzeros.
916 
917  adjwgt = new Metis::idx_t[noNeighbours];
918  Metis::idx_t* adjwp=adjwgt;
919 #endif
920 
921  for(auto n= oocomm.remoteIndices().begin(); n != oocomm.remoteIndices().end();
922  ++n)
923  if(n->first != rank) {
924  *adjp=n->first;
925  ++adjp;
926 #ifdef USE_WEIGHTS
927  *adjwp=1; //edgecount[n->first];
928  ++adjwp;
929 #endif
930  }
931  assert(isValidGraph(vtxdist[rank+1]-vtxdist[rank],
932  vtxdist[oocomm.communicator().size()],
933  noNeighbours, xadj, adjncy, false));
934 
935  [[maybe_unused]] Metis::idx_t wgtflag=0;
936  Metis::idx_t numflag=0;
937  Metis::idx_t edgecut;
938 #ifdef USE_WEIGHTS
939  wgtflag=3;
940 #endif
941  Metis::real_t *tpwgts = new Metis::real_t[nparts];
942  for(int i=0; i<nparts; ++i)
943  tpwgts[i]=1.0/nparts;
944  MPI_Comm comm=oocomm.communicator();
945 
946  Dune::dinfo<<rank<<" vtxdist: ";
947  print_carray(Dune::dinfo, vtxdist, oocomm.communicator().size()+1);
948  Dune::dinfo<<std::endl<<rank<<" xadj: ";
949  print_carray(Dune::dinfo, xadj, 2);
950  Dune::dinfo<<std::endl<<rank<<" adjncy: ";
951  print_carray(Dune::dinfo, adjncy, noNeighbours);
952 
953 #ifdef USE_WEIGHTS
954  Dune::dinfo<<std::endl<<rank<<" vwgt: ";
955  print_carray(Dune::dinfo, vwgt, 1);
956  Dune::dinfo<<std::endl<<rank<<" adwgt: ";
957  print_carray(Dune::dinfo, adjwgt, noNeighbours);
958 #endif
959  Dune::dinfo<<std::endl;
960  oocomm.communicator().barrier();
961  if(verbose && oocomm.communicator().rank()==0)
962  std::cout<<"Creating comm graph took "<<time.elapsed()<<std::endl;
963  time.reset();
964 
965 #ifdef PARALLEL_PARTITION
966  Metis::real_t ubvec = 1.15;
967  int ncon=1;
968  int options[5] ={ 0,1,15,0,0};
969 
970  //=======================================================
971  // ParMETIS_V3_PartKway
972  //=======================================================
973  ParMETIS_V3_PartKway(vtxdist, xadj, adjncy,
974  vwgt, adjwgt, &wgtflag,
975  &numflag, &ncon, &nparts, tpwgts, &ubvec, options, &edgecut, part,
976  &comm);
977  if(verbose && oocomm.communicator().rank()==0)
978  std::cout<<"ParMETIS took "<<time.elapsed()<<std::endl;
979  time.reset();
980 #else
981  Timer time1;
982  std::size_t gnoedges=0;
983  int* noedges = 0;
984  noedges = new int[oocomm.communicator().size()];
985  Dune::dverb<<"noNeighbours: "<<noNeighbours<<std::endl;
986  // gather number of edges for each vertex.
987  MPI_Allgather(&noNeighbours,1,MPI_INT,noedges,1, MPI_INT,oocomm.communicator());
988 
989  if(verbose && oocomm.communicator().rank()==0)
990  std::cout<<"Gathering noedges took "<<time1.elapsed()<<std::endl;
991  time1.reset();
992 
993  Metis::idx_t noVertices = vtxdist[oocomm.communicator().size()];
994  Metis::idx_t *gxadj = 0;
995  Metis::idx_t *gvwgt = 0;
996  Metis::idx_t *gadjncy = 0;
997  Metis::idx_t *gadjwgt = 0;
998  Metis::idx_t *gpart = 0;
999  int* displ = 0;
1000  int* noxs = 0;
1001  int* xdispl = 0; // displacement for xadj
1002  int* novs = 0;
1003  int* vdispl=0; // real vertex displacement
1004 #ifdef USE_WEIGHTS
1005  std::size_t localNoVtx=vtxdist[rank+1]-vtxdist[rank];
1006 #endif
1007  std::size_t gxadjlen = vtxdist[oocomm.communicator().size()]-vtxdist[0]+oocomm.communicator().size();
1008 
1009  {
1010  Dune::dinfo<<"noedges: ";
1011  print_carray(Dune::dinfo, noedges, oocomm.communicator().size());
1012  Dune::dinfo<<std::endl;
1013  displ = new int[oocomm.communicator().size()];
1014  xdispl = new int[oocomm.communicator().size()];
1015  noxs = new int[oocomm.communicator().size()];
1016  vdispl = new int[oocomm.communicator().size()];
1017  novs = new int[oocomm.communicator().size()];
1018 
1019  for(int i=0; i < oocomm.communicator().size(); ++i) {
1020  noxs[i]=vtxdist[i+1]-vtxdist[i]+1;
1021  novs[i]=vtxdist[i+1]-vtxdist[i];
1022  }
1023 
1024  Metis::idx_t *so= vtxdist;
1025  int offset = 0;
1026  for(int *xcurr = xdispl, *vcurr = vdispl, *end=vdispl+oocomm.communicator().size();
1027  vcurr!=end; ++vcurr, ++xcurr, ++so, ++offset) {
1028  *vcurr = *so;
1029  *xcurr = offset + *so;
1030  }
1031 
1032  int *pdispl =displ;
1033  int cdispl = 0;
1034  *pdispl = 0;
1035  for(int *curr=noedges, *end=noedges+oocomm.communicator().size()-1;
1036  curr!=end; ++curr) {
1037  ++pdispl; // next displacement
1038  cdispl += *curr; // next value
1039  *pdispl = cdispl;
1040  }
1041  Dune::dinfo<<"displ: ";
1042  print_carray(Dune::dinfo, displ, oocomm.communicator().size());
1043  Dune::dinfo<<std::endl;
1044 
1045  // calculate global number of edges
1046  // It is bigger than the actual one as we habe size-1 additional end entries
1047  for(int *curr=noedges, *end=noedges+oocomm.communicator().size();
1048  curr!=end; ++curr)
1049  gnoedges += *curr;
1050 
1051  // alocate gobal graph
1052  Dune::dinfo<<"gxadjlen: "<<gxadjlen<<" noVertices: "<<noVertices
1053  <<" gnoedges: "<<gnoedges<<std::endl;
1054  gxadj = new Metis::idx_t[gxadjlen];
1055  gpart = new Metis::idx_t[noVertices];
1056 #ifdef USE_WEIGHTS
1057  gvwgt = new Metis::idx_t[noVertices];
1058  gadjwgt = new Metis::idx_t[gnoedges];
1059 #endif
1060  gadjncy = new Metis::idx_t[gnoedges];
1061  }
1062 
1063  if(verbose && oocomm.communicator().rank()==0)
1064  std::cout<<"Preparing global graph took "<<time1.elapsed()<<std::endl;
1065  time1.reset();
1066  // Communicate data
1067 
1068  MPI_Allgatherv(xadj,2,MPITraits<Metis::idx_t>::getType(),
1069  gxadj,noxs,xdispl,MPITraits<Metis::idx_t>::getType(),
1070  comm);
1071  MPI_Allgatherv(adjncy,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1072  gadjncy,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1073  comm);
1074 #ifdef USE_WEIGHTS
1075  MPI_Allgatherv(adjwgt,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1076  gadjwgt,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1077  comm);
1078  MPI_Allgatherv(vwgt,localNoVtx,MPITraits<Metis::idx_t>::getType(),
1079  gvwgt,novs,vdispl,MPITraits<Metis::idx_t>::getType(),
1080  comm);
1081 #endif
1082  if(verbose && oocomm.communicator().rank()==0)
1083  std::cout<<"Gathering global graph data took "<<time1.elapsed()<<std::endl;
1084  time1.reset();
1085 
1086  {
1087  // create the real gxadj array
1088  // i.e. shift entries and add displacements.
1089 
1090  print_carray(Dune::dinfo, gxadj, gxadjlen);
1091 
1092  int offset = 0;
1093  Metis::idx_t increment = vtxdist[1];
1094  Metis::idx_t *start=gxadj+1;
1095  for(int i=1; i<oocomm.communicator().size(); ++i) {
1096  offset+=1;
1097  int lprev = vtxdist[i]-vtxdist[i-1];
1098  int l = vtxdist[i+1]-vtxdist[i];
1099  start+=lprev;
1100  assert((start+l+offset)-gxadj<=static_cast<Metis::idx_t>(gxadjlen));
1101  increment = *(start-1);
1102  std::transform(start+offset, start+l+offset, start, std::bind(std::plus<Metis::idx_t>(), std::placeholders::_1, increment));
1103  }
1104  Dune::dinfo<<std::endl<<"shifted xadj:";
1105  print_carray(Dune::dinfo, gxadj, noVertices+1);
1106  Dune::dinfo<<std::endl<<" gadjncy: ";
1107  print_carray(Dune::dinfo, gadjncy, gnoedges);
1108 #ifdef USE_WEIGHTS
1109  Dune::dinfo<<std::endl<<" gvwgt: ";
1110  print_carray(Dune::dinfo, gvwgt, noVertices);
1111  Dune::dinfo<<std::endl<<"adjwgt: ";
1112  print_carray(Dune::dinfo, gadjwgt, gnoedges);
1113  Dune::dinfo<<std::endl;
1114 #endif
1115  // everything should be fine now!!!
1116  if(verbose && oocomm.communicator().rank()==0)
1117  std::cout<<"Postprocesing global graph data took "<<time1.elapsed()<<std::endl;
1118  time1.reset();
1119 #ifndef NDEBUG
1120  assert(isValidGraph(noVertices, noVertices, gnoedges,
1121  gxadj, gadjncy, true));
1122 #endif
1123 
1124  if(verbose && oocomm.communicator().rank()==0)
1125  std::cout<<"Creating grah one 1 process took "<<time.elapsed()<<std::endl;
1126  time.reset();
1127 #if METIS_VER_MAJOR >= 5
1128  Metis::idx_t ncon = 1;
1129  Metis::idx_t moptions[METIS_NOPTIONS];
1130  METIS_SetDefaultOptions(moptions);
1131  moptions[METIS_OPTION_NUMBERING] = numflag;
1132  METIS_PartGraphRecursive(&noVertices, &ncon, gxadj, gadjncy, gvwgt, NULL, gadjwgt,
1133  &nparts, NULL, NULL, moptions, &edgecut, gpart);
1134 #else
1135  int options[5] = {0, 1, 1, 3, 3};
1136  // Call metis
1137  METIS_PartGraphRecursive(&noVertices, gxadj, gadjncy, gvwgt, gadjwgt, &wgtflag,
1138  &numflag, &nparts, options, &edgecut, gpart);
1139 #endif
1140 
1141  if(verbose && oocomm.communicator().rank()==0)
1142  std::cout<<"METIS took "<<time.elapsed()<<std::endl;
1143  time.reset();
1144 
1145  Dune::dinfo<<std::endl<<"part:";
1146  print_carray(Dune::dinfo, gpart, noVertices);
1147 
1148  delete[] gxadj;
1149  delete[] gadjncy;
1150 #ifdef USE_WEIGHTS
1151  delete[] gvwgt;
1152  delete[] gadjwgt;
1153 #endif
1154  }
1155  // Scatter result
1156  MPI_Scatter(gpart, 1, MPITraits<Metis::idx_t>::getType(), part, 1,
1157  MPITraits<Metis::idx_t>::getType(), 0, comm);
1158 
1159  {
1160  // release remaining memory
1161  delete[] gpart;
1162  delete[] noedges;
1163  delete[] displ;
1164  }
1165 
1166 
1167 #endif
1168  delete[] xadj;
1169  delete[] vtxdist;
1170  delete[] adjncy;
1171 #ifdef USE_WEIGHTS
1172  delete[] vwgt;
1173  delete[] adjwgt;
1174 #endif
1175  delete[] tpwgts;
1176  }
1177  }else{
1178  part[0]=0;
1179  }
1180 #endif
1181  Dune::dinfo<<" repart "<<rank <<" -> "<< part[0]<<std::endl;
1182 
1183  std::vector<int> realpart(mat.N(), part[0]);
1184  delete[] part;
1185 
1186  oocomm.copyOwnerToAll(realpart, realpart);
1187 
1188  if(verbose && oocomm.communicator().rank()==0)
1189  std::cout<<"Scattering repartitioning took "<<time.elapsed()<<std::endl;
1190  time.reset();
1191 
1192 
1193  oocomm.buildGlobalLookup(mat.N());
1194  Dune::Amg::MatrixGraph<M> graph(const_cast<M&>(mat));
1195  fillIndexSetHoles(graph, oocomm);
1196  if(verbose && oocomm.communicator().rank()==0)
1197  std::cout<<"Filling index set took "<<time.elapsed()<<std::endl;
1198  time.reset();
1199 
1200  if(verbose) {
1201  int noNeighbours=oocomm.remoteIndices().neighbours();
1202  noNeighbours = oocomm.communicator().sum(noNeighbours)
1203  / oocomm.communicator().size();
1204  if(oocomm.communicator().rank()==0)
1205  std::cout<<"Average no neighbours was "<<noNeighbours<<std::endl;
1206  }
1207  bool ret = buildCommunication(graph, realpart, oocomm, outcomm, redistInf,
1208  verbose);
1209  if(verbose && oocomm.communicator().rank()==0)
1210  std::cout<<"Building index sets took "<<time.elapsed()<<std::endl;
1211  time.reset();
1212 
1213 
1214  return ret;
1215 
1216  }
1217 
1232  template<class G, class T1, class T2>
1234  std::shared_ptr<Dune::OwnerOverlapCopyCommunication<T1,T2>>& outcomm,
1235  RedistributeInterface& redistInf,
1236  bool verbose=false)
1237  {
1238  Timer time;
1239 
1240  MPI_Comm comm=oocomm.communicator();
1241  oocomm.buildGlobalLookup(graph.noVertices());
1242  fillIndexSetHoles(graph, oocomm);
1243 
1244  if(verbose && oocomm.communicator().rank()==0)
1245  std::cout<<"Filling holes took "<<time.elapsed()<<std::endl;
1246  time.reset();
1247 
1248  // simple precondition checks
1249 
1250 #ifdef PERF_REPART
1251  // Profiling variables
1252  double t1=0.0, t2=0.0, t3=0.0, t4=0.0, tSum=0.0;
1253 #endif
1254 
1255 
1256  // MPI variables
1257  int mype = oocomm.communicator().rank();
1258 
1259  assert(nparts<=static_cast<Metis::idx_t>(oocomm.communicator().size()));
1260 
1261  int myDomain = -1;
1262 
1263  //
1264  // 1) Prepare the required parameters for using ParMETIS
1265  // Especially the arrays that represent the graph must be
1266  // generated by the DUNE Graph and IndexSet input variables.
1267  // These are the arrays:
1268  // - vtxdist
1269  // - xadj
1270  // - adjncy
1271  //
1272  //
1273 #ifdef PERF_REPART
1274  // reset timer for step 1)
1275  t1=MPI_Wtime();
1276 #endif
1277 
1278 
1279  typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2> OOComm;
1280  typedef typename OOComm::OwnerSet OwnerSet;
1281 
1282  // Create the vtxdist array and parmetisVtxMapping.
1283  // Global communications are necessary
1284  // The parmetis global identifiers for the owner vertices.
1285  ParmetisDuneIndexMap indexMap(graph,oocomm);
1286  Metis::idx_t *part = new Metis::idx_t[indexMap.numOfOwnVtx()];
1287  for(std::size_t i=0; i < indexMap.numOfOwnVtx(); ++i)
1288  part[i]=mype;
1289 
1290 #if !HAVE_PARMETIS
1291  if(oocomm.communicator().rank()==0 && nparts>1)
1292  std::cerr<<"ParMETIS not activated. Will repartition to 1 domain instead of requested "
1293  <<nparts<<" domains."<<std::endl;
1294  nparts=1; // No parmetis available, fallback to agglomerating to 1 process
1295 
1296 #else
1297 
1298  if(nparts>1) {
1299  // Create the xadj and adjncy arrays
1300  Metis::idx_t *xadj = new Metis::idx_t[indexMap.numOfOwnVtx()+1];
1301  Metis::idx_t *adjncy = new Metis::idx_t[graph.noEdges()];
1302  EdgeFunctor<G> ef(adjncy, indexMap, graph.noEdges());
1303  getAdjArrays<OwnerSet>(graph, oocomm.globalLookup(), xadj, ef);
1304 
1305  //
1306  // 2) Call ParMETIS
1307  //
1308  //
1309  Metis::idx_t numflag=0, wgtflag=0, options[3], edgecut=0, ncon=1;
1310  //float *tpwgts = NULL;
1311  Metis::real_t *tpwgts = new Metis::real_t[nparts];
1312  for(int i=0; i<nparts; ++i)
1313  tpwgts[i]=1.0/nparts;
1314  Metis::real_t ubvec[1];
1315  options[0] = 0; // 0=default, 1=options are defined in [1]+[2]
1316 #ifdef DEBUG_REPART
1317  options[1] = 3; // show info: 0=no message
1318 #else
1319  options[1] = 0; // show info: 0=no message
1320 #endif
1321  options[2] = 1; // random number seed, default is 15
1322  wgtflag = (ef.getWeights()!=NULL) ? 1 : 0;
1323  numflag = 0;
1324  edgecut = 0;
1325  ncon=1;
1326  ubvec[0]=1.05; // recommended by ParMETIS
1327 
1328 #ifdef DEBUG_REPART
1329  if (mype == 0) {
1330  std::cout<<std::endl;
1331  std::cout<<"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1332  <<options[1]<<" "<<options[2]<<"}, Ncon: "
1333  <<ncon<<", Nparts: "<<nparts<<std::endl;
1334  }
1335 #endif
1336 #ifdef PERF_REPART
1337  // stop the time for step 1)
1338  t1=MPI_Wtime()-t1;
1339  // reset timer for step 2)
1340  t2=MPI_Wtime();
1341 #endif
1342 
1343  if(verbose) {
1344  oocomm.communicator().barrier();
1345  if(oocomm.communicator().rank()==0)
1346  std::cout<<"Preparing for parmetis took "<<time.elapsed()<<std::endl;
1347  }
1348  time.reset();
1349 
1350  //=======================================================
1351  // ParMETIS_V3_PartKway
1352  //=======================================================
1353  ParMETIS_V3_PartKway(indexMap.vtxDist(), xadj, adjncy,
1354  NULL, ef.getWeights(), &wgtflag,
1355  &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &const_cast<MPI_Comm&>(comm));
1356 
1357 
1358  delete[] xadj;
1359  delete[] adjncy;
1360  delete[] tpwgts;
1361 
1362  ef.free();
1363 
1364 #ifdef DEBUG_REPART
1365  if (mype == 0) {
1366  std::cout<<std::endl;
1367  std::cout<<"ParMETIS_V3_PartKway reported a cut of "<<edgecut<<std::endl;
1368  std::cout<<std::endl;
1369  }
1370  std::cout<<mype<<": PARMETIS-Result: ";
1371  for(int i=0; i < indexMap.vtxDist()[mype+1]-indexMap.vtxDist()[mype]; ++i) {
1372  std::cout<<part[i]<<" ";
1373  }
1374  std::cout<<std::endl;
1375  std::cout<<"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1376  <<options[1]<<" "<<options[2]<<"}, Ncon: "
1377  <<ncon<<", Nparts: "<<nparts<<std::endl;
1378 #endif
1379 #ifdef PERF_REPART
1380  // stop the time for step 2)
1381  t2=MPI_Wtime()-t2;
1382  // reset timer for step 3)
1383  t3=MPI_Wtime();
1384 #endif
1385 
1386 
1387  if(verbose) {
1388  oocomm.communicator().barrier();
1389  if(oocomm.communicator().rank()==0)
1390  std::cout<<"Parmetis took "<<time.elapsed()<<std::endl;
1391  }
1392  time.reset();
1393  }else
1394 #endif
1395  {
1396  // Everything goes to process 0!
1397  for(std::size_t i=0; i<indexMap.numOfOwnVtx(); ++i)
1398  part[i]=0;
1399  }
1400 
1401 
1402  //
1403  // 3) Find a optimal domain based on the ParMETIS repartitioning
1404  // result
1405  //
1406 
1407  std::vector<int> domainMapping(nparts);
1408  if(nparts>1)
1409  getDomain(comm, part, indexMap.numOfOwnVtx(), nparts, &myDomain, domainMapping);
1410  else
1411  domainMapping[0]=0;
1412 
1413 #ifdef DEBUG_REPART
1414  std::cout<<mype<<": myDomain: "<<myDomain<<std::endl;
1415  std::cout<<mype<<": DomainMapping: ";
1416  for(auto j : range(nparts)) {
1417  std::cout<<" do: "<<j<<" pe: "<<domainMapping[j]<<" ";
1418  }
1419  std::cout<<std::endl;
1420 #endif
1421 
1422  // Make a domain mapping for the indexset and translate
1423  //domain number to real process number
1424  // domainMapping is the one of parmetis, that is without
1425  // the overlap/copy vertices
1426  std::vector<int> setPartition(oocomm.indexSet().size(), -1);
1427 
1428  std::size_t i=0; // parmetis index
1429  for(auto index = oocomm.indexSet().begin(); index != oocomm.indexSet().end(); ++index)
1430  if(OwnerSet::contains(index->local().attribute())) {
1431  setPartition[index->local()]=domainMapping[part[i++]];
1432  }
1433 
1434  delete[] part;
1435  oocomm.copyOwnerToAll(setPartition, setPartition);
1436  // communication only needed for ALU
1437  // (ghosts with same global id as owners on the same process)
1438  if (SolverCategory::category(oocomm) ==
1439  static_cast<int>(SolverCategory::nonoverlapping))
1440  oocomm.copyCopyToAll(setPartition, setPartition);
1441  bool ret = buildCommunication(graph, setPartition, oocomm, outcomm, redistInf,
1442  verbose);
1443  if(verbose) {
1444  oocomm.communicator().barrier();
1445  if(oocomm.communicator().rank()==0)
1446  std::cout<<"Creating indexsets took "<<time.elapsed()<<std::endl;
1447  }
1448  return ret;
1449  }
1450 
1451 
1452 
1453  template<class G, class T1, class T2>
1454  bool buildCommunication(const G& graph,
1455  std::vector<int>& setPartition, Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm,
1456  std::shared_ptr<Dune::OwnerOverlapCopyCommunication<T1,T2>>& outcomm,
1457  RedistributeInterface& redistInf,
1458  bool verbose)
1459  {
1460  typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2> OOComm;
1461  typedef typename OOComm::OwnerSet OwnerSet;
1462 
1463  Timer time;
1464 
1465  // Build the send interface
1466  redistInf.buildSendInterface<OwnerSet>(setPartition, oocomm.indexSet());
1467 
1468 #ifdef PERF_REPART
1469  // stop the time for step 3)
1470  t3=MPI_Wtime()-t3;
1471  // reset timer for step 4)
1472  t4=MPI_Wtime();
1473 #endif
1474 
1475 
1476  //
1477  // 4) Create the output IndexSet and RemoteIndices
1478  // 4.1) Determine the "send to" and "receive from" relation
1479  // according to the new partition using a MPI ring
1480  // communication.
1481  //
1482  // 4.2) Depends on the "send to" and "receive from" vector,
1483  // the processes will exchange the vertices each other
1484  //
1485  // 4.3) Create the IndexSet, RemoteIndices and the new MPI
1486  // communicator
1487  //
1488 
1489  //
1490  // 4.1) Let's start...
1491  //
1492  int npes = oocomm.communicator().size();
1493  int *sendTo = 0;
1494  int noSendTo = 0;
1495  std::set<int> recvFrom;
1496 
1497  // the max number of vertices is stored in the sendTo buffer,
1498  // not the number of vertices to send! Because the max number of Vtx
1499  // is used as the fixed buffer size by the MPI send/receive calls
1500 
1501  int mype = oocomm.communicator().rank();
1502 
1503  {
1504  std::set<int> tsendTo;
1505  for(auto i=setPartition.begin(), iend = setPartition.end(); i!=iend; ++i)
1506  tsendTo.insert(*i);
1507 
1508  noSendTo = tsendTo.size();
1509  sendTo = new int[noSendTo];
1510  int idx=0;
1511  for(auto i=tsendTo.begin(); i != tsendTo.end(); ++i, ++idx)
1512  sendTo[idx]=*i;
1513  }
1514 
1515  //
1516  int* gnoSend= new int[oocomm.communicator().size()];
1517  int* gsendToDispl = new int[oocomm.communicator().size()+1];
1518 
1519  MPI_Allgather(&noSendTo, 1, MPI_INT, gnoSend, 1,
1520  MPI_INT, oocomm.communicator());
1521 
1522  // calculate total receive message size
1523  int totalNoRecv = 0;
1524  for(int i=0; i<npes; ++i)
1525  totalNoRecv += gnoSend[i];
1526 
1527  int *gsendTo = new int[totalNoRecv];
1528 
1529  // calculate displacement for allgatherv
1530  gsendToDispl[0]=0;
1531  for(int i=0; i<npes; ++i)
1532  gsendToDispl[i+1]=gsendToDispl[i]+gnoSend[i];
1533 
1534  // gather the data
1535  MPI_Allgatherv(sendTo, noSendTo, MPI_INT, gsendTo, gnoSend, gsendToDispl,
1536  MPI_INT, oocomm.communicator());
1537 
1538  // Extract from which processes we will receive data
1539  for(int proc=0; proc < npes; ++proc)
1540  for(int i=gsendToDispl[proc]; i < gsendToDispl[proc+1]; ++i)
1541  if(gsendTo[i]==mype)
1542  recvFrom.insert(proc);
1543 
1544  bool existentOnNextLevel = recvFrom.size()>0;
1545 
1546  // Delete memory
1547  delete[] gnoSend;
1548  delete[] gsendToDispl;
1549  delete[] gsendTo;
1550 
1551 
1552 #ifdef DEBUG_REPART
1553  if(recvFrom.size()) {
1554  std::cout<<mype<<": recvFrom: ";
1555  for(auto i=recvFrom.begin(); i!= recvFrom.end(); ++i) {
1556  std::cout<<*i<<" ";
1557  }
1558  }
1559 
1560  std::cout<<std::endl<<std::endl;
1561  std::cout<<mype<<": sendTo: ";
1562  for(int i=0; i<noSendTo; i++) {
1563  std::cout<<sendTo[i]<<" ";
1564  }
1565  std::cout<<std::endl<<std::endl;
1566 #endif
1567 
1568  if(verbose)
1569  if(oocomm.communicator().rank()==0)
1570  std::cout<<" Communicating the receive information took "<<
1571  time.elapsed()<<std::endl;
1572  time.reset();
1573 
1574  //
1575  // 4.2) Start the communication
1576  //
1577 
1578  // Get all the owner and overlap vertices for myself ans save
1579  // it in the vectors myOwnerVec and myOverlapVec.
1580  // The received vertices from the other processes are simple
1581  // added to these vector.
1582  //
1583 
1584 
1585  typedef typename OOComm::ParallelIndexSet::GlobalIndex GI;
1586  typedef std::vector<GI> GlobalVector;
1587  std::vector<std::pair<GI,int> > myOwnerVec;
1588  std::set<GI> myOverlapSet;
1589  GlobalVector sendOwnerVec;
1590  std::set<GI> sendOverlapSet;
1591  std::set<int> myNeighbors;
1592 
1593  // getOwnerOverlapVec<OwnerSet>(graph, setPartition, oocomm.globalLookup(),
1594  // mype, mype, myOwnerVec, myOverlapSet, redistInf, myNeighbors);
1595 
1596  char **sendBuffers=new char*[noSendTo];
1597  MPI_Request *requests = new MPI_Request[noSendTo];
1598 
1599  // Create all messages to be sent
1600  for(int i=0; i < noSendTo; ++i) {
1601  // clear the vector for sending
1602  sendOwnerVec.clear();
1603  sendOverlapSet.clear();
1604  // get all owner and overlap vertices for process j and save these
1605  // in the vectors sendOwnerVec and sendOverlapSet
1606  std::set<int> neighbors;
1607  getOwnerOverlapVec<OwnerSet>(graph, setPartition, oocomm.globalLookup(),
1608  mype, sendTo[i], sendOwnerVec, sendOverlapSet, redistInf,
1609  neighbors);
1610  // +2, we need 2 integer more for the length of each part
1611  // (owner/overlap) of the array
1612  int buffersize=0;
1613  int tsize;
1614  MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &buffersize);
1615  MPI_Pack_size(sendOwnerVec.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1616  buffersize +=tsize;
1617  MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1618  buffersize +=tsize;
1619  MPI_Pack_size(sendOverlapSet.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1620  buffersize += tsize;
1621  MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1622  buffersize += tsize;
1623  MPI_Pack_size(neighbors.size(), MPI_INT, oocomm.communicator(), &tsize);
1624  buffersize += tsize;
1625 
1626  sendBuffers[i] = new char[buffersize];
1627 
1628 #ifdef DEBUG_REPART
1629  std::cout<<mype<<" sending "<<sendOwnerVec.size()<<" owner and "<<
1630  sendOverlapSet.size()<<" overlap to "<<sendTo[i]<<" buffersize="<<buffersize<<std::endl;
1631 #endif
1632  createSendBuf(sendOwnerVec, sendOverlapSet, neighbors, sendBuffers[i], buffersize, oocomm.communicator());
1633  MPI_Issend(sendBuffers[i], buffersize, MPI_PACKED, sendTo[i], 99, oocomm.communicator(), requests+i);
1634  }
1635 
1636  if(verbose) {
1637  oocomm.communicator().barrier();
1638  if(oocomm.communicator().rank()==0)
1639  std::cout<<" Creating sends took "<<
1640  time.elapsed()<<std::endl;
1641  }
1642  time.reset();
1643 
1644  // Receive Messages
1645  int noRecv = recvFrom.size();
1646  int oldbuffersize=0;
1647  char* recvBuf = 0;
1648  while(noRecv>0) {
1649  // probe for an incoming message
1650  MPI_Status stat;
1651  MPI_Probe(MPI_ANY_SOURCE, 99, oocomm.communicator(), &stat);
1652  int buffersize;
1653  MPI_Get_count(&stat, MPI_PACKED, &buffersize);
1654 
1655  if(oldbuffersize<buffersize) {
1656  // buffer too small, reallocate
1657  delete[] recvBuf;
1658  recvBuf = new char[buffersize];
1659  oldbuffersize = buffersize;
1660  }
1661  MPI_Recv(recvBuf, buffersize, MPI_PACKED, stat.MPI_SOURCE, 99, oocomm.communicator(), &stat);
1662  saveRecvBuf(recvBuf, buffersize, myOwnerVec, myOverlapSet, myNeighbors, redistInf,
1663  stat.MPI_SOURCE, oocomm.communicator());
1664  --noRecv;
1665  }
1666 
1667  if(recvBuf)
1668  delete[] recvBuf;
1669 
1670  time.reset();
1671  // Wait for sending messages to complete
1672  MPI_Status *statuses = new MPI_Status[noSendTo];
1673  int send = MPI_Waitall(noSendTo, requests, statuses);
1674 
1675  // check for errors
1676  if(send==MPI_ERR_IN_STATUS) {
1677  std::cerr<<mype<<": Error in sending :"<<std::endl;
1678  // Search for the error
1679  for(int i=0; i< noSendTo; i++)
1680  if(statuses[i].MPI_ERROR!=MPI_SUCCESS) {
1681  char message[300];
1682  int messageLength;
1683  MPI_Error_string(statuses[i].MPI_ERROR, message, &messageLength);
1684  std::cerr<<" source="<<statuses[i].MPI_SOURCE<<" message: ";
1685  for(int j = 0; j < messageLength; j++)
1686  std::cout<<message[j];
1687  }
1688  std::cerr<<std::endl;
1689  }
1690 
1691  if(verbose) {
1692  oocomm.communicator().barrier();
1693  if(oocomm.communicator().rank()==0)
1694  std::cout<<" Receiving and saving took "<<
1695  time.elapsed()<<std::endl;
1696  }
1697  time.reset();
1698 
1699  for(int i=0; i < noSendTo; ++i)
1700  delete[] sendBuffers[i];
1701 
1702  delete[] sendBuffers;
1703  delete[] statuses;
1704  delete[] requests;
1705 
1706  redistInf.setCommunicator(oocomm.communicator());
1707 
1708  //
1709  // 4.2) Create the IndexSet etc.
1710  //
1711 
1712  // build the new outputIndexSet
1713 
1714 
1715  int color=0;
1716 
1717  if (!existentOnNextLevel) {
1718  // this process is not used anymore
1719  color= MPI_UNDEFINED;
1720  }
1721  MPI_Comm outputComm;
1722 
1723  MPI_Comm_split(oocomm.communicator(), color, oocomm.communicator().rank(), &outputComm);
1724  outcomm = std::make_shared<OOComm>(outputComm,SolverCategory::category(oocomm),true);
1725 
1726  // translate neighbor ranks.
1727  int newrank=outcomm->communicator().rank();
1728  int *newranks=new int[oocomm.communicator().size()];
1729  std::vector<int> tneighbors;
1730  tneighbors.reserve(myNeighbors.size());
1731 
1732  typename OOComm::ParallelIndexSet& outputIndexSet = outcomm->indexSet();
1733 
1734  MPI_Allgather(&newrank, 1, MPI_INT, newranks, 1,
1735  MPI_INT, oocomm.communicator());
1736 
1737 #ifdef DEBUG_REPART
1738  std::cout<<oocomm.communicator().rank()<<" ";
1739  for(auto i=myNeighbors.begin(), end=myNeighbors.end();
1740  i!=end; ++i) {
1741  assert(newranks[*i]>=0);
1742  std::cout<<*i<<"->"<<newranks[*i]<<" ";
1743  tneighbors.push_back(newranks[*i]);
1744  }
1745  std::cout<<std::endl;
1746 #else
1747  for(auto i=myNeighbors.begin(), end=myNeighbors.end();
1748  i!=end; ++i) {
1749  tneighbors.push_back(newranks[*i]);
1750  }
1751 #endif
1752  delete[] newranks;
1753  myNeighbors.clear();
1754 
1755  if(verbose) {
1756  oocomm.communicator().barrier();
1757  if(oocomm.communicator().rank()==0)
1758  std::cout<<" Calculating new neighbours ("<<tneighbors.size()<<") took "<<
1759  time.elapsed()<<std::endl;
1760  }
1761  time.reset();
1762 
1763 
1764  outputIndexSet.beginResize();
1765  // 1) add the owner vertices
1766  // Sort the owners
1767  std::sort(myOwnerVec.begin(), myOwnerVec.end(), SortFirst());
1768  // The owners are sorted according to there global index
1769  // Therefore the entries of ownerVec are the same as the
1770  // ones in the resulting index set.
1771  int i=0;
1772  using LocalIndexT = typename OOComm::ParallelIndexSet::LocalIndex;
1773  for(auto g=myOwnerVec.begin(), end =myOwnerVec.end(); g!=end; ++g, ++i ) {
1774  outputIndexSet.add(g->first,LocalIndexT(i, OwnerOverlapCopyAttributeSet::owner, true));
1775  redistInf.addReceiveIndex(g->second, i);
1776  }
1777 
1778  if(verbose) {
1779  oocomm.communicator().barrier();
1780  if(oocomm.communicator().rank()==0)
1781  std::cout<<" Adding owner indices took "<<
1782  time.elapsed()<<std::endl;
1783  }
1784  time.reset();
1785 
1786 
1787  // After all the vertices are received, the vectors must
1788  // be "merged" together to create the final vectors.
1789  // Because some vertices that are sent as overlap could now
1790  // already included as owner vertiecs in the new partition
1791  mergeVec(myOwnerVec, myOverlapSet);
1792 
1793  // Trick to free memory
1794  myOwnerVec.clear();
1795  myOwnerVec.swap(myOwnerVec);
1796 
1797  if(verbose) {
1798  oocomm.communicator().barrier();
1799  if(oocomm.communicator().rank()==0)
1800  std::cout<<" Merging indices took "<<
1801  time.elapsed()<<std::endl;
1802  }
1803  time.reset();
1804 
1805 
1806  // 2) add the overlap vertices
1807  for(auto g=myOverlapSet.begin(), end=myOverlapSet.end(); g!=end; ++g, i++) {
1808  outputIndexSet.add(*g,LocalIndexT(i, OwnerOverlapCopyAttributeSet::copy, true));
1809  }
1810  myOverlapSet.clear();
1811  outputIndexSet.endResize();
1812 
1813 #ifdef DUNE_ISTL_WITH_CHECKING
1814  int numOfOwnVtx =0;
1815  auto end = outputIndexSet.end();
1816  for(auto index = outputIndexSet.begin(); index != end; ++index) {
1817  if (OwnerSet::contains(index->local().attribute())) {
1818  numOfOwnVtx++;
1819  }
1820  }
1821  numOfOwnVtx = oocomm.communicator().sum(numOfOwnVtx);
1822  // if(numOfOwnVtx!=indexMap.globalOwnerVertices)
1823  // {
1824  // std::cerr<<numOfOwnVtx<<"!="<<indexMap.globalOwnerVertices<<" owners missing or additional ones!"<<std::endl;
1825  // DUNE_THROW(ISTLError, numOfOwnVtx<<"!="<<indexMap.globalOwnerVertices<<" owners missing or additional ones"
1826  // <<" during repartitioning.");
1827  // }
1828  std::is_sorted(outputIndexSet.begin(), outputIndexSet.end(),
1829  [](const auto& v1, const auto& v2){ return v1.global() < v2.global();});
1830 #endif
1831  if(verbose) {
1832  oocomm.communicator().barrier();
1833  if(oocomm.communicator().rank()==0)
1834  std::cout<<" Adding overlap indices took "<<
1835  time.elapsed()<<std::endl;
1836  }
1837  time.reset();
1838 
1839 
1840  if(color != MPI_UNDEFINED) {
1841  outcomm->remoteIndices().setNeighbours(tneighbors);
1842  outcomm->remoteIndices().template rebuild<true>();
1843 
1844  }
1845 
1846  // release the memory
1847  delete[] sendTo;
1848 
1849  if(verbose) {
1850  oocomm.communicator().barrier();
1851  if(oocomm.communicator().rank()==0)
1852  std::cout<<" Storing indexsets took "<<
1853  time.elapsed()<<std::endl;
1854  }
1855 
1856 #ifdef PERF_REPART
1857  // stop the time for step 4) and print the results
1858  t4=MPI_Wtime()-t4;
1859  tSum = t1 + t2 + t3 + t4;
1860  std::cout<<std::endl
1861  <<mype<<": WTime for step 1): "<<t1
1862  <<" 2): "<<t2
1863  <<" 3): "<<t3
1864  <<" 4): "<<t4
1865  <<" total: "<<tSum
1866  <<std::endl;
1867 #endif
1868 
1869  return color!=MPI_UNDEFINED;
1870 
1871  }
1872 #else
1873  template<class G, class P,class T1, class T2, class R>
1874  bool graphRepartition(const G& graph, P& oocomm, int nparts,
1875  std::shared_ptr<P>& outcomm,
1876  R& redistInf,
1877  bool v=false)
1878  {
1879  if(nparts!=oocomm.size())
1880  DUNE_THROW(NotImplemented, "only available for MPI programs");
1881  }
1882 
1883 
1884  template<class G, class P,class T1, class T2, class R>
1885  bool commGraphRepartition(const G& graph, P& oocomm, int nparts,
1886  std::shared_ptr<P>& outcomm,
1887  R& redistInf,
1888  bool v=false)
1889  {
1890  if(nparts!=oocomm.size())
1891  DUNE_THROW(NotImplemented, "only available for MPI programs");
1892  }
1893 #endif // HAVE_MPI
1894 } // end of namespace Dune
1895 #endif
Classes providing communication interfaces for overlapping Schwarz methods.
Provides classes for building the matrix graph.
int globalOwnerVertices
Definition: repartition.hh:173
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
void fillIndexSetHoles(const G &graph, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm)
Fills the holes in an index set.
Definition: repartition.hh:81
bool graphRepartition(const G &graph, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, Metis::idx_t nparts, std::shared_ptr< Dune::OwnerOverlapCopyCommunication< T1, T2 >> &outcomm, RedistributeInterface &redistInf, bool verbose=false)
execute a graph repartition for a giving graph and indexset.
Definition: repartition.hh:1233
void print_carray(S &os, T *array, std::size_t l)
Definition: repartition.hh:776
bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T *xadj, T *adjncy, bool checkSymmetry)
Definition: repartition.hh:783
bool commGraphRepartition(const M &mat, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, Metis::idx_t nparts, std::shared_ptr< Dune::OwnerOverlapCopyCommunication< T1, T2 >> &outcomm, RedistributeInterface &redistInf, bool verbose=false)
Definition: repartition.hh:827
bool buildCommunication(const G &graph, std::vector< int > &realparts, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, std::shared_ptr< Dune::OwnerOverlapCopyCommunication< T1, T2 >> &outcomm, RedistributeInterface &redistInf, bool verbose=false)
Definition: repartition.hh:1454
float real_t
Definition: repartition.hh:51
std::size_t idx_t
Definition: repartition.hh:61
@ owner
Definition: owneroverlapcopy.hh:59
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:172
const CollectiveCommunication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:297
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:460
void copyCopyToAll(const T &source, T &dest) const
Communicate values from copy data points to all other data points.
Definition: owneroverlapcopy.hh:326
Dune::GlobalLookupIndexSet< ParallelIndexSet > GlobalLookupIndexSet
The type of the reverse lookup of indices.
Definition: owneroverlapcopy.hh:454
void buildGlobalLookup()
Definition: owneroverlapcopy.hh:493
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:469
void copyOwnerToAll(const T &source, T &dest) const
Communicate values from owner data points to all other data points.
Definition: owneroverlapcopy.hh:309
void freeGlobalLookup()
Definition: owneroverlapcopy.hh:518
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:447
const GlobalLookupIndexSet & globalLookup() const
Definition: owneroverlapcopy.hh:524
The (undirected) graph of a matrix.
Definition: graph.hh:49
Definition: repartition.hh:258
void reserveSpaceForReceiveInterface(int proc, int size)
Definition: repartition.hh:282
void buildReceiveInterface(std::vector< std::pair< TG, int > > &indices)
Definition: repartition.hh:291
~RedistributeInterface()
Definition: repartition.hh:299
void setCommunicator(MPI_Comm comm)
Definition: repartition.hh:259
void buildSendInterface(const std::vector< int > &toPart, const IS &idxset)
Definition: repartition.hh:264
void addReceiveIndex(int proc, std::size_t idx)
Definition: repartition.hh:286