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