46 #ifndef MUELU_REPARTITIONFACTORY_DEF_HPP 47 #define MUELU_REPARTITIONFACTORY_DEF_HPP 56 #include <Teuchos_DefaultMpiComm.hpp> 57 #include <Teuchos_CommHelpers.hpp> 58 #include <Teuchos_Details_MpiTypeTraits.hpp> 60 #include <Xpetra_Map.hpp> 61 #include <Xpetra_MapFactory.hpp> 62 #include <Xpetra_MultiVectorFactory.hpp> 63 #include <Xpetra_VectorFactory.hpp> 64 #include <Xpetra_Import.hpp> 65 #include <Xpetra_ImportFactory.hpp> 66 #include <Xpetra_Export.hpp> 67 #include <Xpetra_ExportFactory.hpp> 68 #include <Xpetra_Matrix.hpp> 69 #include <Xpetra_MatrixFactory.hpp> 71 #include "MueLu_Utilities.hpp" 73 #include "MueLu_CloneRepartitionInterface.hpp" 78 #include "MueLu_PerfUtils.hpp" 82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 RCP<ParameterList> validParamList = rcp(
new ParameterList());
86 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 92 #undef SET_VALID_ENTRY 94 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Factory of the matrix A");
95 validParamList->set< RCP<const FactoryBase> >(
"number of partitions", Teuchos::null,
"Instance of RepartitionHeuristicFactory.");
96 validParamList->set< RCP<const FactoryBase> >(
"Partition", Teuchos::null,
"Factory of the partition");
98 return validParamList;
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 Input(currentLevel,
"A");
104 Input(currentLevel,
"number of partitions");
105 Input(currentLevel,
"Partition");
108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 const Teuchos::ParameterList & pL = GetParameterList();
115 bool remapPartitions = pL.get<
bool> (
"repartition: remap parts");
118 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
119 RCP<const Map> rowMap = A->getRowMap();
120 GO indexBase = rowMap->getIndexBase();
121 Xpetra::UnderlyingLib lib = rowMap->lib();
123 RCP<const Teuchos::Comm<int> > origComm = rowMap->getComm();
124 RCP<const Teuchos::Comm<int> > comm = origComm;
126 int myRank = comm->getRank();
127 int numProcs = comm->getSize();
129 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
130 TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null,
Exceptions::RuntimeError,
"Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
131 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
134 int numPartitions = Get<int>(currentLevel,
"number of partitions");
139 RCP<GOVector> decomposition = Get<RCP<GOVector> >(currentLevel,
"Partition");
142 if(remapPartitions ==
true && Teuchos::rcp_dynamic_cast<const CloneRepartitionInterface>(GetFactory(
"Partition")) != Teuchos::null) {
146 remapPartitions =
false;
150 if (numPartitions == 1) {
155 GetOStream(
Warnings0) <<
"Only one partition: Skip call to the repartitioner." << std::endl;
156 }
else if (numPartitions == -1) {
158 GetOStream(
Warnings0) <<
"No repartitioning necessary: partitions were left unchanged by the repartitioner" << std::endl;
159 Set<RCP<const Import> >(currentLevel,
"Importer", Teuchos::null);
164 const int nodeRepartLevel = pL.get<
int> (
"repartition: node repartition level");
165 if(currentLevel.GetLevelID() == nodeRepartLevel) {
170 remapPartitions =
false;
210 if (remapPartitions) {
213 bool acceptPartition = pL.get<
bool>(
"repartition: remap accept partition");
214 bool allSubdomainsAcceptPartitions;
215 int localNumAcceptPartition = acceptPartition;
216 int globalNumAcceptPartition;
217 MueLu_sumAll(comm, localNumAcceptPartition, globalNumAcceptPartition);
218 GetOStream(
Statistics2) <<
"Number of ranks that accept partitions: " << globalNumAcceptPartition << std::endl;
219 if (globalNumAcceptPartition < numPartitions) {
220 GetOStream(
Warnings0) <<
"Not enough ranks are willing to accept a partition, allowing partitions on all ranks." << std::endl;
221 acceptPartition =
true;
222 allSubdomainsAcceptPartitions =
true;
223 }
else if (numPartitions > numProcs) {
225 allSubdomainsAcceptPartitions =
true;
227 allSubdomainsAcceptPartitions =
false;
230 DeterminePartitionPlacement(*A, *decomposition, numPartitions, acceptPartition, allSubdomainsAcceptPartitions);
241 ArrayRCP<const GO> decompEntries;
242 if (decomposition->getLocalLength() > 0)
243 decompEntries = decomposition->getData(0);
245 #ifdef HAVE_MUELU_DEBUG 247 int incorrectRank = -1;
248 for (
int i = 0; i < decompEntries.size(); i++)
249 if (decompEntries[i] >= numProcs || decompEntries[i] < 0) {
250 incorrectRank = myRank;
254 int incorrectGlobalRank = -1;
260 myGIDs.reserve(decomposition->getLocalLength());
265 typedef std::map<GO, Array<GO> > map_type;
267 for (LO i = 0; i < decompEntries.size(); i++) {
268 GO
id = decompEntries[i];
269 GO GID = rowMap->getGlobalElement(i);
272 myGIDs .push_back(GID);
274 sendMap[id].push_back(GID);
276 decompEntries = Teuchos::null;
279 GO numLocalKept = myGIDs.size(), numGlobalKept, numGlobalRows = A->getGlobalNumRows();
281 GetOStream(
Statistics2) <<
"Unmoved rows: " << numGlobalKept <<
" / " << numGlobalRows <<
" (" << 100*Teuchos::as<double>(numGlobalKept)/numGlobalRows <<
"%)" << std::endl;
284 int numSend = sendMap.size(), numRecv;
287 Array<GO> myParts(numSend), myPart(1);
290 for (
typename map_type::const_iterator it = sendMap.begin(); it != sendMap.end(); it++)
291 myParts[cnt++] = it->first;
295 GO partsIndexBase = 0;
296 RCP<Map> partsIHave = MapFactory ::Build(lib, Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), myParts(), partsIndexBase, comm);
297 RCP<Map> partsIOwn = MapFactory ::Build(lib, numProcs, myPart(), partsIndexBase, comm);
298 RCP<Export> partsExport = ExportFactory::Build(partsIHave, partsIOwn);
300 RCP<GOVector> partsISend = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIHave);
301 RCP<GOVector> numPartsIRecv = Xpetra::VectorFactory<GO, LO, GO, NO>::Build(partsIOwn);
303 ArrayRCP<GO> partsISendData = partsISend->getDataNonConst(0);
304 for (
int i = 0; i < numSend; i++)
305 partsISendData[i] = 1;
307 (numPartsIRecv->getDataNonConst(0))[0] = 0;
309 numPartsIRecv->doExport(*partsISend, *partsExport, Xpetra::ADD);
310 numRecv = (numPartsIRecv->getData(0))[0];
313 MPI_Datatype MpiType = Teuchos::Details::MpiTypeTraits<GO>::getType();
317 Array<MPI_Request> sendReqs(numSend);
319 for (
typename map_type::iterator it = sendMap.begin(); it != sendMap.end(); it++)
320 MPI_Isend(static_cast<void*>(it->second.getRawPtr()), it->second.size(), MpiType, Teuchos::as<GO>(it->first), msgTag, *rawMpiComm, &sendReqs[cnt++]);
323 size_t totalGIDs = myGIDs.size();
324 for (
int i = 0; i < numRecv; i++) {
326 MPI_Probe(MPI_ANY_SOURCE, msgTag, *rawMpiComm, &status);
329 int fromRank = status.MPI_SOURCE, count;
330 MPI_Get_count(&status, MpiType, &count);
332 recvMap[fromRank].resize(count);
333 MPI_Recv(static_cast<void*>(recvMap[fromRank].getRawPtr()), count, MpiType, fromRank, msgTag, *rawMpiComm, &status);
340 Array<MPI_Status> sendStatuses(numSend);
341 MPI_Waitall(numSend, sendReqs.getRawPtr(), sendStatuses.getRawPtr());
345 myGIDs.reserve(totalGIDs);
346 for (
typename map_type::const_iterator it = recvMap.begin(); it != recvMap.end(); it++) {
347 int offset = myGIDs.size(), len = it->second.size();
349 myGIDs.resize(offset + len);
350 memcpy(myGIDs.getRawPtr() + offset, it->second.getRawPtr(), len*
sizeof(GO));
355 std::sort(myGIDs.begin(), myGIDs.end());
358 RCP<Map> newRowMap = MapFactory ::Build(lib, rowMap->getGlobalNumElements(), myGIDs(), indexBase, origComm);
359 RCP<const Import> rowMapImporter;
361 RCP<const BlockedMap> blockedRowMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
366 if(blockedRowMap.is_null())
367 rowMapImporter = ImportFactory::Build(rowMap, newRowMap);
369 rowMapImporter = ImportFactory::Build(blockedRowMap->getMap(), newRowMap);
374 if(!blockedRowMap.is_null()) {
380 size_t numBlocks = blockedRowMap->getNumMaps();
381 std::vector<RCP<const Import> > subImports(numBlocks);
383 for(
size_t i=0; i<numBlocks; i++) {
384 RCP<const Map> source = blockedRowMap->getMap(i);
385 RCP<const Map> target = blockedTargetMap->getMap(i);
386 subImports[i] = ImportFactory::Build(source,target);
388 Set(currentLevel,
"SubImporters",subImports);
392 Set(currentLevel,
"Importer", rowMapImporter);
397 if (!rowMapImporter.is_null() && IsPrint(
Statistics2)) {
402 if (pL.get<
bool>(
"repartition: print partition distribution") && IsPrint(
Statistics2)) {
404 GetOStream(
Statistics2) <<
"Partition distribution over cores (ownership is indicated by '+')" << std::endl;
406 char amActive = (myGIDs.size() ? 1 : 0);
407 std::vector<char> areActive(numProcs, 0);
408 MPI_Gather(&amActive, 1, MPI_CHAR, &areActive[0], 1, MPI_CHAR, 0, *rawMpiComm);
410 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
411 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
412 for (
int j = 0; j < rowWidth; j++)
413 if (proc + j < numProcs)
414 GetOStream(
Statistics2) << (areActive[proc + j] ?
"+" :
".");
418 GetOStream(
Statistics2) <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
425 template<
typename T,
typename W>
430 template<
typename T,
typename W>
435 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
437 DeterminePartitionPlacement(
const Matrix& A, GOVector& decomposition, GO numPartitions,
bool willAcceptPartition,
bool allSubdomainsAcceptPartitions)
const {
438 RCP<const Map> rowMap = A.getRowMap();
440 RCP<const Teuchos::Comm<int> > comm = rowMap->getComm()->duplicate();
441 int numProcs = comm->getSize();
443 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
444 TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null,
Exceptions::RuntimeError,
"Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
445 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
447 const Teuchos::ParameterList& pL = GetParameterList();
453 const int maxLocal = pL.get<
int>(
"repartition: remap num values");
454 const int dataSize = 2*maxLocal;
456 ArrayRCP<GO> decompEntries;
457 if (decomposition.getLocalLength() > 0)
458 decompEntries = decomposition.getDataNonConst(0);
470 std::map<GO,GO> lEdges;
471 if (willAcceptPartition)
472 for (LO i = 0; i < decompEntries.size(); i++)
473 lEdges[decompEntries[i]] += A.getNumEntriesInLocalRow(i);
477 std::multimap<GO,GO> revlEdges;
478 for (
typename std::map<GO,GO>::const_iterator it = lEdges.begin(); it != lEdges.end(); it++)
479 revlEdges.insert(std::make_pair(it->second, it->first));
484 Array<GO> lData(dataSize, -1), gData(numProcs * dataSize);
486 for (
typename std::multimap<GO,GO>::reverse_iterator rit = revlEdges.rbegin(); rit != revlEdges.rend() && numEdges < maxLocal; rit++) {
487 lData[2*numEdges+0] = rit->second;
488 lData[2*numEdges+1] = rit->first;
494 MPI_Datatype MpiType = Teuchos::Details::MpiTypeTraits<GO>::getType();
495 MPI_Allgather(static_cast<void*>(lData.getRawPtr()), dataSize, MpiType, static_cast<void*>(gData.getRawPtr()), dataSize, MpiType, *rawMpiComm);
500 Teuchos::Array<Triplet<int,int> > gEdges(numProcs * maxLocal);
501 Teuchos::Array<bool> procWillAcceptPartition(numProcs, allSubdomainsAcceptPartitions);
503 for (LO i = 0; i < gData.size(); i += 2) {
504 int procNo = i/dataSize;
505 GO part = gData[i+0];
506 GO weight = gData[i+1];
508 gEdges[k].i = procNo;
510 gEdges[k].v = weight;
511 procWillAcceptPartition[procNo] =
true;
519 std::sort(gEdges.begin(), gEdges.end(), compareTriplets<int,int>);
522 std::map<int,int> match;
523 Teuchos::Array<char> matchedRanks(numProcs, 0), matchedParts(numPartitions, 0);
525 for (
typename Teuchos::Array<
Triplet<int,int> >::const_iterator it = gEdges.begin(); it != gEdges.end(); it++) {
528 if (matchedRanks[rank] == 0 && matchedParts[part] == 0) {
529 matchedRanks[rank] = 1;
530 matchedParts[part] = 1;
535 GetOStream(
Statistics1) <<
"Number of unassigned partitions before cleanup stage: " << (numPartitions - numMatched) <<
" / " << numPartitions << std::endl;
542 if (numPartitions - numMatched > 0) {
543 Teuchos::Array<char> partitionCounts(numPartitions, 0);
544 for (
typename std::map<int,int>::const_iterator it = match.begin(); it != match.end(); it++)
545 partitionCounts[it->first] += 1;
546 for (
int part = 0, matcher = 0; part < numPartitions; part++) {
547 if (partitionCounts[part] == 0) {
549 while (matchedRanks[matcher] || !procWillAcceptPartition[matcher])
552 match[part] = matcher++;
558 TEUCHOS_TEST_FOR_EXCEPTION(numMatched != numPartitions,
Exceptions::RuntimeError,
"MueLu::RepartitionFactory::DeterminePartitionPlacement: Only " << numMatched <<
" partitions out of " << numPartitions <<
" got assigned to ranks.");
561 for (LO i = 0; i < decompEntries.size(); i++)
562 decompEntries[i] = match[decompEntries[i]];
567 #endif //ifdef HAVE_MPI 569 #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)
#define MueLu_sumAll(rcpComm, in, out)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
#define MueLu_maxAll(rcpComm, in, out)
Timer to be used in factories. Similar to Monitor but with additional timers.
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, bool willAcceptPartition=true, bool allSubdomainsAcceptPartitions=true) const
Determine which process should own each partition.
#define SET_VALID_ENTRY(name)
static std::string PrintImporterInfo(RCP< const Import > importer, const std::string &msgTag)
Print even more statistics.
void Build(Level ¤tLevel) const
Build an object with this factory.
static bool compareTriplets(const Triplet< T, W > &a, const Triplet< T, W > &b)
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void DeclareInput(Level ¤tLevel) const
Determines the data that RepartitionFactory needs, and the factories that generate that data...
Exception throws to report errors in the internal logical of the program.