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