46 #ifndef MUELU_REPARTITIONFACTORY_DEF_HPP 47 #define MUELU_REPARTITIONFACTORY_DEF_HPP 57 #include <Teuchos_CommHelpers.hpp> 59 #include <Xpetra_Map.hpp> 60 #include <Xpetra_MapFactory.hpp> 61 #include <Xpetra_MultiVectorFactory.hpp> 62 #include <Xpetra_VectorFactory.hpp> 63 #include <Xpetra_Import.hpp> 64 #include <Xpetra_ImportFactory.hpp> 65 #include <Xpetra_Export.hpp> 66 #include <Xpetra_ExportFactory.hpp> 67 #include <Xpetra_Matrix.hpp> 68 #include <Xpetra_MatrixFactory.hpp> 70 #include "MueLu_Utilities.hpp" 72 #include "MueLu_CloneRepartitionInterface.hpp" 80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 88 #undef SET_VALID_ENTRY 91 validParamList->
set<
RCP<const FactoryBase> >(
"number of partitions", Teuchos::null,
"Instance of RepartitionHeuristicFactory.");
94 return validParamList;
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 Input(currentLevel,
"A");
100 Input(currentLevel,
"number of partitions");
101 Input(currentLevel,
"Partition");
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 bool remapPartitions = pL.
get<
bool> (
"repartition: remap parts");
121 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
123 GO indexBase = rowMap->getIndexBase();
124 Xpetra::UnderlyingLib lib = rowMap->lib();
129 int myRank = comm->getRank();
130 int numProcs = comm->getSize();
137 int numPartitions = Get<int>(currentLevel,
"number of partitions");
142 RCP<GOVector> decomposition = Get<RCP<GOVector> >(currentLevel,
"Partition");
145 if(remapPartitions ==
true && Teuchos::rcp_dynamic_cast<const CloneRepartitionInterface>(GetFactory(
"Partition")) != Teuchos::null) {
149 remapPartitions =
false;
153 if (numPartitions == 1) {
158 GetOStream(
Warnings0) <<
"Only one partition: Skip call to the repartitioner." << std::endl;
159 }
else if (numPartitions == -1) {
161 GetOStream(
Warnings0) <<
"No repartitioning necessary: partitions were left unchanged by the repartitioner" << std::endl;
162 Set<RCP<const Import> >(currentLevel,
"Importer", Teuchos::null);
203 if (remapPartitions) {
206 DeterminePartitionPlacement(*A, *decomposition, numPartitions);
218 if (decomposition->getLocalLength() > 0)
219 decompEntries = decomposition->getData(0);
221 #ifdef HAVE_MUELU_DEBUG 223 int incorrectRank = -1;
224 for (
int i = 0; i < decompEntries.
size(); i++)
225 if (decompEntries[i] >= numProcs || decompEntries[i] < 0) {
226 incorrectRank = myRank;
230 int incorrectGlobalRank = -1;
236 myGIDs.
reserve(decomposition->getLocalLength());
241 typedef std::map<GO, Array<GO> > map_type;
243 for (LO i = 0; i < decompEntries.
size(); i++) {
244 GO
id = decompEntries[i];
245 GO GID = rowMap->getGlobalElement(i);
250 sendMap[id].push_back(GID);
252 decompEntries = Teuchos::null;
255 GO numLocalKept = myGIDs.
size(), numGlobalKept, numGlobalRows = A->getGlobalNumRows();
257 GetOStream(
Statistics2) <<
"Unmoved rows: " << numGlobalKept <<
" / " << numGlobalRows <<
" (" << 100*Teuchos::as<double>(numGlobalKept)/numGlobalRows <<
"%)" << std::endl;
260 int numSend = sendMap.size(), numRecv;
266 for (
typename map_type::const_iterator it = sendMap.begin(); it != sendMap.end(); it++)
267 myParts[cnt++] = it->first;
271 GO partsIndexBase = 0;
273 RCP<Map> partsIOwn = MapFactory ::Build(lib, numProcs, myPart(), partsIndexBase, comm);
274 RCP<Export> partsExport = ExportFactory::Build(partsIHave, partsIOwn);
276 RCP<GOVector> partsISend = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIHave);
277 RCP<GOVector> numPartsIRecv = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIOwn);
279 ArrayRCP<GO> partsISendData = partsISend->getDataNonConst(0);
280 for (
int i = 0; i < numSend; i++)
281 partsISendData[i] = 1;
283 (numPartsIRecv->getDataNonConst(0))[0] = 0;
285 numPartsIRecv->doExport(*partsISend, *partsExport, Xpetra::ADD);
286 numRecv = (numPartsIRecv->getData(0))[0];
295 for (
typename map_type::iterator it = sendMap.begin(); it != sendMap.end(); it++)
296 MPI_Isend(static_cast<void*>(it->second.getRawPtr()), it->second.size(), MpiType, Teuchos::as<GO>(it->first), msgTag, *rawMpiComm, &sendReqs[cnt++]);
299 size_t totalGIDs = myGIDs.
size();
300 for (
int i = 0; i < numRecv; i++) {
302 MPI_Probe(MPI_ANY_SOURCE, msgTag, *rawMpiComm, &status);
305 int fromRank = status.MPI_SOURCE, count;
306 MPI_Get_count(&status, MpiType, &count);
308 recvMap[fromRank].resize(count);
309 MPI_Recv(static_cast<void*>(recvMap[fromRank].getRawPtr()), count, MpiType, fromRank, msgTag, *rawMpiComm, &status);
322 for (
typename map_type::const_iterator it = recvMap.begin(); it != recvMap.end(); it++) {
323 int offset = myGIDs.
size(), len = it->second.size();
325 myGIDs.
resize(offset + len);
326 memcpy(myGIDs.
getRawPtr() + offset, it->second.getRawPtr(), len*
sizeof(GO));
331 std::sort(myGIDs.
begin(), myGIDs.
end());
334 RCP<Map> newRowMap = MapFactory ::Build(lib, rowMap->getGlobalNumElements(), myGIDs(), indexBase, origComm);
343 rowMapImporter = ImportFactory::Build(rowMap, newRowMap);
345 rowMapImporter = ImportFactory::Build(blockedRowMap->getMap(), newRowMap);
356 size_t numBlocks = blockedRowMap->getNumMaps();
357 std::vector<RCP<const Import> > subImports(numBlocks);
359 for(
size_t i=0; i<numBlocks; i++) {
362 subImports[i] = ImportFactory::Build(source,target);
364 Set(currentLevel,
"SubImporters",subImports);
368 Set(currentLevel,
"Importer", rowMapImporter);
373 if (pL.
get<
bool>(
"repartition: print partition distribution") && IsPrint(
Statistics2)) {
375 GetOStream(
Statistics2) <<
"Partition distribution over cores (ownership is indicated by '+')" << std::endl;
377 char amActive = (myGIDs.
size() ? 1 : 0);
378 std::vector<char> areActive(numProcs, 0);
379 MPI_Gather(&amActive, 1, MPI_CHAR, &areActive[0], 1, MPI_CHAR, 0, *rawMpiComm);
381 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
382 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
383 for (
int j = 0; j < rowWidth; j++)
384 if (proc + j < numProcs)
385 GetOStream(
Statistics2) << (areActive[proc + j] ?
"+" :
".");
389 GetOStream(
Statistics2) <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
396 template<
typename T,
typename W>
401 template<
typename T,
typename W>
406 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
412 int numProcs = comm->getSize();
424 const int maxLocal = pL.
get<
int>(
"repartition: remap num values");
425 const int dataSize = 2*maxLocal;
428 if (decomposition.getLocalLength() > 0)
429 decompEntries = decomposition.getDataNonConst(0);
441 std::map<GO,GO> lEdges;
442 for (LO i = 0; i < decompEntries.
size(); i++)
443 lEdges[decompEntries[i]] += A.getNumEntriesInLocalRow(i);
447 std::multimap<GO,GO> revlEdges;
448 for (
typename std::map<GO,GO>::const_iterator it = lEdges.begin(); it != lEdges.end(); it++)
449 revlEdges.insert(std::make_pair(it->second, it->first));
454 Array<GO> lData(dataSize, -1), gData(numProcs * dataSize);
456 for (
typename std::multimap<GO,GO>::reverse_iterator rit = revlEdges.rbegin(); rit != revlEdges.rend() && numEdges < maxLocal; rit++) {
457 lData[2*numEdges+0] = rit->second;
458 lData[2*numEdges+1] = rit->first;
465 MPI_Allgather(static_cast<void*>(lData.getRawPtr()), dataSize, MpiType, static_cast<void*>(gData.
getRawPtr()), dataSize, MpiType, *rawMpiComm);
470 std::vector<Triplet<int,int> > gEdges(numProcs * maxLocal);
472 for (LO i = 0; i < gData.
size(); i += 2) {
473 GO part = gData[i+0];
474 GO weight = gData[i+1];
476 gEdges[k].i = i/dataSize;
478 gEdges[k].v = weight;
486 std::sort(gEdges.begin(), gEdges.end(), compareTriplets<int,int>);
489 std::map<int,int> match;
490 std::vector<char> matchedRanks(numProcs, 0), matchedParts(numProcs, 0);
492 for (
typename std::vector<
Triplet<int,int> >::const_iterator it = gEdges.begin(); it != gEdges.end(); it++) {
495 if (matchedRanks[rank] == 0 && matchedParts[part] == 0) {
496 matchedRanks[rank] = 1;
497 matchedParts[part] = 1;
502 GetOStream(
Statistics1) <<
"Number of unassigned partitions before cleanup stage: " << (numPartitions - numMatched) <<
" / " << numPartitions << std::endl;
509 if (numPartitions - numMatched > 0) {
510 for (
int part = 0, matcher = 0; part < numProcs; part++) {
511 if (match.count(part) == 0) {
513 while (matchedRanks[matcher])
516 match[part] = matcher++;
522 for (LO i = 0; i < decompEntries.
size(); i++)
523 decompEntries[i] = match[decompEntries[i]];
528 #endif //ifdef HAVE_MPI 530 #endif // MUELU_REPARTITIONFACTORY_DEF_HPP RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Important warning messages (one line)
void reserve(size_type n)
#define MueLu_sumAll(rcpComm, in, out)
static MPI_Datatype getType()
#define MueLu_maxAll(rcpComm, in, out)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static MPI_Datatype getType()
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
Namespace for MueLu classes and methods.
void DeterminePartitionPlacement(const Matrix &A, GOVector &decomposition, GO numPartitions) const
Determine which process should own each partition.
#define SET_VALID_ENTRY(name)
Print even more statistics.
static MPI_Datatype getType()
void Build(Level ¤tLevel) const
Build an object with this factory.
static bool compareTriplets(const Triplet< T, W > &a, const Triplet< T, W > &b)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static MPI_Datatype getType()
void DeclareInput(Level ¤tLevel) const
Determines the data that RepartitionFactory needs, and the factories that generate that data...
void resize(size_type new_size, const value_type &x=value_type())
static MPI_Datatype getType()
static MPI_Datatype getType()
void push_back(const value_type &x)
Exception throws to report errors in the internal logical of the program.
std::string toString(const T &t)