DUNE-FEM (unstable)

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