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