46 #ifndef MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
47 #define MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
49 #include <Teuchos_Tuple.hpp>
67 #include "MueLu_HierarchyUtils.hpp"
70 #include "MueLu_PerfUtils.hpp"
74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 validParamList->
set<
RCP<const FactoryBase> >(
"P", Teuchos::null,
"Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
79 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Factory for generating the non-rebalanced coarse level A. We need this to make sure the non-rebalanced coarse A is calculated first before rebalancing takes place.");
81 validParamList->
set<
RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Factory for generating the non-rebalanced Coordinates.");
82 validParamList->
set<
RCP<const FactoryBase> >(
"Importer", Teuchos::null,
"Generating factory of the matrix Importer for rebalancing");
83 validParamList->
set<
RCP<const FactoryBase> >(
"SubImporters", Teuchos::null,
"Generating factory of the matrix sub-Importers for rebalancing");
85 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
87 #undef SET_VALID_ENTRY
93 return validParamList;
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 FactManager_.push_back(FactManager);
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
103 Input(coarseLevel,
"P");
104 Input(coarseLevel,
"A");
106 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
107 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
113 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
114 if((*it)->hasFactory(
"Coordinates") ==
true)
115 coarseLevel.
DeclareInput(
"Coordinates",(*it)->GetFactory(
"Coordinates").get(),
this);
119 if(FactManager_.size() == 0) {
120 Input(coarseLevel,
"Importer");
121 Input(coarseLevel,
"SubImporters");
122 Input(coarseLevel,
"Coordinates");
128 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
134 bool UseSingleSource = FactManager_.size() == 0;
140 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"P");
156 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
157 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
160 std::vector<GO> fullRangeMapVector;
161 std::vector<GO> fullDomainMapVector;
162 std::vector<RCP<const Map> > subBlockPRangeMaps;
163 std::vector<RCP<const Map> > subBlockPDomainMaps;
164 subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows());
165 subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols());
167 std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
168 subBlockRebP.reserve(bOriginalTransferOp->Rows());
171 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
172 std::vector<RCP<xdMV> > newCoordinates(bOriginalTransferOp->Rows());
174 if(UseSingleSource) {
175 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
176 oldCoordinates = Get<RCP<xdBV> >(coarseLevel,
"Coordinates");
181 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
182 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
188 if(UseSingleSource) rebalanceImporter = importers[curBlockId];
194 TEUCHOS_TEST_FOR_EXCEPTION(Pwii == Teuchos::null,
Xpetra::Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: block " << curBlockId <<
" is not of type CrsMatrixWrap. We need an underlying CsrMatrix to replace domain map and importer!");
198 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Pii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
201 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Pii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
204 if(rebalanceImporter != Teuchos::null) {
205 std::stringstream ss; ss <<
"Rebalancing prolongator block P(" << curBlockId <<
"," << curBlockId <<
")";
221 SubFactoryMonitor subM(*
this,
"Rebalancing prolongator -- fast map replacement", coarseLevel);
222 newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
223 Pwii->getCrsMatrix()->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
227 params->
set(
"printLoadBalancingInfo",
true);
228 std::stringstream ss2; ss2 <<
"P(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
232 subBlockRebP.push_back(Pii);
237 if(UseSingleSource) {
238 RCP<xdMV> localCoords = oldCoordinates->getMultiVector(curBlockId);
244 permutedLocalCoords->doImport(*localCoords, *coordImporter,
Xpetra::INSERT);
246 newCoordinates[curBlockId] = permutedLocalCoords;
248 else if ( (*it)->hasFactory(
"Coordinates") ==
true && coarseLevel.IsAvailable(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()) ==
true) {
249 RCP<xdMV> coords = coarseLevel.Get<
RCP<xdMV> >(
"Coordinates",(*it)->GetFactory(
"Coordinates").get());
254 LO nodeNumElts = coords->getMap()->getNodeNumElements();
257 LO myBlkSize = 0, blkSize = 0;
259 if (nodeNumElts > 0) {
262 "MueLu::RebalanceBlockInterpolationFactory::Build: block size. " << rebalanceImporter->getSourceMap()->getNodeNumElements() <<
" not divisable by " << nodeNumElts);
263 myBlkSize = rebalanceImporter->getSourceMap()->getNodeNumElements() / nodeNumElts;
266 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
270 coordImporter = rebalanceImporter;
276 GO indexBase = origMap->getIndexBase();
279 LO numEntries = OEntries.
size()/blkSize;
281 for (LO i = 0; i < numEntries; i++)
282 Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
284 RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
285 coordImporter = ImportFactory::Build(origMap, targetMap);
289 permutedCoords->doImport(*coords, *coordImporter,
Xpetra::INSERT);
292 if (pL.
isParameter(
"repartition: use subcommunicators") ==
true && pL.
get<
bool>(
"repartition: use subcommunicators") ==
true)
293 permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
295 Set(coarseLevel,
"Coordinates", permutedCoords);
300 params->
set(
"printLoadBalancingInfo",
true);
301 std::stringstream ss; ss <<
"P(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
304 subBlockRebP.push_back(Pii);
309 if((*it)->hasFactory(
"Coordinates") ==
true && coarseLevel.IsAvailable(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()) ==
true) {
310 coarseLevel.Set(
"Coordinates", coarseLevel.Get<
RCP<xdMV> >(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()),
this);
316 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
318 if(orig_stridedRgMap != Teuchos::null) {
319 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
321 stridedRgMap = StridedMapFactory::Build(
322 originalTransferOp->getRangeMap()->lib(),
325 Pii->getRangeMap()->getIndexBase(),
327 originalTransferOp->getRangeMap()->getComm(),
328 orig_stridedRgMap->getStridedBlockId(),
329 orig_stridedRgMap->getOffset());
330 }
else stridedRgMap = Pii->getRangeMap();
333 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
336 if(orig_stridedDoMap != Teuchos::null) {
337 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
339 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
342 Pii->getDomainMap()->getIndexBase(),
344 originalTransferOp->getDomainMap()->getComm(),
345 orig_stridedDoMap->getStridedBlockId(),
346 orig_stridedDoMap->getOffset());
347 }
else stridedDoMap = Pii->getDomainMap();
349 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
350 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
353 if(Pii->IsView(
"stridedMaps")) Pii->RemoveView(
"stridedMaps");
354 Pii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
358 subBlockPRangeMaps.push_back(rangeMapii);
361 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.
begin(), nodeRangeMapii.
end());
362 if(bThyraRangeGIDs ==
false)
363 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
367 subBlockPDomainMaps.push_back(domainMapii);
370 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.
begin(), nodeDomainMapii.
end());
371 if(bThyraDomainGIDs ==
false)
372 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
379 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
380 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
383 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
386 if(stridedRgFullMap != Teuchos::null) {
387 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
389 StridedMapFactory::Build(
390 originalTransferOp->getRangeMap()->lib(),
395 originalTransferOp->getRangeMap()->getComm(),
396 stridedRgFullMap->getStridedBlockId(),
397 stridedRgFullMap->getOffset());
401 originalTransferOp->getRangeMap()->lib(),
405 originalTransferOp->getRangeMap()->getComm());
409 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
412 if(stridedDoFullMap != Teuchos::null) {
413 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
414 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
416 StridedMapFactory::Build(
417 originalTransferOp->getDomainMap()->lib(),
422 originalTransferOp->getDomainMap()->getComm(),
423 stridedDoFullMap->getStridedBlockId(),
424 stridedDoFullMap->getOffset());
429 originalTransferOp->getDomainMap()->lib(),
433 originalTransferOp->getDomainMap()->getComm());
438 MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
440 MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
443 for(
size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
445 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null,
Xpetra::Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: block P" << i <<
" is not of type CrsMatrixWrap.");
446 bRebP->setMatrix(i,i,crsOpii);
448 bRebP->fillComplete();
450 Set(coarseLevel,
"P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
454 if(UseSingleSource) {
455 RCP<xdBV> bcoarseCoords =
rcp(
new xdBV(rebrangeMapExtractor->getBlockedMap(),newCoordinates));
456 Set(coarseLevel,
"Coordinates",bcoarseCoords);
Exception indicating invalid cast attempted.
#define MueLu_maxAll(rcpComm, in, out)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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)
bool isParameter(const std::string &name) const
Print statistics that do not involve significant additional computation.
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 std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Exception throws to report errors in the internal logical of the program.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
An exception safe way to call the method 'Level::SetFactoryManager()'.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()