Dune Core Modules (2.4.2)

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