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> 57 template<
class G,
class T1,
class T2>
61 typedef typename IndexSet::LocalIndex::Attribute Attribute;
63 IndexSet& indexSet = oocomm.
indexSet();
67 typedef typename G::ConstVertexIterator VertexIterator;
70 std::size_t sum=0, needed = graph.noVertices()-indexSet.size();
71 std::vector<std::size_t> neededall(oocomm.
communicator().size(), 0);
73 MPI_Allgather(&needed, 1, MPITraits<std::size_t>::getType() , &(neededall[0]), 1, MPITraits<std::size_t>::getType(), oocomm.
communicator());
83 typedef typename IndexSet::const_iterator Iterator;
86 for(Iterator it = indexSet.begin(); it != end; ++it)
87 maxgi=std::max(maxgi,it->global());
96 maxgi=maxgi+neededall[i];
99 std::map<int,SLList<std::pair<T1,Attribute> > > globalIndices;
100 storeGlobalIndicesOfRemoteIndices(globalIndices, oocomm.
remoteIndices());
101 indexSet.beginResize();
103 for(VertexIterator vertex = graph.begin(), vend=graph.end(); vertex != vend; ++vertex) {
104 const typename IndexSet::IndexPair* pair=lookup.pair(*vertex);
112 indexSet.endResize();
114 repairLocalIndexPointers(globalIndices, oocomm.
remoteIndices(), indexSet);
119 std::cout<<
"Holes are filled!"<<std::endl;
127 class ParmetisDuneIndexMap
132 #if PARMETIS_MAJOR_VERSION > 3 133 typedef idx_t idxtype;
136 #endif // PARMETIS_MAJOR_VERSION > 3 138 typedef std::size_t idxtype;
139 #endif // #if HAVE_PARMETIS 141 template<
class Graph,
class OOComm>
142 ParmetisDuneIndexMap(
const Graph& graph,
const OOComm& com);
143 int toParmetis(
int i)
const 145 return duneToParmetis[i];
147 int toLocalParmetis(
int i)
const 149 return duneToParmetis[i]-base_;
151 int operator[](
int i)
const 153 return duneToParmetis[i];
155 int toDune(
int i)
const 157 return parmetisToDune[i];
159 std::vector<int>::size_type numOfOwnVtx()
const 161 return parmetisToDune.size();
170 std::vector<int> duneToParmetis;
171 std::vector<int> parmetisToDune;
173 std::vector<idxtype> vtxDist_;
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)
180 int npes=oocomm.communicator().size(), mype=oocomm.communicator().rank();
182 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
183 typedef typename OOComm::OwnerSet OwnerSet;
186 Iterator end = oocomm.indexSet().end();
187 for(Iterator index = oocomm.indexSet().begin(); index != end; ++index) {
188 if (OwnerSet::contains(index->local().attribute())) {
192 parmetisToDune.resize(numOfOwnVtx);
193 std::vector<int> globalNumOfVtx(npes);
195 MPI_Allgather(&numOfOwnVtx, 1, MPI_INT, &(globalNumOfVtx[0]), 1, MPI_INT, oocomm.communicator());
199 for(
int i=0; i<npes; i++) {
201 base += globalNumOfVtx[i];
203 vtxDist_[i+1] = vtxDist_[i] + globalNumOfVtx[i];
209 typedef typename G::ConstVertexIterator VertexIterator;
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;
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);
225 if (OwnerSet::contains(index->local().attribute())) {
227 parmetisToDune[base-base_]=index->local();
228 duneToParmetis[index->local()] = base++;
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;
243 oocomm.copyOwnerToAll(duneToParmetis,duneToParmetis);
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;
260 template<
class Flags,
class IS>
263 std::map<int,int> sizes;
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()]];
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);
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());
284 interfaces()[proc].second.reserve(size);
288 interfaces()[proc].second.add(idx);
290 template<
typename TG>
293 typedef typename std::vector<std::pair<TG,int> >::const_iterator VIter;
295 for(VIter idx=indices.begin(); idx!= indices.end(); ++idx) {
296 interfaces()[idx->second].second.add(i++);
308 typedef ParmetisDuneIndexMap :: idxtype idxtype ;
320 void createSendBuf(std::vector<GI>& ownerVec, std::set<GI>& overlapVec, std::set<int>& neighbors,
char *sendBuf,
int buffersize, MPI_Comm comm) {
322 std::size_t s=ownerVec.size();
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);
335 MPI_Pack(&s, 1, MPITraits<std::size_t>::getType(), sendBuf, buffersize, &pos, comm);
336 typedef typename std::set<int>::iterator IIter;
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);
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) {
355 MPI_Unpack(recvBuf, bufferSize, &pos, &size, 1, MPITraits<std::size_t>::getType(), comm);
357 ownerVec.reserve(ownerVec.size()+size);
358 for(; size!=0; --size) {
360 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
361 ownerVec.push_back(std::make_pair(gi,from));
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) {
369 MPI_Unpack(recvBuf, bufferSize, &pos, &gi, 1, MPITraits<GI>::getType(), comm);
370 ipos=overlapVec.insert(ipos, gi);
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) {
378 MPI_Unpack(recvBuf, bufferSize, &pos, &n, 1, MPI_INT, comm);
379 npos=neighbors.insert(npos, n);
397 void getDomain(
const MPI_Comm& comm, T *part,
int numOfOwnVtx,
int nparts,
int *myDomain, std::vector<int> &domainMapping) {
399 MPI_Comm_size(comm, &npes);
400 MPI_Comm_rank(comm, &mype);
407 std::vector<int> domain(nparts, 0);
408 std::vector<int> assigned(npes, 0);
410 domainMapping.assign(domainMapping.size(), -1);
413 for (i=0; i<numOfOwnVtx; i++) {
417 std::vector<int> domainMatrix(npes * nparts, -1);
420 int *buf =
new int[nparts];
421 for (i=0; i<nparts; i++) {
423 domainMatrix[mype*nparts+i] = domain[i];
426 int src = (mype-1+npes)%npes;
427 int dest = (mype+1)%npes;
429 for (i=0; i<npes-1; i++) {
430 MPI_Sendrecv_replace(buf, nparts, MPI_INT, dest, 0, src, 0, comm, &status);
432 pe = ((mype-1-i)+npes)%npes;
433 for(j=0; j<nparts; j++) {
435 domainMatrix[pe*nparts+j] = buf[j];
443 int maxOccurance = 0;
445 std::set<std::size_t> unassigned;
447 for(i=0; i<nparts; i++) {
448 for(j=0; j<npes; j++) {
450 if (assigned[j]==0) {
451 if (maxOccurance < domainMatrix[j*nparts+i]) {
452 maxOccurance = domainMatrix[j*nparts+i];
460 domainMapping[i] = pe;
470 unassigned.insert(i);
475 typename std::vector<int>::iterator next_free = assigned.begin();
477 for(
typename std::set<std::size_t>::iterator domain = unassigned.begin(),
478 end = unassigned.end(); domain != end; ++domain)
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();
490 bool operator()(
const T& t1,
const T& t2)
const 508 void mergeVec(std::vector<std::pair<GI, int> >& ownerVec, std::set<GI>& overlapSet) {
510 typedef typename std::vector<std::pair<GI,int> >::const_iterator VIter;
513 if(ownerVec.size()>0)
515 VIter old=ownerVec.begin();
516 for(VIter i=old+1, end=ownerVec.end(); i != end; old=i++)
518 if(i->first==old->first)
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;
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;)
534 while(v!=vend && v->first<*s) ++v;
535 if(v!=vend && v->first==*s) {
540 overlapSet.erase(tmp);
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)
567 const typename IS::IndexPair* pindex = indexSet.pair(edge.target());
569 if(part[pindex->local()]!=toPe || !OwnerSet::contains(pindex->local().attribute()))
572 neighbor.insert(pindex->global());
573 neighborProcs.insert(part[pindex->local()]);
578 template<
class T,
class I>
579 void my_push_back(std::vector<T>& ownerVec,
const I& index,
int proc)
581 DUNE_UNUSED_PARAMETER(proc);
582 ownerVec.push_back(index);
585 template<
class T,
class I>
586 void my_push_back(std::vector<std::pair<T,int> >& ownerVec,
const I& index,
int proc)
588 ownerVec.push_back(std::make_pair(index,proc));
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,
621 DUNE_UNUSED_PARAMETER(myPe);
623 typedef typename IS::const_iterator Iterator;
624 for(Iterator index = indexSet.begin(); index != indexSet.end(); ++index) {
626 if(OwnerSet::contains(index->local().attribute()))
628 if(part[index->local()]==toPe)
630 getNeighbor<OwnerSet>(graph, part, index->local(), indexSet,
631 toPe, overlapSet, neighborProcs);
632 my_push_back(ownerVec, index->global(), toPe);
636 reserve(ownerVec, redist, toPe);
647 template<
class F,
class IS>
648 inline bool isOwner(IS& indexSet,
int index) {
650 const typename IS::IndexPair* pindex=indexSet.pair(index);
653 return F::contains(pindex->local().attribute());
657 class BaseEdgeFunctor
660 BaseEdgeFunctor(idxtype* adj,
const ParmetisDuneIndexMap& data)
661 : i_(), adj_(adj), data_(data)
665 void operator()(
const T& edge)
669 adj_[i_] = data_.toParmetis(edge.target());
680 const ParmetisDuneIndexMap& data_;
685 :
public BaseEdgeFunctor
687 EdgeFunctor(idxtype* adj,
const ParmetisDuneIndexMap& data, std::size_t)
688 : BaseEdgeFunctor(adj, data)
691 idxtype* getWeights()
698 template<
class G,
class V,
class E,
class VM,
class EM>
700 :
public BaseEdgeFunctor
703 EdgeFunctor(idxtype* adj,
const ParmetisDuneIndexMap& data, std::size_t s)
704 : BaseEdgeFunctor(adj, data)
706 weight_=
new idxtype[s];
710 void operator()(
const T& edge)
712 weight_[index()]=edge.properties().depends() ? 3 : 1;
713 BaseEdgeFunctor::operator()(edge);
715 idxtype* getWeights()
744 template<
class F,
class G,
class IS,
class EW>
745 void getAdjArrays(G& graph, IS& indexSet, idxtype *xadj,
751 typedef typename G::ConstVertexIterator VertexIterator;
753 typedef typename IS::const_iterator Iterator;
755 VertexIterator vend = graph.end();
758 for(VertexIterator vertex = graph.begin(); vertex != vend; ++vertex) {
759 if (isOwner<F>(indexSet,*vertex)) {
761 typedef typename G::ConstEdgeIterator EdgeIterator;
762 EdgeIterator eend = vertex.end();
763 xadj[j] = ew.index();
765 for(EdgeIterator edge = vertex.begin(); edge != eend; ++edge) {
770 xadj[j] = ew.index();
774 template<
class G,
class T1,
class T2>
781 #ifndef METIS_VER_MAJOR 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);
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);
794 #endif // HAVE_PARMETIS 796 template<
class S,
class T>
799 for(T *cur=array, *end=array+l; cur!=end; ++cur)
803 template<
class S,
class T>
804 inline bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T* xadj,
805 T* adjncy,
bool checkSymmetry)
809 for(idxtype vtx=0; vtx<(idxtype)noVtx; ++vtx) {
810 if(xadj[vtx]>noEdges||xadj[vtx]<0) {
811 std::cerr <<
"Check graph: xadj["<<vtx<<
"]="<<xadj[vtx]<<
" (>" 812 <<noEdges<<
") out of range!"<<std::endl;
815 if(xadj[vtx+1]>noEdges||xadj[vtx+1]<0) {
816 std::cerr <<
"Check graph: xadj["<<vtx+1<<
"]="<<xadj[vtx+1]<<
" (>" 817 <<noEdges<<
") out of range!"<<std::endl;
821 for(idxtype i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
822 if(adjncy[i]<0||((std::size_t)adjncy[i])>gnoVtx) {
823 std::cerr<<
" Edge "<<adjncy[i]<<
" out of range ["<<0<<
","<<noVtx<<
")" 829 for(idxtype i=xadj[vtx]; i< xadj[vtx+1]; ++i) {
830 idxtype target=adjncy[i];
833 for(idxtype j=xadj[target]; j< xadj[target+1]; ++j)
837 std::cerr<<
"Edge ("<<target<<
","<<vtx<<
") "<<i<<
" time"<<std::endl;
846 template<
class M,
class T1,
class T2>
854 std::cout<<
"Repartitioning from "<<oocomm.
communicator().size()
855 <<
" to "<<nparts<<
" parts"<<std::endl;
859 int* part =
new int[1];
862 idxtype* part =
new idxtype[1];
874 typedef typename RemoteIndices::const_iterator
889 idxtype *xadj=
new idxtype[2];
890 idxtype *vtxdist=
new idxtype[oocomm.
communicator().size()+1];
891 idxtype *adjncy=
new idxtype[noNeighbours];
903 xadj[1]=noNeighbours;
933 typedef typename RemoteIndices::const_iterator NeighbourIterator;
935 idxtype* adjp=adjncy;
938 vwgt =
new idxtype[1];
941 adjwgt =
new idxtype[noNeighbours];
942 idxtype* adjwp=adjwgt;
947 if(n->first != rank) {
957 noNeighbours, xadj, adjncy,
false));
959 idxtype wgtflag=0, numflag=0;
964 float *tpwgts =
new float[nparts];
965 for(
int i=0; i<nparts; ++i)
966 tpwgts[i]=1.0/nparts;
967 int options[5] ={ 0,1,15,0,0};
970 Dune::dinfo<<rank<<
" vtxdist: ";
972 Dune::dinfo<<std::endl<<rank<<
" xadj: ";
974 Dune::dinfo<<std::endl<<rank<<
" adjncy: ";
978 Dune::dinfo<<std::endl<<rank<<
" vwgt: ";
980 Dune::dinfo<<std::endl<<rank<<
" adwgt: ";
983 Dune::dinfo<<std::endl;
986 std::cout<<
"Creating comm graph took "<<time.elapsed()<<std::endl;
989 #ifdef PARALLEL_PARTITION 996 ParMETIS_V3_PartKway(vtxdist, xadj, adjncy,
997 vwgt, adjwgt, &wgtflag,
998 &numflag, &ncon, &nparts, tpwgts, &ubvec, options, &edgecut, part,
1001 std::cout<<
"ParMETIS took "<<time.elapsed()<<std::endl;
1005 std::size_t gnoedges=0;
1008 Dune::dverb<<
"noNeighbours: "<<noNeighbours<<std::endl;
1010 MPI_Allgather(&noNeighbours,1,MPI_INT,noedges,1, MPI_INT,oocomm.
communicator());
1013 std::cout<<
"Gathering noedges took "<<time1.elapsed()<<std::endl;
1016 idxtype noVertices = vtxdist[oocomm.
communicator().size()];
1019 idxtype *gadjncy = 0;
1020 idxtype *gadjwgt = 0;
1028 std::size_t localNoVtx=vtxdist[rank+1]-vtxdist[rank];
1033 Dune::dinfo<<
"noedges: ";
1035 Dune::dinfo<<std::endl;
1043 noxs[i]=vtxdist[i+1]-vtxdist[i]+1;
1044 novs[i]=vtxdist[i+1]-vtxdist[i];
1047 idxtype *so= vtxdist;
1049 for(
int *xcurr = xdispl, *vcurr = vdispl, *end=vdispl+oocomm.
communicator().size();
1050 vcurr!=end; ++vcurr, ++xcurr, ++so, ++offset) {
1052 *xcurr = offset + *so;
1058 for(
int *curr=noedges, *end=noedges+oocomm.
communicator().size()-1;
1059 curr!=end; ++curr) {
1064 Dune::dinfo<<
"displ: ";
1066 Dune::dinfo<<std::endl;
1070 for(
int *curr=noedges, *end=noedges+oocomm.
communicator().size();
1075 Dune::dinfo<<
"gxadjlen: "<<gxadjlen<<
" noVertices: "<<noVertices
1076 <<
" gnoedges: "<<gnoedges<<std::endl;
1077 gxadj =
new idxtype[gxadjlen];
1078 gpart =
new idxtype[noVertices];
1080 gvwgt =
new idxtype[noVertices];
1081 gadjwgt =
new idxtype[gnoedges];
1083 gadjncy =
new idxtype[gnoedges];
1087 std::cout<<
"Preparing global graph took "<<time1.elapsed()<<std::endl;
1091 MPI_Allgatherv(xadj,2,MPITraits<idxtype>::getType(),
1092 gxadj,noxs,xdispl,MPITraits<idxtype>::getType(),
1094 MPI_Allgatherv(adjncy,noNeighbours,MPITraits<idxtype>::getType(),
1095 gadjncy,noedges,displ,MPITraits<idxtype>::getType(),
1098 MPI_Allgatherv(adjwgt,noNeighbours,MPITraits<idxtype>::getType(),
1099 gadjwgt,noedges,displ,MPITraits<idxtype>::getType(),
1101 MPI_Allgatherv(vwgt,localNoVtx,MPITraits<idxtype>::getType(),
1102 gvwgt,novs,vdispl,MPITraits<idxtype>::getType(),
1106 std::cout<<
"Gathering global graph data took "<<time1.elapsed()<<std::endl;
1116 idxtype increment = vtxdist[1];
1117 idxtype *start=gxadj+1;
1120 int lprev = vtxdist[i]-vtxdist[i-1];
1121 int l = vtxdist[i+1]-vtxdist[i];
1123 assert((start+l+offset)-gxadj<=static_cast<idxtype>(gxadjlen));
1124 increment = *(start-1);
1125 std::transform(start+offset, start+l+offset, start, std::bind2nd(std::plus<idxtype>(), increment));
1127 Dune::dinfo<<std::endl<<
"shifted xadj:";
1129 Dune::dinfo<<std::endl<<
" gadjncy: ";
1132 Dune::dinfo<<std::endl<<
" gvwgt: ";
1134 Dune::dinfo<<std::endl<<
"adjwgt: ";
1136 Dune::dinfo<<std::endl;
1140 std::cout<<
"Postprocesing global graph data took "<<time1.elapsed()<<std::endl;
1144 gxadj, gadjncy,
true));
1148 std::cout<<
"Creating grah one 1 process took "<<time.elapsed()<<std::endl;
1150 options[0]=0; options[1]=1; options[2]=1; options[3]=3; options[4]=3;
1151 #if METIS_VER_MAJOR >= 5 1153 idxtype moptions[METIS_NOPTIONS];
1154 METIS_SetDefaultOptions(moptions);
1155 moptions[METIS_OPTION_NUMBERING] = numflag;
1156 METIS_PartGraphRecursive(&noVertices, &ncon, gxadj, gadjncy, gvwgt, NULL, gadjwgt,
1157 &nparts, NULL, NULL, moptions, &edgecut, gpart);
1160 METIS_PartGraphRecursive(&noVertices, gxadj, gadjncy, gvwgt, gadjwgt, &wgtflag,
1161 &numflag, &nparts, options, &edgecut, gpart);
1165 std::cout<<
"METIS took "<<time.elapsed()<<std::endl;
1168 Dune::dinfo<<std::endl<<
"part:";
1179 MPI_Scatter(gpart, 1, MPITraits<idxtype>::getType(), part, 1,
1180 MPITraits<idxtype>::getType(), 0, comm);
1204 Dune::dinfo<<
" repart "<<rank <<
" -> "<< part[0]<<std::endl;
1206 std::vector<int> realpart(mat.N(), part[0]);
1212 std::cout<<
"Scattering repartitioning took "<<time.elapsed()<<std::endl;
1220 std::cout<<
"Filling index set took "<<time.elapsed()<<std::endl;
1228 std::cout<<
"Average no neighbours was "<<noNeighbours<<std::endl;
1233 std::cout<<
"Building index sets took "<<time.elapsed()<<std::endl;
1255 template<
class G,
class T1,
class T2>
1268 std::cout<<
"Filling holes took "<<time.elapsed()<<std::endl;
1275 double t1=0.0, t2=0.0, t3=0.0, t4=0.0, tSum=0.0;
1303 typedef typename OOComm::OwnerSet OwnerSet;
1308 ParmetisDuneIndexMap indexMap(graph,oocomm);
1310 idxtype *part =
new idxtype[indexMap.numOfOwnVtx()];
1312 std::size_t *part =
new std::size_t[indexMap.numOfOwnVtx()];
1314 for(std::size_t i=0; i < indexMap.numOfOwnVtx(); ++i)
1319 std::cerr<<
"ParMETIS not activated. Will repartition to 1 domain instead of requested " 1320 <<nparts<<
" domains."<<std::endl;
1327 idxtype *xadj =
new idxtype[indexMap.numOfOwnVtx()+1];
1328 idxtype *adjncy =
new idxtype[graph.noEdges()];
1329 EdgeFunctor<G> ef(adjncy, indexMap, graph.noEdges());
1330 getAdjArrays<OwnerSet>(graph, oocomm.
globalLookup(), xadj, ef);
1336 idxtype numflag=0, wgtflag=0, options[3], edgecut=0, ncon=1;
1338 float *tpwgts =
new float[nparts];
1339 for(
int i=0; i<nparts; ++i)
1340 tpwgts[i]=1.0/nparts;
1349 wgtflag = (ef.getWeights()!=NULL) ? 1 : 0;
1357 std::cout<<std::endl;
1358 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {" 1359 <<options[1]<<
" "<<options[2]<<
"}, Ncon: " 1360 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1373 std::cout<<
"Preparing for parmetis took "<<time.elapsed()<<std::endl;
1380 ParMETIS_V3_PartKway(indexMap.vtxDist(), xadj, adjncy,
1381 NULL, ef.getWeights(), &wgtflag,
1382 &numflag, &ncon, &nparts, tpwgts, ubvec, options, &edgecut, part, &
const_cast<MPI_Comm&
>(comm));
1393 std::cout<<std::endl;
1394 std::cout<<
"ParMETIS_V3_PartKway reported a cut of "<<edgecut<<std::endl;
1395 std::cout<<std::endl;
1397 std::cout<<mype<<
": PARMETIS-Result: ";
1398 for(
int i=0; i < indexMap.vtxDist()[mype+1]-indexMap.vtxDist()[mype]; ++i) {
1399 std::cout<<part[i]<<
" ";
1401 std::cout<<std::endl;
1402 std::cout<<
"Testing ParMETIS_V3_PartKway with options[1-2] = {" 1403 <<options[1]<<
" "<<options[2]<<
"}, Ncon: " 1404 <<ncon<<
", Nparts: "<<nparts<<std::endl;
1417 std::cout<<
"Parmetis took "<<time.elapsed()<<std::endl;
1424 for(std::size_t i=0; i<indexMap.numOfOwnVtx(); ++i)
1434 std::vector<int> domainMapping(nparts);
1436 getDomain(comm, part, indexMap.numOfOwnVtx(), nparts, &myDomain, domainMapping);
1441 std::cout<<mype<<
": myDomain: "<<myDomain<<std::endl;
1442 std::cout<<mype<<
": DomainMapping: ";
1443 for(
int j=0; j<nparts; j++) {
1444 std::cout<<
" do: "<<j<<
" pe: "<<domainMapping[j]<<
" ";
1446 std::cout<<std::endl;
1453 std::vector<int> setPartition(oocomm.
indexSet().size(), -1);
1455 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
1457 for(Iterator index = oocomm.
indexSet().begin(); index != oocomm.
indexSet().end(); ++index)
1458 if(OwnerSet::contains(index->local().attribute())) {
1459 setPartition[index->local()]=domainMapping[part[i++]];
1467 static_cast<int>(SolverCategory::nonoverlapping))
1474 std::cout<<
"Creating indexsets took "<<time.elapsed()<<std::endl;
1481 template<
class G,
class T1,
class T2>
1489 typedef typename OOComm::OwnerSet OwnerSet;
1523 std::set<int> recvFrom;
1529 typedef typename std::vector<int>::const_iterator VIter;
1533 std::set<int> tsendTo;
1534 for(VIter i=setPartition.begin(), iend = setPartition.end(); i!=iend; ++i)
1537 noSendTo = tsendTo.size();
1538 sendTo =
new int[noSendTo];
1539 typedef std::set<int>::const_iterator iterator;
1541 for(iterator i=tsendTo.begin(); i != tsendTo.end(); ++i, ++idx)
1547 int* gsendToDispl =
new int[oocomm.
communicator().size()+1];
1549 MPI_Allgather(&noSendTo, 1, MPI_INT, gnoSend, 1,
1553 int totalNoRecv = 0;
1554 for(
int i=0; i<npes; ++i)
1555 totalNoRecv += gnoSend[i];
1557 int *gsendTo =
new int[totalNoRecv];
1561 for(
int i=0; i<npes; ++i)
1562 gsendToDispl[i+1]=gsendToDispl[i]+gnoSend[i];
1565 MPI_Allgatherv(sendTo, noSendTo, MPI_INT, gsendTo, gnoSend, gsendToDispl,
1569 for(
int proc=0; proc < npes; ++proc)
1570 for(
int i=gsendToDispl[proc]; i < gsendToDispl[proc+1]; ++i)
1571 if(gsendTo[i]==mype)
1572 recvFrom.insert(proc);
1574 bool existentOnNextLevel = recvFrom.size()>0;
1578 delete[] gsendToDispl;
1583 if(recvFrom.size()) {
1584 std::cout<<mype<<
": recvFrom: ";
1585 typedef typename std::set<int>::const_iterator siter;
1586 for(siter i=recvFrom.begin(); i!= recvFrom.end(); ++i) {
1591 std::cout<<std::endl<<std::endl;
1592 std::cout<<mype<<
": sendTo: ";
1593 for(
int i=0; i<noSendTo; i++) {
1594 std::cout<<sendTo[i]<<
" ";
1596 std::cout<<std::endl<<std::endl;
1601 std::cout<<
" Communicating the receive information took "<<
1602 time.elapsed()<<std::endl;
1616 typedef typename OOComm::ParallelIndexSet::GlobalIndex GI;
1617 typedef std::vector<GI> GlobalVector;
1618 std::vector<std::pair<GI,int> > myOwnerVec;
1619 std::set<GI> myOverlapSet;
1620 GlobalVector sendOwnerVec;
1621 std::set<GI> sendOverlapSet;
1622 std::set<int> myNeighbors;
1627 char **sendBuffers=
new char*[noSendTo];
1628 MPI_Request *requests =
new MPI_Request[noSendTo];
1631 for(
int i=0; i < noSendTo; ++i) {
1633 sendOwnerVec.clear();
1634 sendOverlapSet.clear();
1637 std::set<int> neighbors;
1638 getOwnerOverlapVec<OwnerSet>(graph, setPartition, oocomm.
globalLookup(),
1639 mype, sendTo[i], sendOwnerVec, sendOverlapSet, redistInf,
1645 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.
communicator(), &buffersize);
1646 MPI_Pack_size(sendOwnerVec.size(), MPITraits<GI>::getType(), oocomm.
communicator(), &tsize);
1648 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.
communicator(), &tsize);
1650 MPI_Pack_size(sendOverlapSet.size(), MPITraits<GI>::getType(), oocomm.
communicator(), &tsize);
1651 buffersize += tsize;
1652 MPI_Pack_size(1, MPITraits<std::size_t>::getType(), oocomm.
communicator(), &tsize);
1653 buffersize += tsize;
1654 MPI_Pack_size(neighbors.size(), MPI_INT, oocomm.
communicator(), &tsize);
1655 buffersize += tsize;
1657 sendBuffers[i] =
new char[buffersize];
1660 std::cout<<mype<<
" sending "<<sendOwnerVec.size()<<
" owner and "<<
1661 sendOverlapSet.size()<<
" overlap to "<<sendTo[i]<<
" buffersize="<<buffersize<<std::endl;
1663 createSendBuf(sendOwnerVec, sendOverlapSet, neighbors, sendBuffers[i], buffersize, oocomm.
communicator());
1664 MPI_Issend(sendBuffers[i], buffersize, MPI_PACKED, sendTo[i], 99, oocomm.
communicator(), requests+i);
1670 std::cout<<
" Creating sends took "<<
1671 time.elapsed()<<std::endl;
1676 int noRecv = recvFrom.size();
1677 int oldbuffersize=0;
1682 MPI_Probe(MPI_ANY_SOURCE, 99, oocomm.
communicator(), &stat);
1684 MPI_Get_count(&stat, MPI_PACKED, &buffersize);
1686 if(oldbuffersize<buffersize) {
1689 recvBuf =
new char[buffersize];
1690 oldbuffersize = buffersize;
1692 MPI_Recv(recvBuf, buffersize, MPI_PACKED, stat.MPI_SOURCE, 99, oocomm.
communicator(), &stat);
1693 saveRecvBuf(recvBuf, buffersize, myOwnerVec, myOverlapSet, myNeighbors, redistInf,
1703 MPI_Status *statuses =
new MPI_Status[noSendTo];
1704 int send = MPI_Waitall(noSendTo, requests, statuses);
1707 if(send==MPI_ERR_IN_STATUS) {
1708 std::cerr<<mype<<
": Error in sending :"<<std::endl;
1710 for(
int i=0; i< noSendTo; i++)
1711 if(statuses[i].MPI_ERROR!=MPI_SUCCESS) {
1714 MPI_Error_string(statuses[i].MPI_ERROR, message, &messageLength);
1715 std::cerr<<
" source="<<statuses[i].MPI_SOURCE<<
" message: ";
1716 for(
int j = 0; j < messageLength; j++)
1717 std::cout<<message[j];
1719 std::cerr<<std::endl;
1725 std::cout<<
" Receiving and saving took "<<
1726 time.elapsed()<<std::endl;
1730 for(
int i=0; i < noSendTo; ++i)
1731 delete[] sendBuffers[i];
1733 delete[] sendBuffers;
1748 if (!existentOnNextLevel) {
1750 color= MPI_UNDEFINED;
1752 MPI_Comm outputComm;
1760 std::vector<int> tneighbors;
1761 tneighbors.reserve(myNeighbors.size());
1763 typename OOComm::ParallelIndexSet& outputIndexSet = outcomm->indexSet();
1765 MPI_Allgather(&newrank, 1, MPI_INT, newranks, 1,
1767 typedef typename std::set<int>::const_iterator IIter;
1771 for(IIter i=myNeighbors.begin(), end=myNeighbors.end();
1773 assert(newranks[*i]>=0);
1774 std::cout<<*i<<
"->"<<newranks[*i]<<
" ";
1775 tneighbors.push_back(newranks[*i]);
1777 std::cout<<std::endl;
1779 for(IIter i=myNeighbors.begin(), end=myNeighbors.end();
1781 tneighbors.push_back(newranks[*i]);
1785 myNeighbors.clear();
1790 std::cout<<
" Calculating new neighbours ("<<tneighbors.size()<<
") took "<<
1791 time.elapsed()<<std::endl;
1796 outputIndexSet.beginResize();
1799 std::sort(myOwnerVec.begin(), myOwnerVec.end(), SortFirst());
1803 typedef typename OOComm::ParallelIndexSet::LocalIndex LocalIndex;
1804 typedef typename std::vector<std::pair<GI,int> >::const_iterator VPIter;
1806 for(VPIter g=myOwnerVec.begin(), end =myOwnerVec.end(); g!=end; ++g, ++i ) {
1807 outputIndexSet.add(g->first,LocalIndex(i, OwnerOverlapCopyAttributeSet::owner,
true));
1814 std::cout<<
" Adding owner indices took "<<
1815 time.elapsed()<<std::endl;
1824 mergeVec(myOwnerVec, myOverlapSet);
1828 myOwnerVec.swap(myOwnerVec);
1833 std::cout<<
" Merging indices took "<<
1834 time.elapsed()<<std::endl;
1840 typedef typename std::set<GI>::const_iterator SIter;
1841 for(SIter g=myOverlapSet.begin(), end=myOverlapSet.end(); g!=end; ++g, i++) {
1842 outputIndexSet.add(*g,LocalIndex(i, OwnerOverlapCopyAttributeSet::copy,
true));
1844 myOverlapSet.clear();
1845 outputIndexSet.endResize();
1847 #ifdef DUNE_ISTL_WITH_CHECKING 1849 typedef typename OOComm::ParallelIndexSet::const_iterator Iterator;
1850 Iterator end = outputIndexSet.end();
1851 for(Iterator index = outputIndexSet.begin(); index != end; ++index) {
1852 if (OwnerSet::contains(index->local().attribute())) {
1863 Iterator index=outputIndexSet.begin();
1866 for(Iterator old = outputIndexSet.begin(); index != end; old=index++) {
1867 if(old->global()>index->global())
1868 DUNE_THROW(
ISTLError,
"Index set's globalindex not sorted correctly");
1875 std::cout<<
" Adding overlap indices took "<<
1876 time.elapsed()<<std::endl;
1881 if(color != MPI_UNDEFINED) {
1882 outcomm->remoteIndices().setNeighbours(tneighbors);
1883 outcomm->remoteIndices().template rebuild<true>();
1893 std::cout<<
" Storing indexsets took "<<
1894 time.elapsed()<<std::endl;
1900 tSum = t1 + t2 + t3 + t4;
1901 std::cout<<std::endl
1902 <<mype<<
": WTime for step 1): "<<t1
1910 return color!=MPI_UNDEFINED;
1914 template<
class G,
class P,
class T1,
class T2,
class R>
1920 if(nparts!=oocomm.size())
1921 DUNE_THROW(NotImplemented,
"only available for MPI programs");
1925 template<
class G,
class P,
class T1,
class T2,
class R>
1931 if(nparts!=oocomm.size())
1932 DUNE_THROW(NotImplemented,
"only available for MPI programs");
Attaches properties to the edges and vertices of a graph.
Definition: graph.hh:975
Definition: owneroverlapcopy.hh:59
Dune::GlobalLookupIndexSet< ParallelIndexSet > GlobalLookupIndexSet
The type of the reverse lookup of indices.
Definition: owneroverlapcopy.hh:457
const GlobalLookupIndexSet & globalLookup() const
Definition: owneroverlapcopy.hh:527
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:1256
void addReceiveIndex(int proc, std::size_t idx)
Definition: repartition.hh:286
Definition: example.cc:35
SolverCategory::Category getSolverCategory() const
Get Solver Category.
Definition: owneroverlapcopy.hh:299
The (undirected) graph of a matrix.
Definition: graph.hh:48
Dune::RemoteIndices< PIS > RemoteIndices
The type of the remote indices.
Definition: owneroverlapcopy.hh:453
void buildReceiveInterface(std::vector< std::pair< TG, int > > &indices)
Definition: repartition.hh:291
void copyOwnerToAll(const T &source, T &dest) const
Communicate values from owner data points to all other data points.
Definition: owneroverlapcopy.hh:315
void copyCopyToAll(const T &source, T &dest) const
Communicate values from copy data points to all other data points.
Definition: owneroverlapcopy.hh:332
const ParallelIndexSet & indexSet() const
Get the underlying parallel index set.
Definition: owneroverlapcopy.hh:463
void setCommunicator(MPI_Comm comm)
Definition: repartition.hh:256
void print_carray(S &os, T *array, std::size_t l)
Definition: repartition.hh:797
void buildGlobalLookup()
Definition: owneroverlapcopy.hh:496
Classes providing communication interfaces for overlapping Schwarz methods.
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: basearray.hh:19
bool isValidGraph(std::size_t noVtx, std::size_t gnoVtx, S noEdges, T *xadj, T *adjncy, bool checkSymmetry)
Definition: repartition.hh:804
const CollectiveCommunication< MPI_Comm > & communicator() const
Definition: owneroverlapcopy.hh:303
void fillIndexSetHoles(const G &graph, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm)
Fills the holes in an index set.
Definition: repartition.hh:58
Definition: repartition.hh:253
Provides classes for building the matrix graph.
void reserveSpaceForReceiveInterface(int proc, int size)
Definition: repartition.hh:282
~RedistributeInterface()
Definition: repartition.hh:300
A class setting up standard communication for a two-valued attribute set with owner/overlap/copy sema...
Definition: owneroverlapcopy.hh:171
Dune::ParallelIndexSet< GlobalIdType, LI, 512 > ParallelIndexSet
The type of the parallel index set.
Definition: owneroverlapcopy.hh:450
derive error class from the base class in common
Definition: istlexception.hh:16
int globalOwnerVertices
Definition: repartition.hh:167
bool commGraphRepartition(const M &mat, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, idxtype nparts, Dune::OwnerOverlapCopyCommunication< T1, T2 > *&outcomm, RedistributeInterface &redistInf, bool verbose=false)
Definition: repartition.hh:847
bool buildCommunication(const G &graph, std::vector< int > &realparts, Dune::OwnerOverlapCopyCommunication< T1, T2 > &oocomm, Dune::OwnerOverlapCopyCommunication< T1, T2 > *&outcomm, RedistributeInterface &redistInf, bool verbose=false)
Definition: repartition.hh:1482
void buildSendInterface(const std::vector< int > &toPart, const IS &idxset)
Definition: repartition.hh:261
const RemoteIndices & remoteIndices() const
Get the underlying remote indices.
Definition: owneroverlapcopy.hh:472
void freeGlobalLookup()
Definition: owneroverlapcopy.hh:521