46 #ifndef MUELU_REBALANCETRANSFERFACTORY_DEF_HPP 47 #define MUELU_REBALANCETRANSFERFACTORY_DEF_HPP 50 #include <Teuchos_Tuple.hpp> 52 #include "Xpetra_MultiVector.hpp" 53 #include "Xpetra_MultiVectorFactory.hpp" 54 #include "Xpetra_Vector.hpp" 55 #include "Xpetra_VectorFactory.hpp" 56 #include <Xpetra_Matrix.hpp> 57 #include <Xpetra_MapFactory.hpp> 58 #include <Xpetra_MatrixFactory.hpp> 59 #include <Xpetra_Import.hpp> 60 #include <Xpetra_ImportFactory.hpp> 61 #include <Xpetra_IO.hpp> 68 #include "MueLu_PerfUtils.hpp" 72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 RCP<ParameterList> validParamList = rcp(
new ParameterList());
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 81 #undef SET_VALID_ENTRY 84 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
85 RCP<validatorType> typeValidator = rcp (
new validatorType(Teuchos::tuple<std::string>(
"Interpolation",
"Restriction"),
"type"));
86 validParamList->set(
"type",
"Interpolation",
"Type of the transfer operator that need to be rebalanced (Interpolation or Restriction)", typeValidator);
89 validParamList->set< RCP<const FactoryBase> >(
"P", null,
"Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
90 validParamList->set< RCP<const FactoryBase> >(
"R", null,
"Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
91 validParamList->set< RCP<const FactoryBase> >(
"Nullspace", null,
"Factory of the nullspace that need to be rebalanced (only used if type=Interpolation)");
92 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", null,
"Factory of the coordinates that need to be rebalanced (only used if type=Interpolation)");
93 validParamList->set< RCP<const FactoryBase> >(
"BlockNumber", null,
"Factory of the block ids that need to be rebalanced (only used if type=Interpolation)");
94 validParamList->set< RCP<const FactoryBase> >(
"Importer", null,
"Factory of the importer object used for the rebalancing");
95 validParamList->set<
int > (
"write start", -1,
"First level at which coordinates should be written to file");
96 validParamList->set<
int > (
"write end", -1,
"Last level at which coordinates should be written to file");
102 return validParamList;
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 const ParameterList& pL = GetParameterList();
109 if (pL.get<std::string>(
"type") ==
"Interpolation") {
110 Input(coarseLevel,
"P");
111 if (pL.get<
bool>(
"repartition: rebalance Nullspace"))
112 Input(coarseLevel,
"Nullspace");
113 if (pL.get< RCP<const FactoryBase> >(
"Coordinates") != Teuchos::null)
114 Input(coarseLevel,
"Coordinates");
115 if (pL.get< RCP<const FactoryBase> >(
"BlockNumber") != Teuchos::null)
116 Input(coarseLevel,
"BlockNumber");
119 if (pL.get<
bool>(
"transpose: use implicit") ==
false)
120 Input(coarseLevel,
"R");
123 Input(coarseLevel,
"Importer");
126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
131 const ParameterList& pL = GetParameterList();
133 int implicit = !pL.get<
bool>(
"repartition: rebalance P and R");
134 int writeStart = pL.get<
int> (
"write start");
135 int writeEnd = pL.get<
int> (
"write end");
137 if (writeStart == 0 && fineLevel.
GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel,
"Coordinates")) {
138 std::string fileName =
"coordinates_level_0.m";
139 RCP<xdMV> fineCoords = fineLevel.
Get< RCP<xdMV> >(
"Coordinates");
140 if (fineCoords != Teuchos::null)
141 Xpetra::IO<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *fineCoords);
144 if (writeStart == 0 && fineLevel.
GetLevelID() == 0 && writeStart <= writeEnd && IsAvailable(fineLevel,
"BlockNumber")) {
145 std::string fileName =
"BlockNumber_level_0.m";
146 RCP<LocalOrdinalVector> fineBlockNumber = fineLevel.
Get< RCP<LocalOrdinalVector> >(
"BlockNumber");
147 if (fineBlockNumber != Teuchos::null)
148 Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *fineBlockNumber);
151 RCP<const Import> importer = Get<RCP<const Import> >(coarseLevel,
"Importer");
157 RCP<ParameterList> params = rcp(
new ParameterList());
159 params->set(
"printLoadBalancingInfo",
true);
160 params->set(
"printCommInfo",
true);
163 std::string transferType = pL.get<std::string>(
"type");
164 if (transferType ==
"Interpolation") {
165 RCP<Matrix> originalP = Get< RCP<Matrix> >(coarseLevel,
"P");
171 if (implicit || importer.is_null()) {
172 GetOStream(
Runtime0) <<
"Using original prolongator" << std::endl;
173 Set(coarseLevel,
"P", originalP);
191 RCP<Matrix> rebalancedP = originalP;
192 RCP<const CrsMatrixWrap> crsOp = rcp_dynamic_cast<
const CrsMatrixWrap>(originalP);
193 TEUCHOS_TEST_FOR_EXCEPTION(crsOp == Teuchos::null,
Exceptions::BadCast,
"Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
195 RCP<CrsMatrix> rebalancedP2 = crsOp->getCrsMatrix();
196 TEUCHOS_TEST_FOR_EXCEPTION(rebalancedP2 == Teuchos::null, std::runtime_error,
"Xpetra::CrsMatrixWrap doesn't have a CrsMatrix");
199 SubFactoryMonitor subM(*
this,
"Rebalancing prolongator -- fast map replacement", coarseLevel);
201 RCP<const Import> newImporter;
204 newImporter = ImportFactory::Build(importer->getTargetMap(), rebalancedP->getColMap());
206 rebalancedP2->replaceDomainMapAndImporter(importer->getTargetMap(), newImporter);
215 if(!rebalancedP.is_null()) {std::ostringstream oss; oss <<
"P_" << coarseLevel.GetLevelID(); rebalancedP->setObjectLabel(oss.str());}
216 Set(coarseLevel,
"P", rebalancedP);
223 if (importer.is_null()) {
224 if (IsAvailable(coarseLevel,
"Nullspace"))
225 Set(coarseLevel,
"Nullspace", Get<RCP<MultiVector> >(coarseLevel,
"Nullspace"));
227 if (pL.isParameter(
"Coordinates") && pL.get< RCP<const FactoryBase> >(
"Coordinates") != Teuchos::null)
228 if (IsAvailable(coarseLevel,
"Coordinates"))
229 Set(coarseLevel,
"Coordinates", Get< RCP<xdMV> >(coarseLevel,
"Coordinates"));
231 if (pL.isParameter(
"BlockNumber") && pL.get< RCP<const FactoryBase> >(
"BlockNumber") != Teuchos::null)
232 if (IsAvailable(coarseLevel,
"BlockNumber"))
233 Set(coarseLevel,
"BlockNumber", Get< RCP<LocalOrdinalVector> >(coarseLevel,
"BlockNumber"));
238 if (pL.isParameter(
"Coordinates") &&
239 pL.get< RCP<const FactoryBase> >(
"Coordinates") != Teuchos::null &&
240 IsAvailable(coarseLevel,
"Coordinates")) {
241 RCP<xdMV> coords = Get<RCP<xdMV> >(coarseLevel,
"Coordinates");
246 LO nodeNumElts = coords->getMap()->getNodeNumElements();
249 LO myBlkSize = 0, blkSize = 0;
251 myBlkSize = importer->getSourceMap()->getNodeNumElements() / nodeNumElts;
252 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
254 RCP<const Import> coordImporter;
256 coordImporter = importer;
262 RCP<const Map> origMap = coords->getMap();
263 GO indexBase = origMap->getIndexBase();
265 ArrayView<const GO> OEntries = importer->getTargetMap()->getNodeElementList();
266 LO numEntries = OEntries.size()/blkSize;
267 ArrayRCP<GO> Entries(numEntries);
268 for (LO i = 0; i < numEntries; i++)
269 Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
271 RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
272 coordImporter = ImportFactory::Build(origMap, targetMap);
275 RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
276 permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
278 if (pL.isParameter(
"repartition: use subcommunicators") ==
true && pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
279 permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
281 if (permutedCoords->getMap() == Teuchos::null)
282 permutedCoords = Teuchos::null;
284 Set(coarseLevel,
"Coordinates", permutedCoords);
286 std::string fileName =
"rebalanced_coordinates_level_" +
toString(coarseLevel.GetLevelID()) +
".m";
287 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedCoords->getMap() != Teuchos::null)
288 Xpetra::IO<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Write(fileName, *permutedCoords);
291 if (pL.isParameter(
"BlockNumber") &&
292 pL.get< RCP<const FactoryBase> >(
"BlockNumber") != Teuchos::null &&
293 IsAvailable(coarseLevel,
"BlockNumber")) {
294 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(coarseLevel,
"BlockNumber");
299 RCP<LocalOrdinalVector> permutedBlockNumber = LocalOrdinalVectorFactory::Build(importer->getTargetMap(),
false);
300 permutedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
302 if (pL.isParameter(
"repartition: use subcommunicators") ==
true && pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
303 permutedBlockNumber->replaceMap(permutedBlockNumber->getMap()->removeEmptyProcesses());
305 if (permutedBlockNumber->getMap() == Teuchos::null)
306 permutedBlockNumber = Teuchos::null;
308 Set(coarseLevel,
"BlockNumber", permutedBlockNumber);
310 std::string fileName =
"rebalanced_BlockNumber_level_" +
toString(coarseLevel.GetLevelID()) +
".m";
311 if (writeStart <= coarseLevel.GetLevelID() && coarseLevel.GetLevelID() <= writeEnd && permutedBlockNumber->getMap() != Teuchos::null)
312 Xpetra::IO<LO,LO,GO,NO>::Write(fileName, *permutedBlockNumber);
315 if (IsAvailable(coarseLevel,
"Nullspace")) {
316 RCP<MultiVector> nullspace = Get< RCP<MultiVector> >(coarseLevel,
"Nullspace");
321 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(importer->getTargetMap(), nullspace->getNumVectors());
322 permutedNullspace->doImport(*nullspace, *importer, Xpetra::INSERT);
324 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
325 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
327 if (permutedNullspace->getMap() == Teuchos::null)
328 permutedNullspace = Teuchos::null;
330 Set(coarseLevel,
"Nullspace", permutedNullspace);
334 if (pL.get<
bool>(
"transpose: use implicit") ==
false) {
335 RCP<Matrix> originalR = Get< RCP<Matrix> >(coarseLevel,
"R");
339 if (implicit || importer.is_null()) {
340 GetOStream(
Runtime0) <<
"Using original restrictor" << std::endl;
341 Set(coarseLevel,
"R", originalR);
344 RCP<Matrix> rebalancedR;
346 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
349 Teuchos::ParameterList listLabel;
350 listLabel.set(
"Timer Label",
"MueLu::RebalanceR-" +
Teuchos::toString(coarseLevel.GetLevelID()));
351 rebalancedR = MatrixFactory::Build(originalR, *importer, dummy, importer->getTargetMap(),Teuchos::rcp(&listLabel,
false));
353 if(!rebalancedR.is_null()) {std::ostringstream oss; oss <<
"R_" << coarseLevel.GetLevelID(); rebalancedR->setObjectLabel(oss.str());}
354 Set(coarseLevel,
"R", rebalancedR);
372 #endif // MUELU_REBALANCETRANSFERFACTORY_DEF_HPP Exception indicating invalid cast attempted.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
#define MueLu_maxAll(rcpComm, in, out)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
static const NoFactory * get()
Print even more statistics.
int GetLevelID() const
Return level number.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
#define SET_VALID_ENTRY(name)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)