Dune Core Modules (2.8.0)

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