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