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