3 #ifndef DUNE_ISTL_REPARTITION_HH
4 #define DUNE_ISTL_REPARTITION_HH
20 #include <dune/common/timer.hh>
21 #include <dune/common/unused.hh>
22 #include <dune/common/enumset.hh>
23 #include <dune/common/stdstreams.hh>
24 #include <dune/common/parallel/mpitraits.hh>
25 #include <dune/common/parallel/communicator.hh>
26 #include <dune/common/parallel/indexset.hh>
27 #include <dune/common/parallel/indicessyncer.hh>
28 #include <dune/common/parallel/remoteindices.hh>
29 #include <dune/common/rangeutilities.hh>
48 #if HAVE_PARMETIS && defined(REALTYPEWIDTH)
54 #if HAVE_PARMETIS && defined(IDXTYPEWIDTH)
56 #elif HAVE_PARMETIS && defined(HAVE_SCOTCH_NUM_TYPE)
57 using idx_t = SCOTCH_Num;
80 template<
class G,
class T1,
class T2>
81 void fillIndexSetHoles(
const G& graph, Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm)
83 typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2>::ParallelIndexSet IndexSet;
84 typedef typename IndexSet::LocalIndex::Attribute Attribute;
86 IndexSet& indexSet = oocomm.indexSet();
87 const typename Dune::OwnerOverlapCopyCommunication<T1,T2>::GlobalLookupIndexSet& lookup =oocomm.globalLookup();
90 typedef typename G::ConstVertexIterator VertexIterator;
93 std::size_t sum=0, needed = graph.noVertices()-indexSet.size();
94 std::vector<std::size_t> neededall(oocomm.communicator().size(), 0);
96 MPI_Allgather(&needed, 1, MPITraits<std::size_t>::getType() , &(neededall[0]), 1, MPITraits<std::size_t>::getType(), oocomm.communicator());
97 for(
int i=0; i<oocomm.communicator().size(); ++i)
106 typedef typename IndexSet::const_iterator Iterator;
108 end = indexSet.end();
109 for(Iterator it = indexSet.begin(); it != end; ++it)
110 maxgi=std::max(maxgi,it->global());
115 maxgi=oocomm.communicator().max(maxgi);
118 for(
int i=0; i<oocomm.communicator().rank(); ++i)
119 maxgi=maxgi+neededall[i];
122 std::map<int,SLList<std::pair<T1,Attribute> > > globalIndices;
123 storeGlobalIndicesOfRemoteIndices(globalIndices, oocomm.remoteIndices());
124 indexSet.beginResize();
126 for(VertexIterator vertex = graph.begin(), vend=graph.end(); vertex != vend; ++vertex) {
127 const typename IndexSet::IndexPair* pair=lookup.pair(*vertex);
135 indexSet.endResize();
137 repairLocalIndexPointers(globalIndices, oocomm.remoteIndices(), indexSet);
139 oocomm.freeGlobalLookup();
140 oocomm.buildGlobalLookup();
142 std::cout<<
"Holes are filled!"<<std::endl;
143 std::cout<<oocomm.communicator().rank()<<
": "<<oocomm.indexSet()<<std::endl;
150 class ParmetisDuneIndexMap
153 template<
class Graph,
class OOComm>
154 ParmetisDuneIndexMap(
const Graph& graph,
const OOComm& com);
155 int toParmetis(
int i)
const
157 return duneToParmetis[i];
159 int toLocalParmetis(
int i)
const
161 return duneToParmetis[i]-base_;
163 int operator[](
int i)
const
165 return duneToParmetis[i];
167 int toDune(
int i)
const
169 return parmetisToDune[i];
171 std::vector<int>::size_type numOfOwnVtx()
const
173 return parmetisToDune.size();
179 int globalOwnerVertices;
182 std::vector<int> duneToParmetis;
183 std::vector<int> parmetisToDune;
185 std::vector<Metis::idx_t> vtxDist_;
188 template<
class G,
class OOComm>
189 ParmetisDuneIndexMap::ParmetisDuneIndexMap(
const G& graph,
const OOComm& oocomm)
190 : duneToParmetis(graph.noVertices(), -1), vtxDist_(oocomm.communicator().size()+1)
192 int npes=oocomm.communicator().size(), mype=oocomm.communicator().rank();
194 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
195 typedef typename OOComm::OwnerSet OwnerSet;
198 Iterator end = oocomm.indexSet().end();
199 for(Iterator index = oocomm.indexSet().begin(); index != end; ++index) {
200 if (OwnerSet::contains(index->local().attribute())) {
204 parmetisToDune.resize(numOfOwnVtx);
205 std::vector<int> globalNumOfVtx(npes);
207 MPI_Allgather(&numOfOwnVtx, 1, MPI_INT, &(globalNumOfVtx[0]), 1, MPI_INT, oocomm.communicator());
211 for(
int i=0; i<npes; i++) {
213 base += globalNumOfVtx[i];
215 vtxDist_[i+1] = vtxDist_[i] + globalNumOfVtx[i];
217 globalOwnerVertices=vtxDist_[npes];
221 typedef typename G::ConstVertexIterator VertexIterator;
223 std::cout << oocomm.communicator().rank()<<
" vtxDist: ";
224 for(
int i=0; i<= npes; ++i)
225 std::cout << vtxDist_[i]<<
" ";
226 std::cout<<std::endl;
233 VertexIterator vend = graph.end();
234 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
235 const typename OOComm::ParallelIndexSet::IndexPair* index=oocomm.globalLookup().pair(*vertex);
237 if (OwnerSet::contains(index->local().attribute())) {
239 parmetisToDune[base-base_]=index->local();
240 duneToParmetis[index->local()] = base++;
250 std::cout <<oocomm.communicator().rank()<<
": before ";
251 for(std::size_t i=0; i<duneToParmetis.size(); ++i)
252 std::cout<<duneToParmetis[i]<<
" ";
253 std::cout<<std::endl;
255 oocomm.copyOwnerToAll(duneToParmetis,duneToParmetis);
257 std::cout <<oocomm.communicator().rank()<<
": after ";
258 for(std::size_t i=0; i<duneToParmetis.size(); ++i)
259 std::cout<<duneToParmetis[i]<<
" ";
260 std::cout<<std::endl;
265 struct RedistributeInterface
268 void setCommunicator(MPI_Comm comm)
272 template<
class Flags,
class IS>
273 void buildSendInterface(
const std::vector<int>& toPart,
const IS& idxset)
275 std::map<int,int> sizes;
277 typedef typename IS::const_iterator IIter;
278 for(IIter i=idxset.begin(), end=idxset.end(); i!=end; ++i)
279 if(Flags::contains(i->local().attribute()))
280 ++sizes[toPart[i->local()]];
283 typedef std::map<int,int>::const_iterator MIter;
284 for(MIter i=sizes.begin(), end=sizes.end(); i!=end; ++i)
285 interfaces()[i->first].first.reserve(i->second);
288 typedef typename IS::const_iterator IIter;
289 for(IIter i=idxset.begin(), end=idxset.end(); i!=end; ++i)
290 if(Flags::contains(i->local().attribute()))
291 interfaces()[toPart[i->local()]].first.add(i->local());
294 void reserveSpaceForReceiveInterface(
int proc,
int size)
296 interfaces()[proc].second.reserve(size);
298 void addReceiveIndex(
int proc, std::size_t idx)
300 interfaces()[proc].second.add(idx);
302 template<
typename TG>
303 void buildReceiveInterface(std::vector<std::pair<TG,int> >& indices)
305 typedef typename std::vector<std::pair<TG,int> >::const_iterator VIter;
307 for(VIter idx=indices.begin(); idx!= indices.end(); ++idx) {
308 interfaces()[idx->second].second.add(i++);
312 ~RedistributeInterface()
329 void createSendBuf(std::vector<GI>& ownerVec, std::set<GI>& overlapVec, std::set<int>& neighbors,
char *sendBuf,
int buffersize, MPI_Comm comm) {
331 std::size_t s=ownerVec.size();
335 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
336 MPI_Pack(&(ownerVec[0]), s, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
337 s = overlapVec.size();
338 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
339 typedef typename std::set<GI>::iterator Iter;
340 for(Iter i=overlapVec.begin(), end= overlapVec.end(); i != end; ++i)
341 MPI_Pack(
const_cast<GI*
>(&(*i)), 1, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
344 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
345 typedef typename std::set<int>::iterator IIter;
347 for(IIter i=neighbors.begin(), end= neighbors.end(); i != end; ++i)
348 MPI_Pack(
const_cast<int*
>(&(*i)), 1, MPI_INT, sendBuf, buffersize, &pos, comm);
359 void saveRecvBuf(
char *recvBuf,
int bufferSize, std::vector<std::pair<GI,int> >& ownerVec,
360 std::set<GI>& overlapVec, std::set<int>& neighbors, RedistributeInterface& inf,
int from, MPI_Comm comm) {
364 MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
365 inf.reserveSpaceForReceiveInterface(from, size);
366 ownerVec.reserve(ownerVec.size()+size);
367 for(; size!=0; --size) {
369 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
370 ownerVec.push_back(std::make_pair(gi,from));
373 MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
374 typename std::set<GI>::iterator ipos = overlapVec.begin();
375 Dune::dverb <<
"unpacking "<<size<<
" overlap"<<std::endl;
376 for(; size!=0; --size) {
378 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
379 ipos=overlapVec.insert(ipos, gi);
382 MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
383 Dune::dverb <<
"unpacking "<<size<<
" neighbors"<<std::endl;
384 typename std::set<int>::iterator npos = neighbors.begin();
385 for(; size!=0; --size) {
387 MPI_Unpack(recvBuf, bufferSize, &pos, &n, 1, MPI_INT, comm);
388 npos=neighbors.insert(npos, n);
406 void getDomain(
const MPI_Comm& comm, T *part,
int numOfOwnVtx,
int nparts,
int *myDomain, std::vector<int> &domainMapping) {
408 MPI_Comm_size(comm, &npes);
409 MPI_Comm_rank(comm, &mype);
416 std::vector<int> domain(nparts, 0);
417 std::vector<int> assigned(npes, 0);
419 domainMapping.assign(domainMapping.size(), -1);
422 for (i=0; i<numOfOwnVtx; i++) {
426 std::vector<int> domainMatrix(npes * nparts, -1);
429 int *buf =
new int[nparts];
430 for (i=0; i<nparts; i++) {
432 domainMatrix[mype*nparts+i] = domain[i];
435 int src = (mype-1+npes)%npes;
436 int dest = (mype+1)%npes;
438 for (i=0; i<npes-1; i++) {
439 MPI_Sendrecv_replace(buf, nparts, MPI_INT, dest, 0, src, 0, comm, &status);
441 pe = ((mype-1-i)+npes)%npes;
442 for(j=0; j<nparts; j++) {
444 domainMatrix[pe*nparts+j] = buf[j];
452 int maxOccurance = 0;
454 std::set<std::size_t> unassigned;
456 for(i=0; i<nparts; i++) {
457 for(j=0; j<npes; j++) {
459 if (assigned[j]==0) {
460 if (maxOccurance < domainMatrix[j*nparts+i]) {
461 maxOccurance = domainMatrix[j*nparts+i];
469 domainMapping[i] = pe;
479 unassigned.insert(i);
484 typename std::vector<int>::iterator next_free = assigned.begin();
486 for(
typename std::set<std::size_t>::iterator domain = unassigned.begin(),
487 end = unassigned.end(); domain != end; ++domain)
489 next_free = std::find_if(next_free, assigned.end(), std::bind(std::less<int>(), std::placeholders::_1, 1));
490 assert(next_free != assigned.end());
491 domainMapping[*domain] = next_free-assigned.begin();
499 bool operator()(
const T& t1,
const T& t2)
const
517 void mergeVec(std::vector<std::pair<GI, int> >& ownerVec, std::set<GI>& overlapSet) {
519 typedef typename std::vector<std::pair<GI,int> >::const_iterator VIter;
522 if(ownerVec.size()>0)
524 VIter old=ownerVec.begin();
525 for(VIter i=old+1, end=ownerVec.end(); i != end; old=i++)
527 if(i->first==old->first)
529 std::cerr<<
"Value at indes"<<old-ownerVec.begin()<<
" is the same as at index "
530 <<i-ownerVec.begin()<<
" ["<<old->first<<
","<<old->second<<
"]==["
531 <<i->first<<
","<<i->second<<
"]"<<std::endl;
539 typedef typename std::set<GI>::iterator SIter;
540 VIter v=ownerVec.begin(), vend=ownerVec.end();
541 for(SIter s=overlapSet.begin(), send=overlapSet.end(); s!=send;)
543 while(v!=vend && v->first<*s) ++v;
544 if(v!=vend && v->first==*s) {
549 overlapSet.erase(tmp);
569 template<
class OwnerSet,
class Graph,
class IS,
class GI>
570 void getNeighbor(
const Graph& g, std::vector<int>& part,
571 typename Graph::VertexDescriptor vtx,
const IS& indexSet,
572 int toPe, std::set<GI>& neighbor, std::set<int>& neighborProcs) {
573 typedef typename Graph::ConstEdgeIterator Iter;
574 for(Iter edge=g.beginEdges(vtx), end=g.endEdges(vtx); edge!=end; ++edge)
576 const typename IS::IndexPair* pindex = indexSet.pair(edge.target());
578 if(part[pindex->local()]!=toPe || !OwnerSet::contains(pindex->local().attribute()))
581 neighbor.insert(pindex->global());
582 neighborProcs.insert(part[pindex->local()]);
587 template<
class T,
class I>
588 void my_push_back(std::vector<T>& ownerVec,
const I& index,
int proc)
590 DUNE_UNUSED_PARAMETER(proc);
591 ownerVec.push_back(index);
594 template<
class T,
class I>
595 void my_push_back(std::vector<std::pair<T,int> >& ownerVec,
const I& index,
int proc)
597 ownerVec.push_back(std::make_pair(index,proc));
600 void reserve(std::vector<T>&, RedistributeInterface&,
int)
603 void reserve(std::vector<std::pair<T,int> >& ownerVec, RedistributeInterface& redist,
int proc)
605 redist.reserveSpaceForReceiveInterface(proc, ownerVec.size());
626 template<
class OwnerSet,
class G,
class IS,
class T,
class GI>
627 void getOwnerOverlapVec(
const G& graph, std::vector<int>& part, IS& indexSet,
628 int myPe,
int toPe, std::vector<T>& ownerVec, std::set<GI>& overlapSet,
629 RedistributeInterface& redist, std::set<int>& neighborProcs) {
630 DUNE_UNUSED_PARAMETER(myPe);
632 typedef typename IS::const_iterator Iterator;
633 for(Iterator index = indexSet.begin(); index != indexSet.end(); ++index) {
635 if(OwnerSet::contains(index->local().attribute()))
637 if(part[index->local()]==toPe)
639 getNeighbor<OwnerSet>(graph, part, index->local(), indexSet,
640 toPe, overlapSet, neighborProcs);
641 my_push_back(ownerVec, index->global(), toPe);
645 reserve(ownerVec, redist, toPe);
656 template<
class F,
class IS>
657 inline bool isOwner(IS& indexSet,
int index) {
659 const typename IS::IndexPair* pindex=indexSet.pair(index);
662 return F::contains(pindex->local().attribute());
666 class BaseEdgeFunctor
669 BaseEdgeFunctor(
Metis::idx_t* adj,
const ParmetisDuneIndexMap& data)
670 : i_(), adj_(adj), data_(data)
674 void operator()(
const T& edge)
678 adj_[i_] = data_.toParmetis(edge.target());
689 const ParmetisDuneIndexMap& data_;
694 :
public BaseEdgeFunctor
696 EdgeFunctor(
Metis::idx_t* adj,
const ParmetisDuneIndexMap& data, std::size_t)
697 : BaseEdgeFunctor(adj, data)
707 template<
class G,
class V,
class E,
class VM,
class EM>
708 class EdgeFunctor<
Dune::Amg::PropertiesGraph<G,V,E,VM,EM> >
709 :
public BaseEdgeFunctor
712 EdgeFunctor(
Metis::idx_t* adj,
const ParmetisDuneIndexMap& data, std::size_t s)
713 : BaseEdgeFunctor(adj, data)
719 void operator()(
const T& edge)
721 weight_[index()]=edge.properties().depends() ? 3 : 1;
722 BaseEdgeFunctor::operator()(edge);
753 template<
class F,
class G,
class IS,
class EW>
754 void getAdjArrays(G& graph, IS& indexSet,
Metis::idx_t *xadj,
760 typedef typename G::ConstVertexIterator VertexIterator;
762 typedef typename IS::const_iterator Iterator;
764 VertexIterator vend = graph.end();
767 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
768 if (isOwner<F>(indexSet,*vertex)) {
770 typedef typename G::ConstEdgeIterator EdgeIterator;
771 EdgeIterator eend = vertex.end();
772 xadj[j] = ew.index();
774 for(EdgeIterator edge = vertex.begin(); edge != eend; ++edge) {
779 xadj[j] = ew.index();
783 template<
class G,
class T1,
class T2>
784 bool buildCommunication(
const G& graph, std::vector<int>& realparts,
785 Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm,
786 Dune::OwnerOverlapCopyCommunication<T1,T2>*& outcomm,
787 RedistributeInterface& redistInf,
790 #ifndef METIS_VER_MAJOR
795 Metis::idx_t *adjwgt,
int *wgtflag,
int *numflag,
int *nparts,
799 Metis::idx_t *adjwgt,
int *wgtflag,
int *numflag,
int *nparts,
803 #endif // HAVE_PARMETIS
805 template<
class S,
class T>
806 inline void print_carray(S& os, T* array, std::size_t l)
808 for(T *cur=array, *end=array+l; cur!=end; ++cur)
812 template<
class S,
class T>
813 inline bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T* xadj,
814 T* adjncy,
bool checkSymmetry)
819 if(xadj[vtx]>noEdges||xadj[vtx]<0) {
820 std::cerr <<
"Check graph: xadj["<<vtx<<
"]="<<xadj[vtx]<<
" (>"
821 <<noEdges<<
") out of range!"<<std::endl;
824 if(xadj[vtx+1]>noEdges||xadj[vtx+1]<0) {
825 std::cerr <<
"Check graph: xadj["<<vtx+1<<
"]="<<xadj[vtx+1]<<
" (>"
826 <<noEdges<<
") out of range!"<<std::endl;
831 if(adjncy[i]<0||((std::size_t)adjncy[i])>gnoVtx) {
832 std::cerr<<
" Edge "<<adjncy[i]<<
" out of range ["<<0<<
","<<noVtx<<
")"
842 for(
Metis::idx_t j=xadj[target]; j< xadj[target+1]; ++j)
846 std::cerr<<
"Edge ("<<target<<
","<<vtx<<
") "<<i<<
" time"<<std::endl;
855 template<
class M,
class T1,
class T2>
858 Dune::OwnerOverlapCopyCommunication<T1,T2>*& outcomm,
859 RedistributeInterface& redistInf,
862 if(verbose && oocomm.communicator().rank()==0)
863 std::cout<<
"Repartitioning from "<<oocomm.communicator().size()
864 <<
" to "<<nparts<<
" parts"<<std::endl;
866 int rank = oocomm.communicator().rank();
868 int* part =
new int[1];
881 int noNeighbours = oocomm.remoteIndices().neighbours();
882 typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2>::RemoteIndices RemoteIndices;
883 typedef typename RemoteIndices::const_iterator
886 for(NeighbourIterator n= oocomm.remoteIndices().begin(); n != oocomm.remoteIndices().end();
907 for(
int i=0; i<oocomm.communicator().size(); ++i)
909 vtxdist[oocomm.communicator().size()]=oocomm.communicator().size();
912 xadj[1]=noNeighbours;
942 typedef typename RemoteIndices::const_iterator NeighbourIterator;
954 for(NeighbourIterator n= oocomm.remoteIndices().begin(); n != oocomm.remoteIndices().end();
956 if(n->first != rank) {
964 assert(isValidGraph(vtxdist[rank+1]-vtxdist[rank],
965 vtxdist[oocomm.communicator().size()],
966 noNeighbours, xadj, adjncy,
false));
975 for(
int i=0; i<nparts; ++i)
976 tpwgts[i]=1.0/nparts;
977 MPI_Comm comm=oocomm.communicator();
979 Dune::dinfo<<rank<<
" vtxdist: ";
980 print_carray(Dune::dinfo, vtxdist, oocomm.communicator().size()+1);
981 Dune::dinfo<<std::endl<<rank<<
" xadj: ";
982 print_carray(Dune::dinfo, xadj, 2);
983 Dune::dinfo<<std::endl<<rank<<
" adjncy: ";
984 print_carray(Dune::dinfo, adjncy, noNeighbours);
987 Dune::dinfo<<std::endl<<rank<<
" vwgt: ";
988 print_carray(Dune::dinfo, vwgt, 1);
989 Dune::dinfo<<std::endl<<rank<<
" adwgt: ";
990 print_carray(Dune::dinfo, adjwgt, noNeighbours);
992 Dune::dinfo<<std::endl;
993 oocomm.communicator().barrier();
994 if(verbose && oocomm.communicator().rank()==0)
995 std::cout<<
"Creating comm graph took "<<time.elapsed()<<std::endl;
998 #ifdef PARALLEL_PARTITION
1001 int options[5] ={ 0,1,15,0,0};
1006 ParMETIS_V3_PartKway(vtxdist, xadj, adjncy,
1007 vwgt, adjwgt, &wgtflag,
1008 &numflag, &ncon, &nparts, tpwgts, &ubvec, options, &edgecut, part,
1010 if(verbose && oocomm.communicator().rank()==0)
1011 std::cout<<
"ParMETIS took "<<time.elapsed()<<std::endl;
1015 std::size_t gnoedges=0;
1017 noedges =
new int[oocomm.communicator().size()];
1018 Dune::dverb<<
"noNeighbours: "<<noNeighbours<<std::endl;
1020 MPI_Allgather(&noNeighbours,1,MPI_INT,noedges,1, MPI_INT,oocomm.communicator());
1022 if(verbose && oocomm.communicator().rank()==0)
1023 std::cout<<
"Gathering noedges took "<<time1.elapsed()<<std::endl;
1026 Metis::idx_t noVertices = vtxdist[oocomm.communicator().size()];
1038 std::size_t localNoVtx=vtxdist[rank+1]-vtxdist[rank];
1040 std::size_t gxadjlen = vtxdist[oocomm.communicator().size()]-vtxdist[0]+oocomm.communicator().size();
1043 Dune::dinfo<<
"noedges: ";
1044 print_carray(Dune::dinfo, noedges, oocomm.communicator().size());
1045 Dune::dinfo<<std::endl;
1046 displ =
new int[oocomm.communicator().size()];
1047 xdispl =
new int[oocomm.communicator().size()];
1048 noxs =
new int[oocomm.communicator().size()];
1049 vdispl =
new int[oocomm.communicator().size()];
1050 novs =
new int[oocomm.communicator().size()];
1052 for(
int i=0; i < oocomm.communicator().size(); ++i) {
1053 noxs[i]=vtxdist[i+1]-vtxdist[i]+1;
1054 novs[i]=vtxdist[i+1]-vtxdist[i];
1059 for(
int *xcurr = xdispl, *vcurr = vdispl, *end=vdispl+oocomm.communicator().size();
1060 vcurr!=end; ++vcurr, ++xcurr, ++so, ++offset) {
1062 *xcurr = offset + *so;
1068 for(
int *curr=noedges, *end=noedges+oocomm.communicator().size()-1;
1069 curr!=end; ++curr) {
1074 Dune::dinfo<<
"displ: ";
1075 print_carray(Dune::dinfo, displ, oocomm.communicator().size());
1076 Dune::dinfo<<std::endl;
1080 for(
int *curr=noedges, *end=noedges+oocomm.communicator().size();
1085 Dune::dinfo<<
"gxadjlen: "<<gxadjlen<<
" noVertices: "<<noVertices
1086 <<
" gnoedges: "<<gnoedges<<std::endl;
1096 if(verbose && oocomm.communicator().rank()==0)
1097 std::cout<<
"Preparing global graph took "<<time1.elapsed()<<std::endl;
1101 MPI_Allgatherv(xadj,2,MPITraits<Metis::idx_t>::getType(),
1102 gxadj,noxs,xdispl,MPITraits<Metis::idx_t>::getType(),
1104 MPI_Allgatherv(adjncy,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1105 gadjncy,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1108 MPI_Allgatherv(adjwgt,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1109 gadjwgt,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1111 MPI_Allgatherv(vwgt,localNoVtx,MPITraits<Metis::idx_t>::getType(),
1112 gvwgt,novs,vdispl,MPITraits<Metis::idx_t>::getType(),
1115 if(verbose && oocomm.communicator().rank()==0)
1116 std::cout<<
"Gathering global graph data took "<<time1.elapsed()<<std::endl;
1123 print_carray(Dune::dinfo, gxadj, gxadjlen);
1128 for(
int i=1; i<oocomm.communicator().size(); ++i) {
1130 int lprev = vtxdist[i]-vtxdist[i-1];
1131 int l = vtxdist[i+1]-vtxdist[i];
1133 assert((start+l+offset)-gxadj<=
static_cast<Metis::idx_t>(gxadjlen));
1134 increment = *(start-1);
1135 std::transform(start+offset, start+l+offset, start, std::bind(std::plus<Metis::idx_t>(), std::placeholders::_1, increment));
1137 Dune::dinfo<<std::endl<<
"shifted xadj:";
1138 print_carray(Dune::dinfo, gxadj, noVertices+1);
1139 Dune::dinfo<<std::endl<<
" gadjncy: ";
1140 print_carray(Dune::dinfo, gadjncy, gnoedges);
1142 Dune::dinfo<<std::endl<<
" gvwgt: ";
1143 print_carray(Dune::dinfo, gvwgt, noVertices);
1144 Dune::dinfo<<std::endl<<
"adjwgt: ";
1145 print_carray(Dune::dinfo, gadjwgt, gnoedges);
1146 Dune::dinfo<<std::endl;
1149 if(verbose && oocomm.communicator().rank()==0)
1150 std::cout<<
"Postprocesing global graph data took "<<time1.elapsed()<<std::endl;
1153 assert(isValidGraph(noVertices, noVertices, gnoedges,
1154 gxadj, gadjncy,
true));
1157 if(verbose && oocomm.communicator().rank()==0)
1158 std::cout<<
"Creating grah one 1 process took "<<time.elapsed()<<std::endl;
1160 #if METIS_VER_MAJOR >= 5
1163 METIS_SetDefaultOptions(moptions);
1164 moptions[METIS_OPTION_NUMBERING] = numflag;
1165 METIS_PartGraphRecursive(&noVertices, &ncon, gxadj, gadjncy, gvwgt, NULL, gadjwgt,
1166 &nparts, NULL, NULL, moptions, &edgecut, gpart);
1168 int options[5] = {0, 1, 1, 3, 3};
1170 METIS_PartGraphRecursive(&noVertices, gxadj, gadjncy, gvwgt, gadjwgt, &wgtflag,
1171 &numflag, &nparts, options, &edgecut, gpart);
1174 if(verbose && oocomm.communicator().rank()==0)
1175 std::cout<<
"METIS took "<<time.elapsed()<<std::endl;
1178 Dune::dinfo<<std::endl<<
"part:";
1179 print_carray(Dune::dinfo, gpart, noVertices);
1189 MPI_Scatter(gpart, 1, MPITraits<Metis::idx_t>::getType(), part, 1,
1190 MPITraits<Metis::idx_t>::getType(), 0, comm);
1214 Dune::dinfo<<
" repart "<<rank <<
" -> "<< part[0]<<std::endl;
1216 std::vector<int> realpart(
mat.
N(), part[0]);
1219 oocomm.copyOwnerToAll(realpart, realpart);
1221 if(verbose && oocomm.communicator().rank()==0)
1222 std::cout<<
"Scattering repartitioning took "<<time.elapsed()<<std::endl;
1226 oocomm.buildGlobalLookup(
mat.
N());
1228 fillIndexSetHoles(graph, oocomm);
1229 if(verbose && oocomm.communicator().rank()==0)
1230 std::cout<<
"Filling index set took "<<time.elapsed()<<std::endl;
1234 int noNeighbours=oocomm.remoteIndices().neighbours();
1235 noNeighbours = oocomm.communicator().sum(noNeighbours)
1236 / oocomm.communicator().size();
1237 if(oocomm.communicator().rank()==0)
1238 std::cout<<
"Average no neighbours was "<<noNeighbours<<std::endl;
1240 bool ret = buildCommunication(graph, realpart, oocomm, outcomm, redistInf,
1242 if(verbose && oocomm.communicator().rank()==0)
1243 std::cout<<
"Building index sets took "<<time.elapsed()<<std::endl;
1265 template<
class G,
class T1,
class T2>
1267 Dune::OwnerOverlapCopyCommunication<T1,T2>*& outcomm,
1268 RedistributeInterface& redistInf,
1273 MPI_Comm comm=oocomm.communicator();
1274 oocomm.buildGlobalLookup(graph.noVertices());
1275 fillIndexSetHoles(graph, oocomm);
1277 if(verbose && oocomm.communicator().rank()==0)
1278 std::cout<<
"Filling holes took "<<time.elapsed()<<std::endl;
1285 double t1=0.0, t2=0.0, t3=0.0, t4=0.0, tSum=0.0;
1290 int mype = oocomm.communicator().rank();
1292 assert(nparts<=
static_cast<Metis::idx_t>(oocomm.communicator().size()));
1312 typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2> OOComm;
1313 typedef typename OOComm::OwnerSet OwnerSet;
1318 ParmetisDuneIndexMap indexMap(graph,oocomm);
1320 for(std::size_t i=0; i < indexMap.numOfOwnVtx(); ++i)
1324 if(oocomm.communicator().rank()==0 && nparts>1)
1325 std::cerr<<
"ParMETIS not activated. Will repartition to 1 domain instead of requested "
1326 <<nparts<<
" domains."<<std::endl;
1335 EdgeFunctor<G> ef(adjncy, indexMap, graph.noEdges());
1336 getAdjArrays<OwnerSet>(graph, oocomm.globalLookup(), xadj, ef);
1342 Metis::idx_t numflag=0, wgtflag=0, options[3], edgecut=0, ncon=1;
1345 for(
int i=0; i<nparts; ++i)
1346 tpwgts[i]=1.0/nparts;
1355 wgtflag = (ef.getWeights()!=NULL) ? 1 : 0;
1363 std::cout<<std::endl;
1364 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1365 <<options[1]<<
" "<<options[2]<<
"}, Ncon: "
1366 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1377 oocomm.communicator().barrier();
1378 if(oocomm.communicator().rank()==0)
1379 std::cout<<
"Preparing for parmetis took "<<time.elapsed()<<std::endl;
1386 ParMETIS_V3_PartKway(indexMap.vtxDist(), xadj, adjncy,
1387 NULL, ef.getWeights(), &wgtflag,
1388 &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &
const_cast<MPI_Comm&
>(comm));
1399 std::cout<<std::endl;
1400 std::cout<<
"ParMETIS_V3_PartKway reported a cut of "<<edgecut<<std::endl;
1401 std::cout<<std::endl;
1403 std::cout<<mype<<
": PARMETIS-Result: ";
1404 for(
int i=0; i < indexMap.vtxDist()[mype+1]-indexMap.vtxDist()[mype]; ++i) {
1405 std::cout<<part[i]<<
" ";
1407 std::cout<<std::endl;
1408 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1409 <<options[1]<<
" "<<options[2]<<
"}, Ncon: "
1410 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1421 oocomm.communicator().barrier();
1422 if(oocomm.communicator().rank()==0)
1423 std::cout<<
"Parmetis took "<<time.elapsed()<<std::endl;
1430 for(std::size_t i=0; i<indexMap.numOfOwnVtx(); ++i)
1440 std::vector<int> domainMapping(nparts);
1442 getDomain(comm, part, indexMap.numOfOwnVtx(), nparts, &myDomain, domainMapping);
1447 std::cout<<mype<<
": myDomain: "<<myDomain<<std::endl;
1448 std::cout<<mype<<
": DomainMapping: ";
1449 for(
auto j : range(nparts)) {
1450 std::cout<<
" do: "<<j<<
" pe: "<<domainMapping[j]<<
" ";
1452 std::cout<<std::endl;
1459 std::vector<int> setPartition(oocomm.indexSet().size(), -1);
1461 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
1463 for(Iterator index = oocomm.indexSet().begin(); index != oocomm.indexSet().end(); ++index)
1464 if(OwnerSet::contains(index->local().attribute())) {
1465 setPartition[index->local()]=domainMapping[part[i++]];
1469 oocomm.copyOwnerToAll(setPartition, setPartition);
1472 if (oocomm.getSolverCategory() ==
1474 oocomm.copyCopyToAll(setPartition, setPartition);
1475 bool ret = buildCommunication(graph, setPartition, oocomm, outcomm, redistInf,
1478 oocomm.communicator().barrier();
1479 if(oocomm.communicator().rank()==0)
1480 std::cout<<
"Creating indexsets took "<<time.elapsed()<<std::endl;
1487 template<
class G,
class T1,
class T2>
1488 bool buildCommunication(
const G& graph,
1489 std::vector<int>& setPartition, Dune::OwnerOverlapCopyCommunication<T1,T2>& oocomm,
1490 Dune::OwnerOverlapCopyCommunication<T1,T2>*& outcomm,
1491 RedistributeInterface& redistInf,
1494 typedef typename Dune::OwnerOverlapCopyCommunication<T1,T2> OOComm;
1495 typedef typename OOComm::OwnerSet OwnerSet;
1500 redistInf.buildSendInterface<OwnerSet>(setPartition, oocomm.indexSet());
1526 int npes = oocomm.communicator().size();
1529 std::set<int> recvFrom;
1535 typedef typename std::vector<int>::const_iterator VIter;
1536 int mype = oocomm.communicator().rank();
1539 std::set<int> tsendTo;
1540 for(VIter i=setPartition.begin(), iend = setPartition.end(); i!=iend; ++i)
1543 noSendTo = tsendTo.size();
1544 sendTo =
new int[noSendTo];
1545 typedef std::set<int>::const_iterator iterator;
1547 for(iterator i=tsendTo.begin(); i != tsendTo.end(); ++i, ++idx)
1552 int* gnoSend=
new int[oocomm.communicator().size()];
1553 int* gsendToDispl =
new int[oocomm.communicator().size()+1];
1555 MPI_Allgather(&noSendTo, 1, MPI_INT, gnoSend, 1,
1556 MPI_INT, oocomm.communicator());
1559 int totalNoRecv = 0;
1560 for(
int i=0; i<npes; ++i)
1561 totalNoRecv += gnoSend[i];
1563 int *gsendTo =
new int[totalNoRecv];
1567 for(
int i=0; i<npes; ++i)
1568 gsendToDispl[i+1]=gsendToDispl[i]+gnoSend[i];
1571 MPI_Allgatherv(sendTo, noSendTo, MPI_INT, gsendTo, gnoSend, gsendToDispl,
1572 MPI_INT, oocomm.communicator());
1575 for(
int proc=0; proc < npes; ++proc)
1576 for(
int i=gsendToDispl[proc]; i < gsendToDispl[proc+1]; ++i)
1577 if(gsendTo[i]==mype)
1578 recvFrom.insert(proc);
1580 bool existentOnNextLevel = recvFrom.size()>0;
1584 delete[] gsendToDispl;
1589 if(recvFrom.size()) {
1590 std::cout<<mype<<
": recvFrom: ";
1591 typedef typename std::set<int>::const_iterator siter;
1592 for(siter i=recvFrom.begin(); i!= recvFrom.end(); ++i) {
1597 std::cout<<std::endl<<std::endl;
1598 std::cout<<mype<<
": sendTo: ";
1599 for(
int i=0; i<noSendTo; i++) {
1600 std::cout<<sendTo[i]<<
" ";
1602 std::cout<<std::endl<<std::endl;
1606 if(oocomm.communicator().rank()==0)
1607 std::cout<<
" Communicating the receive information took "<<
1608 time.elapsed()<<std::endl;
1622 typedef typename OOComm::ParallelIndexSet::GlobalIndex GI;
1623 typedef std::vector<GI> GlobalVector;
1624 std::vector<std::pair<GI,int> > myOwnerVec;
1625 std::set<GI> myOverlapSet;
1626 GlobalVector sendOwnerVec;
1627 std::set<GI> sendOverlapSet;
1628 std::set<int> myNeighbors;
1633 char **sendBuffers=
new char*[noSendTo];
1634 MPI_Request *requests =
new MPI_Request[noSendTo];
1637 for(
int i=0; i < noSendTo; ++i) {
1639 sendOwnerVec.clear();
1640 sendOverlapSet.clear();
1643 std::set<int> neighbors;
1644 getOwnerOverlapVec<OwnerSet>(graph, setPartition, oocomm.globalLookup(),
1645 mype, sendTo[i], sendOwnerVec, sendOverlapSet, redistInf,
1651 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &buffersize);
1652 MPI_Pack_size(sendOwnerVec.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1654 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1656 MPI_Pack_size(sendOverlapSet.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1657 buffersize += tsize;
1658 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1659 buffersize += tsize;
1660 MPI_Pack_size(neighbors.size(), MPI_INT, oocomm.communicator(), &tsize);
1661 buffersize += tsize;
1663 sendBuffers[i] =
new char[buffersize];
1666 std::cout<<mype<<
" sending "<<sendOwnerVec.size()<<
" owner and "<<
1667 sendOverlapSet.size()<<
" overlap to "<<sendTo[i]<<
" buffersize="<<buffersize<<std::endl;
1669 createSendBuf(sendOwnerVec, sendOverlapSet, neighbors, sendBuffers[i], buffersize, oocomm.communicator());
1670 MPI_Issend(sendBuffers[i], buffersize, MPI_PACKED, sendTo[i], 99, oocomm.communicator(), requests+i);
1674 oocomm.communicator().barrier();
1675 if(oocomm.communicator().rank()==0)
1676 std::cout<<
" Creating sends took "<<
1677 time.elapsed()<<std::endl;
1682 int noRecv = recvFrom.size();
1683 int oldbuffersize=0;
1688 MPI_Probe(MPI_ANY_SOURCE, 99, oocomm.communicator(), &stat);
1690 MPI_Get_count(&stat, MPI_PACKED, &buffersize);
1692 if(oldbuffersize<buffersize) {
1695 recvBuf =
new char[buffersize];
1696 oldbuffersize = buffersize;
1698 MPI_Recv(recvBuf, buffersize, MPI_PACKED, stat.MPI_SOURCE, 99, oocomm.communicator(), &stat);
1699 saveRecvBuf(recvBuf, buffersize, myOwnerVec, myOverlapSet, myNeighbors, redistInf,
1700 stat.MPI_SOURCE, oocomm.communicator());
1709 MPI_Status *statuses =
new MPI_Status[noSendTo];
1710 int send = MPI_Waitall(noSendTo, requests, statuses);
1713 if(send==MPI_ERR_IN_STATUS) {
1714 std::cerr<<mype<<
": Error in sending :"<<std::endl;
1716 for(
int i=0; i< noSendTo; i++)
1717 if(statuses[i].MPI_ERROR!=MPI_SUCCESS) {
1720 MPI_Error_string(statuses[i].MPI_ERROR, message, &messageLength);
1721 std::cerr<<
" source="<<statuses[i].MPI_SOURCE<<
" message: ";
1722 for(
int j = 0; j < messageLength; j++)
1723 std::cout<<message[j];
1725 std::cerr<<std::endl;
1729 oocomm.communicator().barrier();
1730 if(oocomm.communicator().rank()==0)
1731 std::cout<<
" Receiving and saving took "<<
1732 time.elapsed()<<std::endl;
1736 for(
int i=0; i < noSendTo; ++i)
1737 delete[] sendBuffers[i];
1739 delete[] sendBuffers;
1743 redistInf.setCommunicator(oocomm.communicator());
1754 if (!existentOnNextLevel) {
1756 color= MPI_UNDEFINED;
1758 MPI_Comm outputComm;
1760 MPI_Comm_split(oocomm.communicator(), color, oocomm.communicator().rank(), &outputComm);
1761 outcomm =
new OOComm(outputComm,oocomm.getSolverCategory(),
true);
1764 int newrank=outcomm->communicator().rank();
1765 int *newranks=
new int[oocomm.communicator().size()];
1766 std::vector<int> tneighbors;
1767 tneighbors.reserve(myNeighbors.size());
1769 typename OOComm::ParallelIndexSet& outputIndexSet = outcomm->indexSet();
1771 MPI_Allgather(&newrank, 1, MPI_INT, newranks, 1,
1772 MPI_INT, oocomm.communicator());
1773 typedef typename std::set<int>::const_iterator IIter;
1776 std::cout<<oocomm.communicator().rank()<<
" ";
1777 for(IIter i=myNeighbors.begin(), end=myNeighbors.end();
1779 assert(newranks[*i]>=0);
1780 std::cout<<*i<<
"->"<<newranks[*i]<<
" ";
1781 tneighbors.push_back(newranks[*i]);
1783 std::cout<<std::endl;
1785 for(IIter i=myNeighbors.begin(), end=myNeighbors.end();
1787 tneighbors.push_back(newranks[*i]);
1791 myNeighbors.clear();
1794 oocomm.communicator().barrier();
1795 if(oocomm.communicator().rank()==0)
1796 std::cout<<
" Calculating new neighbours ("<<tneighbors.size()<<
") took "<<
1797 time.elapsed()<<std::endl;
1802 outputIndexSet.beginResize();
1805 std::sort(myOwnerVec.begin(), myOwnerVec.end(), SortFirst());
1809 typedef typename OOComm::ParallelIndexSet::LocalIndex LocalIndex;
1810 typedef typename std::vector<std::pair<GI,int> >::const_iterator VPIter;
1812 for(VPIter g=myOwnerVec.begin(), end =myOwnerVec.end(); g!=end; ++g, ++i ) {
1814 redistInf.addReceiveIndex(g->second, i);
1818 oocomm.communicator().barrier();
1819 if(oocomm.communicator().rank()==0)
1820 std::cout<<
" Adding owner indices took "<<
1821 time.elapsed()<<std::endl;
1830 mergeVec(myOwnerVec, myOverlapSet);
1834 myOwnerVec.swap(myOwnerVec);
1837 oocomm.communicator().barrier();
1838 if(oocomm.communicator().rank()==0)
1839 std::cout<<
" Merging indices took "<<
1840 time.elapsed()<<std::endl;
1846 typedef typename std::set<GI>::const_iterator SIter;
1847 for(SIter g=myOverlapSet.begin(), end=myOverlapSet.end(); g!=end; ++g, i++) {
1850 myOverlapSet.clear();
1851 outputIndexSet.endResize();
1853 #ifdef DUNE_ISTL_WITH_CHECKING
1855 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
1856 Iterator end = outputIndexSet.end();
1857 for(Iterator index = outputIndexSet.begin(); index != end; ++index) {
1858 if (OwnerSet::contains(index->local().attribute())) {
1862 numOfOwnVtx = oocomm.communicator().sum(numOfOwnVtx);
1869 Iterator index=outputIndexSet.begin();
1872 for(Iterator old = outputIndexSet.begin(); index != end; old=index++) {
1873 if(old->global()>index->global())
1874 DUNE_THROW(ISTLError,
"Index set's globalindex not sorted correctly");
1879 oocomm.communicator().barrier();
1880 if(oocomm.communicator().rank()==0)
1881 std::cout<<
" Adding overlap indices took "<<
1882 time.elapsed()<<std::endl;
1887 if(color != MPI_UNDEFINED) {
1888 outcomm->remoteIndices().setNeighbours(tneighbors);
1889 outcomm->remoteIndices().template rebuild<true>();
1897 oocomm.communicator().barrier();
1898 if(oocomm.communicator().rank()==0)
1899 std::cout<<
" Storing indexsets took "<<
1900 time.elapsed()<<std::endl;
1906 tSum = t1 + t2 + t3 + t4;
1907 std::cout<<std::endl
1908 <<mype<<
": WTime for step 1): "<<t1
1916 return color!=MPI_UNDEFINED;
1920 template<
class G,
class P,
class T1,
class T2,
class R>
1926 if(nparts!=oocomm.size())
1927 DUNE_THROW(NotImplemented,
"only available for MPI programs");
1931 template<
class G,
class P,
class T1,
class T2,
class R>
1937 if(nparts!=oocomm.size())
1938 DUNE_THROW(NotImplemented,
"only available for MPI programs");