5#ifndef DUNE_ISTL_REPARTITION_HH
6#define DUNE_ISTL_REPARTITION_HH
50#if HAVE_PARMETIS && defined(REALTYPEWIDTH)
56#if HAVE_PARMETIS && defined(IDXTYPEWIDTH)
57 using idx_t = ::idx_t;
58#elif HAVE_PARMETIS && defined(HAVE_SCOTCH_NUM_TYPE)
59 using idx_t = SCOTCH_Num;
63 using idx_t = std::size_t;
82 template<
class G,
class T1,
class T2>
86 typedef typename IndexSet::LocalIndex::Attribute
Attribute;
91 std::size_t sum=0, needed = graph.noVertices()-indexSet.
size();
92 std::vector<std::size_t> neededall(oocomm.communicator().
size(), 0);
95 for(
int i=0; i<oocomm.communicator().
size(); ++i)
104 auto end = indexSet.end();
105 for(
auto it = indexSet.begin(); it != end; ++it)
111 maxgi=oocomm.communicator().
max(maxgi);
114 for(
int i=0; i<oocomm.communicator().rank(); ++i)
115 maxgi=maxgi+neededall[i];
118 std::map<int,SLList<std::pair<T1,Attribute> > > globalIndices;
120 indexSet.beginResize();
123 const typename IndexSet::IndexPair* pair=lookup.
pair(*
vertex);
126 indexSet.add(maxgi,
typename IndexSet::LocalIndex(*
vertex, OwnerOverlapCopyAttributeSet::owner,
false));
131 indexSet.endResize();
135 oocomm.freeGlobalLookup();
136 oocomm.buildGlobalLookup();
138 std::cout<<
"Holes are filled!"<<std::endl;
139 std::cout<<oocomm.communicator().
rank()<<
": "<<oocomm.
indexSet()<<std::endl;
146 class ParmetisDuneIndexMap
149 template<
class Graph,
class OOComm>
150 ParmetisDuneIndexMap(
const Graph& graph,
const OOComm& com);
151 int toParmetis(
int i)
const
153 return duneToParmetis[i];
155 int toLocalParmetis(
int i)
const
157 return duneToParmetis[i]-base_;
159 int operator[](
int i)
const
161 return duneToParmetis[i];
163 int toDune(
int i)
const
165 return parmetisToDune[i];
167 std::vector<int>::size_type numOfOwnVtx()
const
169 return parmetisToDune.size();
171 Metis::idx_t* vtxDist()
175 int globalOwnerVertices;
178 std::vector<int> duneToParmetis;
179 std::vector<int> parmetisToDune;
181 std::vector<Metis::idx_t> vtxDist_;
184 template<
class G,
class OOComm>
185 ParmetisDuneIndexMap::ParmetisDuneIndexMap(
const G& graph,
const OOComm& oocomm)
186 : duneToParmetis(graph.noVertices(), -1), vtxDist_(oocomm.communicator().
size()+1)
188 int npes=oocomm.communicator().size(), mype=oocomm.communicator().rank();
190 typedef typename OOComm::OwnerSet OwnerSet;
193 auto end = oocomm.indexSet().end();
194 for(
auto index = oocomm.indexSet().begin(); index !=
end; ++index) {
195 if (OwnerSet::contains(index->local().attribute())) {
199 parmetisToDune.resize(numOfOwnVtx);
200 std::vector<int> globalNumOfVtx(npes);
202 MPI_Allgather(&numOfOwnVtx, 1, MPI_INT, &(globalNumOfVtx[0]), 1, MPI_INT, oocomm.communicator());
206 for(
int i=0; i<npes; i++) {
208 base += globalNumOfVtx[i];
210 vtxDist_[i+1] = vtxDist_[i] + globalNumOfVtx[i];
212 globalOwnerVertices=vtxDist_[npes];
216 std::cout << oocomm.communicator().rank()<<
" vtxDist: ";
217 for(
int i=0; i<= npes; ++i)
218 std::cout << vtxDist_[i]<<
" ";
219 std::cout<<std::endl;
226 auto vend = graph.end();
228 const typename OOComm::ParallelIndexSet::IndexPair* index=oocomm.globalLookup().pair(*
vertex);
230 if (OwnerSet::contains(index->local().attribute())) {
232 parmetisToDune[base-base_]=index->local();
233 duneToParmetis[index->local()] = base++;
243 std::cout <<oocomm.communicator().rank()<<
": before ";
244 for(std::size_t i=0; i<duneToParmetis.size(); ++i)
245 std::cout<<duneToParmetis[i]<<
" ";
246 std::cout<<std::endl;
248 oocomm.copyOwnerToAll(duneToParmetis,duneToParmetis);
250 std::cout <<oocomm.communicator().rank()<<
": after ";
251 for(std::size_t i=0; i<duneToParmetis.size(); ++i)
252 std::cout<<duneToParmetis[i]<<
" ";
253 std::cout<<std::endl;
258 struct RedistributeInterface
261 void setCommunicator(MPI_Comm comm)
265 template<
class Flags,
class IS>
266 void buildSendInterface(
const std::vector<int>& toPart,
const IS& idxset)
268 std::map<int,int> sizes;
270 for(
auto i=idxset.begin(), end=idxset.end(); i!=end; ++i)
271 if(Flags::contains(i->local().attribute()))
272 ++sizes[toPart[i->local()]];
275 for(
auto i=sizes.begin(), end=sizes.end(); i!=end; ++i)
276 interfaces()[i->first].first.reserve(i->second);
279 for(
auto i=idxset.begin(), end=idxset.end(); i!=end; ++i)
280 if(Flags::contains(i->local().attribute()))
281 interfaces()[toPart[i->local()]].first.add(i->local());
284 void reserveSpaceForReceiveInterface(
int proc,
int size)
288 void addReceiveIndex(
int proc, std::size_t idx)
292 template<
typename TG>
293 void buildReceiveInterface(std::vector<std::pair<TG,int> >& indices)
296 for(
auto idx=indices.begin(); idx!= indices.end(); ++idx) {
315 void createSendBuf(std::vector<GI>& ownerVec, std::set<GI>& overlapVec, std::set<int>& neighbors,
char *sendBuf,
int buffersize, MPI_Comm comm) {
317 std::size_t s=ownerVec.size();
321 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
322 MPI_Pack(&(ownerVec[0]), s, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
323 s = overlapVec.size();
324 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
325 for(
auto i=overlapVec.begin(), end= overlapVec.end(); i != end; ++i)
326 MPI_Pack(
const_cast<GI*
>(&(*i)), 1, MPITraits<GI>::getType(), sendBuf, buffersize, &pos, comm);
329 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
331 for(
auto i=neighbors.begin(), end= neighbors.end(); i != end; ++i)
332 MPI_Pack(
const_cast<int*
>(&(*i)), 1, MPI_INT, sendBuf, buffersize, &pos, comm);
343 void saveRecvBuf(
char *recvBuf,
int bufferSize, std::vector<std::pair<GI,int> >& ownerVec,
344 std::set<GI>& overlapVec, std::set<int>& neighbors, RedistributeInterface& inf,
int from, MPI_Comm comm) {
348 MPI_Unpack(recvBuf, bufferSize, &pos, &
size, 1, MPITraits<std::size_t>::getType(), comm);
349 inf.reserveSpaceForReceiveInterface(from,
size);
350 ownerVec.reserve(ownerVec.size()+
size);
353 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
354 ownerVec.push_back(std::make_pair(gi,from));
357 MPI_Unpack(recvBuf, bufferSize, &pos, &
size, 1, MPITraits<std::size_t>::getType(), comm);
358 typename std::set<GI>::iterator ipos = overlapVec.begin();
362 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
363 ipos=overlapVec.insert(ipos, gi);
366 MPI_Unpack(recvBuf, bufferSize, &pos, &
size, 1, MPITraits<std::size_t>::getType(), comm);
368 typename std::set<int>::iterator npos = neighbors.begin();
371 MPI_Unpack(recvBuf, bufferSize, &pos, &n, 1, MPI_INT, comm);
372 npos=neighbors.insert(npos, n);
390 void getDomain(
const MPI_Comm& comm, T *part,
int numOfOwnVtx,
int nparts,
int *myDomain, std::vector<int> &domainMapping) {
392 MPI_Comm_size(comm, &npes);
393 MPI_Comm_rank(comm, &mype);
398 std::vector<int> domain(nparts, 0);
399 std::vector<int> assigned(npes, 0);
401 domainMapping.assign(domainMapping.size(), -1);
404 for (
int i = 0; i < numOfOwnVtx; i++) {
408 std::vector<int> domainMatrix(npes * nparts, -1);
411 int *buf =
new int[nparts];
412 for (
int i = 0; i < nparts; i++) {
414 domainMatrix[mype*nparts+i] = domain[i];
417 int src = (mype-1+npes)%npes;
418 int dest = (mype+1)%npes;
420 for (
int i = 0; i < npes-1; i++) {
421 MPI_Sendrecv_replace(buf, nparts, MPI_INT, dest, 0, src, 0, comm, &status);
423 pe = ((mype-1-i)+npes)%npes;
424 for(
int j = 0; j < nparts; j++) {
426 domainMatrix[pe*nparts+j] = buf[j];
434 int maxOccurance = 0;
436 std::set<std::size_t> unassigned;
438 for (
int i = 0; i < nparts; i++) {
439 for (
int j = 0; j < npes; j++) {
441 if (assigned[j]==0) {
442 if (maxOccurance < domainMatrix[j*nparts+i]) {
443 maxOccurance = domainMatrix[j*nparts+i];
451 domainMapping[i] = pe;
461 unassigned.insert(i);
466 typename std::vector<int>::iterator next_free = assigned.begin();
468 for(
auto udomain = unassigned.begin(),
469 end = unassigned.end(); udomain != end; ++udomain)
471 next_free = std::find_if(next_free, assigned.end(), std::bind(std::less<int>(),
std::placeholders::_1, 1));
472 assert(next_free != assigned.end());
473 domainMapping[*udomain] = next_free-assigned.begin();
481 bool operator()(
const T& t1,
const T& t2)
const
499 void mergeVec(std::vector<std::pair<GI, int> >& ownerVec, std::set<GI>& overlapSet) {
503 if(ownerVec.size()>0)
505 auto old=ownerVec.begin();
506 for(
auto i=old+1, end=ownerVec.end(); i != end; old=i++)
508 if(i->first==old->first)
510 std::cerr<<
"Value at index "<<old-ownerVec.begin()<<
" is the same as at index "
511 <<i-ownerVec.begin()<<
" ["<<old->first<<
","<<old->second<<
"]==["
512 <<i->first<<
","<<i->second<<
"]"<<std::endl;
520 auto v=ownerVec.begin(), vend=ownerVec.end();
521 for(
auto s=overlapSet.begin(), send=overlapSet.end(); s!=send;)
523 while(v!=vend && v->first<*s) ++v;
524 if(v!=vend && v->first==*s) {
529 overlapSet.erase(tmp);
549 template<
class OwnerSet,
class Graph,
class IS,
class GI>
550 void getNeighbor(
const Graph& g, std::vector<int>& part,
551 typename Graph::VertexDescriptor vtx,
const IS& indexSet,
552 int toPe, std::set<GI>& neighbor, std::set<int>& neighborProcs) {
553 for(
auto edge=g.beginEdges(vtx), end=g.endEdges(vtx); edge!=end; ++edge)
555 const typename IS::IndexPair* pindex = indexSet.pair(edge.target());
557 if(part[pindex->local()]!=toPe || !OwnerSet::contains(pindex->local().attribute()))
560 neighbor.insert(pindex->global());
561 neighborProcs.insert(part[pindex->local()]);
566 template<
class T,
class I>
567 void my_push_back(std::vector<T>& ownerVec,
const I& index, [[maybe_unused]]
int proc)
569 ownerVec.push_back(index);
572 template<
class T,
class I>
573 void my_push_back(std::vector<std::pair<T,int> >& ownerVec,
const I& index,
int proc)
575 ownerVec.push_back(std::make_pair(index,proc));
578 void reserve(std::vector<T>&, RedistributeInterface&,
int)
581 void reserve(std::vector<std::pair<T,int> >& ownerVec, RedistributeInterface& redist,
int proc)
583 redist.reserveSpaceForReceiveInterface(proc, ownerVec.size());
604 template<
class OwnerSet,
class G,
class IS,
class T,
class GI>
605 void getOwnerOverlapVec(
const G& graph, std::vector<int>& part, IS& indexSet,
606 [[maybe_unused]]
int myPe,
int toPe, std::vector<T>& ownerVec, std::set<GI>& overlapSet,
607 RedistributeInterface& redist, std::set<int>& neighborProcs) {
608 for(
auto index = indexSet.begin(); index != indexSet.end(); ++index) {
610 if(OwnerSet::contains(index->local().attribute()))
612 if(part[index->local()]==toPe)
614 getNeighbor<OwnerSet>(graph, part, index->local(), indexSet,
615 toPe, overlapSet, neighborProcs);
616 my_push_back(ownerVec, index->global(), toPe);
620 reserve(ownerVec, redist, toPe);
631 template<
class F,
class IS>
632 inline bool isOwner(IS& indexSet,
int index) {
634 const typename IS::IndexPair* pindex=indexSet.pair(index);
637 return F::contains(pindex->local().attribute());
641 class BaseEdgeFunctor
644 BaseEdgeFunctor(Metis::idx_t* adj,
const ParmetisDuneIndexMap& data)
645 : i_(), adj_(adj), data_(data)
649 void operator()(
const T& edge)
653 adj_[i_] = data_.toParmetis(edge.target());
664 const ParmetisDuneIndexMap& data_;
669 :
public BaseEdgeFunctor
671 EdgeFunctor(Metis::idx_t* adj,
const ParmetisDuneIndexMap& data, std::size_t)
672 : BaseEdgeFunctor(adj, data)
675 Metis::idx_t* getWeights()
682 template<
class G,
class V,
class E,
class VM,
class EM>
683 class EdgeFunctor<
Dune::Amg::PropertiesGraph<G,V,E,VM,EM> >
684 :
public BaseEdgeFunctor
687 EdgeFunctor(Metis::idx_t* adj,
const ParmetisDuneIndexMap& data, std::size_t s)
688 : BaseEdgeFunctor(adj, data)
690 weight_=
new Metis::idx_t[s];
694 void operator()(
const T& edge)
696 weight_[index()]=edge.properties().depends() ? 3 : 1;
697 BaseEdgeFunctor::operator()(edge);
699 Metis::idx_t* getWeights()
708 Metis::idx_t* weight_;
726 template<
class F,
class G,
class IS,
class EW>
727 void getAdjArrays(G& graph, IS& indexSet, Metis::idx_t *xadj,
731 auto vend = graph.end();
734 if (isOwner<F>(indexSet,*
vertex)) {
737 xadj[j] = ew.index();
739 for(
auto edge =
vertex.begin(); edge != eend; ++edge) {
744 xadj[j] = ew.index();
748 template<
class G,
class T1,
class T2>
749 bool buildCommunication(
const G& graph, std::vector<int>& realparts,
752 RedistributeInterface& redistInf,
755#ifndef METIS_VER_MAJOR
759 void METIS_PartGraphKway(
int *nvtxs, Metis::idx_t *xadj, Metis::idx_t *adjncy, Metis::idx_t *vwgt,
760 Metis::idx_t *adjwgt,
int *wgtflag,
int *numflag,
int *nparts,
761 int *options,
int *edgecut, Metis::idx_t *part);
763 void METIS_PartGraphRecursive(
int *nvtxs, Metis::idx_t *xadj, Metis::idx_t *adjncy, Metis::idx_t *vwgt,
764 Metis::idx_t *adjwgt,
int *wgtflag,
int *numflag,
int *nparts,
765 int *options,
int *edgecut, Metis::idx_t *part);
770 template<
class S,
class T>
771 inline void print_carray(S& os, T* array, std::size_t l)
773 for(T *cur=array, *end=array+l; cur!=end; ++cur)
777 template<
class S,
class T>
778 inline bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T* xadj,
779 T* adjncy,
bool checkSymmetry)
784 for(Metis::idx_t vtx=0; vtx<(Metis::idx_t)noVtx; ++vtx) {
785 if(
static_cast<S
>(xadj[vtx])>noEdges || signbit(xadj[vtx])) {
786 std::cerr <<
"Check graph: xadj["<<vtx<<
"]="<<xadj[vtx]<<
" (>"
787 <<noEdges<<
") out of range!"<<std::endl;
790 if(
static_cast<S
>(xadj[vtx+1])>noEdges || signbit(xadj[vtx+1])) {
791 std::cerr <<
"Check graph: xadj["<<vtx+1<<
"]="<<xadj[vtx+1]<<
" (>"
792 <<noEdges<<
") out of range!"<<std::endl;
796 for(Metis::idx_t i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
797 if(signbit(adjncy[i]) || ((std::size_t)adjncy[i])>gnoVtx) {
798 std::cerr<<
" Edge "<<adjncy[i]<<
" out of range ["<<0<<
","<<noVtx<<
")"
804 for(Metis::idx_t i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
805 Metis::idx_t target=adjncy[i];
808 for(Metis::idx_t j=xadj[target]; j< xadj[target+1]; ++j)
812 std::cerr<<
"Edge ("<<target<<
","<<vtx<<
") "<<i<<
" time"<<std::endl;
821 template<
class M,
class T1,
class T2>
825 RedistributeInterface& redistInf,
828 if(verbose && oocomm.communicator().
rank()==0)
829 std::cout<<
"Repartitioning from "<<oocomm.communicator().
size()
830 <<
" to "<<nparts<<
" parts"<<std::endl;
832 int rank = oocomm.communicator().
rank();
834 int* part =
new int[1];
837 Metis::idx_t* part =
new Metis::idx_t[1];
861 Metis::idx_t *xadj=
new Metis::idx_t[2];
862 Metis::idx_t *vtxdist=
new Metis::idx_t[oocomm.communicator().
size()+1];
863 Metis::idx_t *adjncy=
new Metis::idx_t[noNeighbours];
865 Metis::idx_t *vwgt = 0;
866 Metis::idx_t *adjwgt = 0;
870 for(
int i=0; i<oocomm.communicator().
size(); ++i)
872 vtxdist[oocomm.communicator().
size()]=oocomm.communicator().
size();
875 xadj[1]=noNeighbours;
906 Metis::idx_t* adjp=adjncy;
909 vwgt =
new Metis::idx_t[1];
912 adjwgt =
new Metis::idx_t[noNeighbours];
913 Metis::idx_t* adjwp=adjwgt;
918 if(n->first != rank) {
926 assert(isValidGraph(vtxdist[rank+1]-vtxdist[rank],
927 vtxdist[oocomm.communicator().
size()],
928 noNeighbours, xadj, adjncy,
false));
930 [[maybe_unused]] Metis::idx_t wgtflag=0;
931 Metis::idx_t numflag=0;
932 Metis::idx_t edgecut;
936 Metis::real_t *tpwgts =
new Metis::real_t[nparts];
937 for(
int i=0; i<nparts; ++i)
938 tpwgts[i]=1.0/nparts;
939 MPI_Comm comm=oocomm.communicator();
955 oocomm.communicator().
barrier();
956 if(verbose && oocomm.communicator().
rank()==0)
957 std::cout<<
"Creating comm graph took "<<time.elapsed()<<std::endl;
960#ifdef PARALLEL_PARTITION
961 Metis::real_t ubvec = 1.15;
963 int options[5] ={ 0,1,15,0,0};
968 ParMETIS_V3_PartKway(vtxdist, xadj, adjncy,
969 vwgt, adjwgt, &wgtflag,
970 &numflag, &ncon, &nparts, tpwgts, &ubvec, options, &edgecut, part,
972 if(verbose && oocomm.communicator().
rank()==0)
973 std::cout<<
"ParMETIS took "<<time.elapsed()<<std::endl;
977 std::size_t gnoedges=0;
979 noedges =
new int[oocomm.communicator().
size()];
980 Dune::dverb<<
"noNeighbours: "<<noNeighbours<<std::endl;
982 MPI_Allgather(&noNeighbours,1,MPI_INT,noedges,1, MPI_INT,oocomm.communicator());
984 if(verbose && oocomm.communicator().
rank()==0)
985 std::cout<<
"Gathering noedges took "<<time1.elapsed()<<std::endl;
988 Metis::idx_t noVertices = vtxdist[oocomm.communicator().
size()];
989 Metis::idx_t *gxadj = 0;
990 Metis::idx_t *gvwgt = 0;
991 Metis::idx_t *gadjncy = 0;
992 Metis::idx_t *gadjwgt = 0;
993 Metis::idx_t *gpart = 0;
1000 std::size_t localNoVtx=vtxdist[rank+1]-vtxdist[rank];
1002 std::size_t gxadjlen = vtxdist[oocomm.communicator().
size()]-vtxdist[0]+oocomm.communicator().
size();
1008 displ =
new int[oocomm.communicator().
size()];
1009 xdispl =
new int[oocomm.communicator().
size()];
1010 noxs =
new int[oocomm.communicator().
size()];
1011 vdispl =
new int[oocomm.communicator().
size()];
1012 novs =
new int[oocomm.communicator().
size()];
1014 for(
int i=0; i < oocomm.communicator().
size(); ++i) {
1015 noxs[i]=vtxdist[i+1]-vtxdist[i]+1;
1016 novs[i]=vtxdist[i+1]-vtxdist[i];
1019 Metis::idx_t *so= vtxdist;
1021 for(
int *xcurr = xdispl, *vcurr = vdispl, *end=vdispl+oocomm.communicator().
size();
1022 vcurr!=end; ++vcurr, ++xcurr, ++so, ++offset) {
1024 *xcurr = offset + *so;
1030 for(
int *curr=noedges, *end=noedges+oocomm.communicator().
size()-1;
1031 curr!=end; ++curr) {
1042 for(
int *curr=noedges, *end=noedges+oocomm.communicator().
size();
1047 Dune::dinfo<<
"gxadjlen: "<<gxadjlen<<
" noVertices: "<<noVertices
1048 <<
" gnoedges: "<<gnoedges<<std::endl;
1049 gxadj =
new Metis::idx_t[gxadjlen];
1050 gpart =
new Metis::idx_t[noVertices];
1052 gvwgt =
new Metis::idx_t[noVertices];
1053 gadjwgt =
new Metis::idx_t[gnoedges];
1055 gadjncy =
new Metis::idx_t[gnoedges];
1058 if(verbose && oocomm.communicator().
rank()==0)
1059 std::cout<<
"Preparing global graph took "<<time1.elapsed()<<std::endl;
1063 MPI_Allgatherv(xadj,2,MPITraits<Metis::idx_t>::getType(),
1064 gxadj,noxs,xdispl,MPITraits<Metis::idx_t>::getType(),
1066 MPI_Allgatherv(adjncy,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1067 gadjncy,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1070 MPI_Allgatherv(adjwgt,noNeighbours,MPITraits<Metis::idx_t>::getType(),
1071 gadjwgt,noedges,displ,MPITraits<Metis::idx_t>::getType(),
1073 MPI_Allgatherv(vwgt,localNoVtx,MPITraits<Metis::idx_t>::getType(),
1074 gvwgt,novs,vdispl,MPITraits<Metis::idx_t>::getType(),
1077 if(verbose && oocomm.communicator().
rank()==0)
1078 std::cout<<
"Gathering global graph data took "<<time1.elapsed()<<std::endl;
1088 Metis::idx_t increment = vtxdist[1];
1089 Metis::idx_t *start=gxadj+1;
1090 for(
int i=1; i<oocomm.communicator().
size(); ++i) {
1092 int lprev = vtxdist[i]-vtxdist[i-1];
1093 int l = vtxdist[i+1]-vtxdist[i];
1095 assert((start+l+offset)-gxadj<=
static_cast<Metis::idx_t
>(gxadjlen));
1096 increment = *(start-1);
1097 std::transform(start+offset, start+l+offset, start, std::bind(std::plus<Metis::idx_t>(),
std::placeholders::_1, increment));
1111 if(verbose && oocomm.communicator().
rank()==0)
1112 std::cout<<
"Postprocessing global graph data took "<<time1.elapsed()<<std::endl;
1115 assert(isValidGraph(noVertices, noVertices, gnoedges,
1116 gxadj, gadjncy,
true));
1119 if(verbose && oocomm.communicator().
rank()==0)
1120 std::cout<<
"Creating grah one 1 process took "<<time.elapsed()<<std::endl;
1122#if METIS_VER_MAJOR >= 5
1123 Metis::idx_t ncon = 1;
1124 Metis::idx_t moptions[METIS_NOPTIONS];
1125 METIS_SetDefaultOptions(moptions);
1126 moptions[METIS_OPTION_NUMBERING] = numflag;
1127 METIS_PartGraphRecursive(&noVertices, &ncon, gxadj, gadjncy, gvwgt, NULL, gadjwgt,
1128 &nparts, NULL, NULL, moptions, &edgecut, gpart);
1130 int options[5] = {0, 1, 1, 3, 3};
1132 METIS_PartGraphRecursive(&noVertices, gxadj, gadjncy, gvwgt, gadjwgt, &wgtflag,
1133 &numflag, &nparts, options, &edgecut, gpart);
1136 if(verbose && oocomm.communicator().
rank()==0)
1137 std::cout<<
"METIS took "<<time.elapsed()<<std::endl;
1151 MPI_Scatter(gpart, 1, MPITraits<Metis::idx_t>::getType(), part, 1,
1152 MPITraits<Metis::idx_t>::getType(), 0, comm);
1176 Dune::dinfo<<
" repart "<<rank <<
" -> "<< part[0]<<std::endl;
1178 std::vector<int> realpart(mat.
N(), part[0]);
1183 if(verbose && oocomm.communicator().
rank()==0)
1184 std::cout<<
"Scattering repartitioning took "<<time.elapsed()<<std::endl;
1188 oocomm.buildGlobalLookup(mat.
N());
1191 if(verbose && oocomm.communicator().
rank()==0)
1192 std::cout<<
"Filling index set took "<<time.elapsed()<<std::endl;
1197 noNeighbours = oocomm.communicator().
sum(noNeighbours)
1198 / oocomm.communicator().
size();
1199 if(oocomm.communicator().
rank()==0)
1200 std::cout<<
"Average no neighbours was "<<noNeighbours<<std::endl;
1202 bool ret = buildCommunication(graph, realpart, oocomm, outcomm, redistInf,
1204 if(verbose && oocomm.communicator().
rank()==0)
1205 std::cout<<
"Building index sets took "<<time.elapsed()<<std::endl;
1227 template<
class G,
class T1,
class T2>
1230 RedistributeInterface& redistInf,
1235 MPI_Comm comm=oocomm.communicator();
1236 oocomm.buildGlobalLookup(graph.noVertices());
1239 if(verbose && oocomm.communicator().
rank()==0)
1240 std::cout<<
"Filling holes took "<<time.
elapsed()<<std::endl;
1247 double t1=0.0, t2=0.0, t3=0.0, t4=0.0, tSum=0.0;
1252 int mype = oocomm.communicator().
rank();
1254 assert(nparts<=
static_cast<Metis::idx_t
>(oocomm.communicator().
size()));
1275 typedef typename OOComm::OwnerSet OwnerSet;
1280 ParmetisDuneIndexMap indexMap(graph,oocomm);
1281 Metis::idx_t *part =
new Metis::idx_t[indexMap.numOfOwnVtx()];
1282 for(std::size_t i=0; i < indexMap.numOfOwnVtx(); ++i)
1286 if(oocomm.communicator().
rank()==0 && nparts>1)
1287 std::cerr<<
"ParMETIS not activated. Will repartition to 1 domain instead of requested "
1288 <<nparts<<
" domains."<<std::endl;
1295 Metis::idx_t *xadj =
new Metis::idx_t[indexMap.numOfOwnVtx()+1];
1296 Metis::idx_t *adjncy =
new Metis::idx_t[graph.noEdges()];
1297 EdgeFunctor<G> ef(adjncy, indexMap, graph.noEdges());
1298 getAdjArrays<OwnerSet>(graph, oocomm.globalLookup(), xadj, ef);
1304 Metis::idx_t numflag=0, wgtflag=0, options[3], edgecut=0, ncon=1;
1306 Metis::real_t *tpwgts =
new Metis::real_t[nparts];
1307 for(
int i=0; i<nparts; ++i)
1308 tpwgts[i]=1.0/nparts;
1309 Metis::real_t ubvec[1];
1317 wgtflag = (ef.getWeights()!=NULL) ? 1 : 0;
1325 std::cout<<std::endl;
1326 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1327 <<options[1]<<
" "<<options[2]<<
"}, Ncon: "
1328 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1339 oocomm.communicator().
barrier();
1340 if(oocomm.communicator().
rank()==0)
1341 std::cout<<
"Preparing for parmetis took "<<time.
elapsed()<<std::endl;
1348 ParMETIS_V3_PartKway(indexMap.vtxDist(), xadj, adjncy,
1349 NULL, ef.getWeights(), &wgtflag,
1350 &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &
const_cast<MPI_Comm&
>(comm));
1361 std::cout<<std::endl;
1362 std::cout<<
"ParMETIS_V3_PartKway reported a cut of "<<edgecut<<std::endl;
1363 std::cout<<std::endl;
1365 std::cout<<mype<<
": PARMETIS-Result: ";
1366 for(
int i=0; i < indexMap.vtxDist()[mype+1]-indexMap.vtxDist()[mype]; ++i) {
1367 std::cout<<part[i]<<
" ";
1369 std::cout<<std::endl;
1370 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {"
1371 <<options[1]<<
" "<<options[2]<<
"}, Ncon: "
1372 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1383 oocomm.communicator().
barrier();
1384 if(oocomm.communicator().
rank()==0)
1385 std::cout<<
"Parmetis took "<<time.
elapsed()<<std::endl;
1392 for(std::size_t i=0; i<indexMap.numOfOwnVtx(); ++i)
1402 std::vector<int> domainMapping(nparts);
1404 getDomain(comm, part, indexMap.numOfOwnVtx(), nparts, &myDomain, domainMapping);
1409 std::cout<<mype<<
": myDomain: "<<myDomain<<std::endl;
1410 std::cout<<mype<<
": DomainMapping: ";
1411 for(
auto j : range(nparts)) {
1412 std::cout<<
" do: "<<j<<
" pe: "<<domainMapping[j]<<
" ";
1414 std::cout<<std::endl;
1421 std::vector<int> setPartition(oocomm.
indexSet().
size(), -1);
1425 if(OwnerSet::contains(index->local().attribute())) {
1426 setPartition[index->local()]=domainMapping[part[i++]];
1436 bool ret = buildCommunication(graph, setPartition, oocomm, outcomm, redistInf,
1439 oocomm.communicator().
barrier();
1440 if(oocomm.communicator().
rank()==0)
1441 std::cout<<
"Creating indexsets took "<<time.
elapsed()<<std::endl;
1448 template<
class G,
class T1,
class T2>
1449 bool buildCommunication(
const G& graph,
1452 RedistributeInterface& redistInf,
1456 typedef typename OOComm::OwnerSet OwnerSet;
1461 redistInf.buildSendInterface<OwnerSet>(setPartition, oocomm.
indexSet());
1487 int npes = oocomm.communicator().
size();
1490 std::set<int> recvFrom;
1496 int mype = oocomm.communicator().
rank();
1499 std::set<int> tsendTo;
1500 for(
auto i=setPartition.begin(), iend = setPartition.end(); i!=iend; ++i)
1503 noSendTo = tsendTo.size();
1504 sendTo =
new int[noSendTo];
1506 for(
auto i=tsendTo.begin(); i != tsendTo.end(); ++i, ++idx)
1511 int* gnoSend=
new int[oocomm.communicator().
size()];
1512 int* gsendToDispl =
new int[oocomm.communicator().
size()+1];
1514 MPI_Allgather(&noSendTo, 1, MPI_INT, gnoSend, 1,
1515 MPI_INT, oocomm.communicator());
1518 int totalNoRecv = 0;
1519 for(
int i=0; i<npes; ++i)
1520 totalNoRecv += gnoSend[i];
1522 int *gsendTo =
new int[totalNoRecv];
1526 for(
int i=0; i<npes; ++i)
1527 gsendToDispl[i+1]=gsendToDispl[i]+gnoSend[i];
1530 MPI_Allgatherv(sendTo, noSendTo, MPI_INT, gsendTo, gnoSend, gsendToDispl,
1531 MPI_INT, oocomm.communicator());
1534 for(
int proc=0; proc < npes; ++proc)
1535 for(
int i=gsendToDispl[proc]; i < gsendToDispl[proc+1]; ++i)
1536 if(gsendTo[i]==mype)
1537 recvFrom.insert(proc);
1539 bool existentOnNextLevel = recvFrom.size()>0;
1543 delete[] gsendToDispl;
1548 if(recvFrom.size()) {
1549 std::cout<<mype<<
": recvFrom: ";
1550 for(
auto i=recvFrom.begin(); i!= recvFrom.end(); ++i) {
1555 std::cout<<std::endl<<std::endl;
1556 std::cout<<mype<<
": sendTo: ";
1557 for(
int i=0; i<noSendTo; i++) {
1558 std::cout<<sendTo[i]<<
" ";
1560 std::cout<<std::endl<<std::endl;
1564 if(oocomm.communicator().
rank()==0)
1565 std::cout<<
" Communicating the receive information took "<<
1566 time.elapsed()<<std::endl;
1580 typedef typename OOComm::ParallelIndexSet::GlobalIndex GI;
1581 typedef std::vector<GI> GlobalVector;
1582 std::vector<std::pair<GI,int> > myOwnerVec;
1583 std::set<GI> myOverlapSet;
1584 GlobalVector sendOwnerVec;
1585 std::set<GI> sendOverlapSet;
1586 std::set<int> myNeighbors;
1591 char **sendBuffers=
new char*[noSendTo];
1592 MPI_Request *requests =
new MPI_Request[noSendTo];
1595 for(
int i=0; i < noSendTo; ++i) {
1597 sendOwnerVec.clear();
1598 sendOverlapSet.clear();
1601 std::set<int> neighbors;
1602 getOwnerOverlapVec<OwnerSet>(graph, setPartition, oocomm.globalLookup(),
1603 mype, sendTo[i], sendOwnerVec, sendOverlapSet, redistInf,
1609 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &buffersize);
1610 MPI_Pack_size(sendOwnerVec.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1612 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1614 MPI_Pack_size(sendOverlapSet.size(), MPITraits<GI>::getType(), oocomm.communicator(), &tsize);
1615 buffersize += tsize;
1616 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.communicator(), &tsize);
1617 buffersize += tsize;
1618 MPI_Pack_size(neighbors.size(), MPI_INT, oocomm.communicator(), &tsize);
1619 buffersize += tsize;
1621 sendBuffers[i] =
new char[buffersize];
1624 std::cout<<mype<<
" sending "<<sendOwnerVec.size()<<
" owner and "<<
1625 sendOverlapSet.size()<<
" overlap to "<<sendTo[i]<<
" buffersize="<<buffersize<<std::endl;
1627 createSendBuf(sendOwnerVec, sendOverlapSet, neighbors, sendBuffers[i], buffersize, oocomm.communicator());
1628 MPI_Issend(sendBuffers[i], buffersize, MPI_PACKED, sendTo[i], 99, oocomm.communicator(), requests+i);
1632 oocomm.communicator().
barrier();
1633 if(oocomm.communicator().
rank()==0)
1634 std::cout<<
" Creating sends took "<<
1635 time.elapsed()<<std::endl;
1640 int noRecv = recvFrom.size();
1641 int oldbuffersize=0;
1646 MPI_Probe(MPI_ANY_SOURCE, 99, oocomm.communicator(), &stat);
1648 MPI_Get_count(&stat, MPI_PACKED, &buffersize);
1650 if(oldbuffersize<buffersize) {
1653 recvBuf =
new char[buffersize];
1654 oldbuffersize = buffersize;
1656 MPI_Recv(recvBuf, buffersize, MPI_PACKED, stat.MPI_SOURCE, 99, oocomm.communicator(), &stat);
1657 saveRecvBuf(recvBuf, buffersize, myOwnerVec, myOverlapSet, myNeighbors, redistInf,
1658 stat.MPI_SOURCE, oocomm.communicator());
1667 MPI_Status *statuses =
new MPI_Status[noSendTo];
1668 int send = MPI_Waitall(noSendTo, requests, statuses);
1671 if(send==MPI_ERR_IN_STATUS) {
1672 std::cerr<<mype<<
": Error in sending :"<<std::endl;
1674 for(
int i=0; i< noSendTo; i++)
1675 if(statuses[i].MPI_ERROR!=MPI_SUCCESS) {
1678 MPI_Error_string(statuses[i].MPI_ERROR, message, &messageLength);
1679 std::cerr<<
" source="<<statuses[i].MPI_SOURCE<<
" message: ";
1680 for(
int j = 0; j < messageLength; j++)
1681 std::cout<<message[j];
1683 std::cerr<<std::endl;
1687 oocomm.communicator().
barrier();
1688 if(oocomm.communicator().
rank()==0)
1689 std::cout<<
" Receiving and saving took "<<
1690 time.elapsed()<<std::endl;
1694 for(
int i=0; i < noSendTo; ++i)
1695 delete[] sendBuffers[i];
1697 delete[] sendBuffers;
1701 redistInf.setCommunicator(oocomm.communicator());
1712 if (!existentOnNextLevel) {
1714 color= MPI_UNDEFINED;
1716 MPI_Comm outputComm;
1718 MPI_Comm_split(oocomm.communicator(), color, oocomm.communicator().
rank(), &outputComm);
1722 int newrank=outcomm->communicator().rank();
1723 int *newranks=
new int[oocomm.communicator().
size()];
1724 std::vector<int> tneighbors;
1725 tneighbors.reserve(myNeighbors.size());
1727 typename OOComm::ParallelIndexSet& outputIndexSet = outcomm->indexSet();
1729 MPI_Allgather(&newrank, 1, MPI_INT, newranks, 1,
1730 MPI_INT, oocomm.communicator());
1733 std::cout<<oocomm.communicator().
rank()<<
" ";
1734 for(
auto i=myNeighbors.begin(), end=myNeighbors.end();
1736 assert(newranks[*i]>=0);
1737 std::cout<<*i<<
"->"<<newranks[*i]<<
" ";
1738 tneighbors.push_back(newranks[*i]);
1740 std::cout<<std::endl;
1742 for(
auto i=myNeighbors.begin(), end=myNeighbors.end();
1744 tneighbors.push_back(newranks[*i]);
1748 myNeighbors.clear();
1751 oocomm.communicator().
barrier();
1752 if(oocomm.communicator().
rank()==0)
1753 std::cout<<
" Calculating new neighbours ("<<tneighbors.size()<<
") took "<<
1754 time.elapsed()<<std::endl;
1759 outputIndexSet.beginResize();
1762 std::sort(myOwnerVec.begin(), myOwnerVec.end(), SortFirst());
1767 using LocalIndexT =
typename OOComm::ParallelIndexSet::LocalIndex;
1768 for(
auto g=myOwnerVec.begin(), end =myOwnerVec.end(); g!=end; ++g, ++i ) {
1769 outputIndexSet.add(g->first,LocalIndexT(i, OwnerOverlapCopyAttributeSet::owner,
true));
1770 redistInf.addReceiveIndex(g->second, i);
1774 oocomm.communicator().
barrier();
1775 if(oocomm.communicator().
rank()==0)
1776 std::cout<<
" Adding owner indices took "<<
1777 time.elapsed()<<std::endl;
1786 mergeVec(myOwnerVec, myOverlapSet);
1790 myOwnerVec.swap(myOwnerVec);
1793 oocomm.communicator().
barrier();
1794 if(oocomm.communicator().
rank()==0)
1795 std::cout<<
" Merging indices took "<<
1796 time.elapsed()<<std::endl;
1802 for(
auto g=myOverlapSet.begin(), end=myOverlapSet.end(); g!=end; ++g, i++) {
1803 outputIndexSet.add(*g,LocalIndexT(i, OwnerOverlapCopyAttributeSet::copy,
true));
1805 myOverlapSet.clear();
1806 outputIndexSet.endResize();
1808#ifdef DUNE_ISTL_WITH_CHECKING
1810 auto end = outputIndexSet.end();
1811 for(
auto index = outputIndexSet.begin(); index != end; ++index) {
1812 if (OwnerSet::contains(index->local().attribute())) {
1816 numOfOwnVtx = oocomm.communicator().
sum(numOfOwnVtx);
1823 std::is_sorted(outputIndexSet.begin(), outputIndexSet.end(),
1824 [](
const auto& v1,
const auto& v2){ return v1.global() < v2.global();});
1827 oocomm.communicator().
barrier();
1828 if(oocomm.communicator().
rank()==0)
1829 std::cout<<
" Adding overlap indices took "<<
1830 time.elapsed()<<std::endl;
1835 if(color != MPI_UNDEFINED) {
1836 outcomm->remoteIndices().setNeighbours(tneighbors);
1837 outcomm->remoteIndices().template rebuild<true>();
1845 oocomm.communicator().
barrier();
1846 if(oocomm.communicator().
rank()==0)
1847 std::cout<<
" Storing indexsets took "<<
1848 time.elapsed()<<std::endl;
1854 tSum = t1 + t2 + t3 + t4;
1855 std::cout<<std::endl
1856 <<mype<<
": WTime for step 1): "<<t1
1864 return color!=MPI_UNDEFINED;
1868 template<
class G,
class P,
class T1,
class T2,
class R>
1870 std::shared_ptr<P>& outcomm,
1874 if(nparts!=oocomm.size())
1875 DUNE_THROW(NotImplemented,
"only available for MPI programs");
1879 template<
class G,
class P,
class T1,
class T2,
class R>
1880 bool commGraphRepartition(
const G& graph, P& oocomm,
int nparts,
1881 std::shared_ptr<P>& outcomm,
1885 if(nparts!=oocomm.size())
1886 DUNE_THROW(NotImplemented,
"only available for MPI programs");
The (undirected) graph of a matrix.
Definition: graph.hh:51
T max(const T &in) const
Compute the maximum of the argument over all processes and return the result in every process....
Definition: mpicommunication.hh:257
int barrier() const
Wait until all processes have arrived at this point in the program.
Definition: mpicommunication.hh:272
int rank() const
Return rank, is between 0 and size()-1.
Definition: mpicommunication.hh:133
T sum(const T &in) const
Compute the sum of the argument over all processes and return the result in every process....
Definition: mpicommunication.hh:208
int size() const
Number of processes in set, is greater than 0.
Definition: mpicommunication.hh:139
Index Set Interface base class.
Definition: indexidset.hh:78
auto size(GeometryType type) const
Return total number of entities of given geometry type in entity set .
Definition: indexidset.hh:223
MPI_Comm communicator_
The MPI communicator we use.
Definition: interface.hh:325
size_type N() const
Return the number of rows.
Definition: matrix.hh:691
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:174
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:462
void copyCopyToAll(const T &source, T &dest) const
Communicate values from copy data points to all other data points.
Definition: owneroverlapcopy.hh:328
void copyOwnerToAll(const T &source, T &dest) const
Communicate values from owner data points to all other data points.
Definition: owneroverlapcopy.hh:311
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:471
LocalIndex::Attribute Attribute
The type of the attribute.
Definition: remoteindices.hh:227
A simple stop watch.
Definition: timer.hh:43
void reset() noexcept
Reset timer while keeping the running/stopped state.
Definition: timer.hh:57
double elapsed() const noexcept
Get elapsed user-time from last reset until now/last stop in seconds.
Definition: timer.hh:77
Provides utility classes for syncing distributed data via MPI communication.
Classes for building sets out of enumeration values.
Provides classes for building the matrix graph.
size_t size() const
Get the total number (public and nonpublic) indices.
void repairLocalIndexPointers(std::map< int, SLList< std::pair< typename T::GlobalIndex, typename T::LocalIndex::Attribute >, A > > &globalMap, RemoteIndices< T, A1 > &remoteIndices, const T &indexSet)
Repair the pointers to the local indices in the remote indices.
Definition: indicessyncer.hh:495
iterator begin()
Get an iterator over the indices positioned at the first index.
iterator end()
Get an iterator over the indices positioned after the last index.
const_iterator end() const
Get an iterator over all remote index lists.
Definition: remoteindices.hh:1528
const InformationMap & interfaces() const
Get information about the interfaces.
Definition: interface.hh:433
int neighbours() const
Get the number of processors we share indices with.
Definition: remoteindices.hh:1445
const IndexPair * pair(const std::size_t &local) const
Get the index pair corresponding to a local index.
void storeGlobalIndicesOfRemoteIndices(std::map< int, SLList< std::pair< typename T::GlobalIndex, typename T::LocalIndex::Attribute >, A > > &globalMap, const RemoteIndices< T, A1 > &remoteIndices)
Stores the corresponding global indices of the remote index information.
Definition: indicessyncer.hh:469
const_iterator begin() const
Get an iterator over all remote index lists.
Definition: remoteindices.hh:1521
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:55
typename FieldTraits< Type >::real_type real_t
Convenient access to FieldTraits<Type>::real_type.
Definition: typetraits.hh:301
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:141
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:117
Provides a map between global and local indices.
Class for adding missing indices of a distributed index set in a local communication.
Traits classes for mapping types onto MPI_Datatype.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
void fillIndexSetHoles(const G &graph, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm)
Fills the holes in an index set.
Definition: repartition.hh:83
bool graphRepartition(const G &graph, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, Metis::idx_t nparts, std::shared_ptr< Dune::OwnerOverlapCopyCommunication< T1, T2 > > &outcomm, RedistributeInterface &redistInf, bool verbose=false)
execute a graph repartition for a giving graph and indexset.
Definition: repartition.hh:1228
Classes providing communication interfaces for overlapping Schwarz methods.
Utilities for reduction like operations on ranges.
Classes describing a distributed indexset.
Standard Dune debug streams.
A traits class describing the mapping of types onto MPI_Datatypes.
Definition: mpitraits.hh:41
@ nonoverlapping
Category for non-overlapping solvers.
Definition: solvercategory.hh:27
static Category category(const OP &op, decltype(op.category()) *=nullptr)
Helperfunction to extract the solver category either from an enum, or from the newly introduced virtu...
Definition: solvercategory.hh:34