46#ifndef MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
47#define MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_
49#include <Teuchos_Tuple.hpp>
51#include "Xpetra_MultiVector.hpp"
52#include "Xpetra_MultiVectorFactory.hpp"
53#include "Xpetra_Vector.hpp"
54#include "Xpetra_VectorFactory.hpp"
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MapFactory.hpp>
58#include <Xpetra_MapExtractor.hpp>
59#include <Xpetra_MapExtractorFactory.hpp>
60#include <Xpetra_MatrixFactory.hpp>
61#include <Xpetra_Import.hpp>
62#include <Xpetra_ImportFactory.hpp>
67#include "MueLu_HierarchyUtils.hpp"
70#include "MueLu_PerfUtils.hpp"
74template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 RCP<ParameterList> validParamList = rcp(
new ParameterList());
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))
93 return validParamList;
96template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 FactManager_.push_back(FactManager);
101template <
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");
128template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
132 typedef Xpetra::BlockedMultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> xdBV;
134 bool UseSingleSource = FactManager_.size() == 0;
137 Teuchos::RCP<Matrix> nonrebCoarseA = Get< RCP<Matrix> >(coarseLevel,
"A");
139 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
140 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"P");
142 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
143 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
144 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
146 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
147 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
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());
173 RCP<xdBV> oldCoordinates;
174 if(UseSingleSource) {
175 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
176 oldCoordinates = Get<RCP<xdBV> >(coarseLevel,
"Coordinates");
180 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
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];
189 else rebalanceImporter = coarseLevel.
Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
192 Teuchos::RCP<Matrix> Pii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
193 Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Pii);
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 <<
")";
219 RCP<const Import> newImporter;
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);
226 RCP<ParameterList> params = rcp(
new ParameterList());
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);
241 RCP<const Import> coordImporter = rebalanceImporter;
243 RCP<xdMV> permutedLocalCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), localCoords->getNumVectors());
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()->getLocalNumElements();
257 LO myBlkSize = 0, blkSize = 0;
259 if (nodeNumElts > 0) {
262 "MueLu::RebalanceBlockInterpolationFactory::Build: block size. " << rebalanceImporter->getSourceMap()->getLocalNumElements() <<
" not divisable by " << nodeNumElts);
263 myBlkSize = rebalanceImporter->getSourceMap()->getLocalNumElements() / nodeNumElts;
266 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
268 RCP<const Import> coordImporter = Teuchos::null;
270 coordImporter = rebalanceImporter;
275 RCP<const Map> origMap = coords->getMap();
276 GO indexBase = origMap->getIndexBase();
278 ArrayView<const GO> OEntries = rebalanceImporter->getTargetMap()->getLocalElementList();
279 LO numEntries = OEntries.size()/blkSize;
280 ArrayRCP<GO> Entries(numEntries);
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);
288 RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
289 permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
291 const ParameterList& pL = GetParameterList();
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);
299 RCP<ParameterList> params = rcp(
new ParameterList());
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()));
317 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
318 if(orig_stridedRgMap != Teuchos::null) {
319 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
320 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
321 stridedRgMap = StridedMapFactory::Build(
322 originalTransferOp->getRangeMap()->lib(),
323 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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()));
335 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
336 if(orig_stridedDoMap != Teuchos::null) {
337 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
338 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
339 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
340 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
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);
357 Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap(
"stridedMaps");
358 subBlockPRangeMaps.push_back(rangeMapii);
359 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
361 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
362 if(bThyraRangeGIDs ==
false)
363 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
366 Teuchos::RCP<const Map> domainMapii = Pii->getColMap(
"stridedMaps");
367 subBlockPDomainMaps.push_back(domainMapii);
368 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
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();
382 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
383 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
384 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangePMapExtractor->getFullMap());
385 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
386 if(stridedRgFullMap != Teuchos::null) {
387 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
389 StridedMapFactory::Build(
390 originalTransferOp->getRangeMap()->lib(),
391 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
395 originalTransferOp->getRangeMap()->getComm(),
396 stridedRgFullMap->getStridedBlockId(),
397 stridedRgFullMap->getOffset());
401 originalTransferOp->getRangeMap()->lib(),
402 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
405 originalTransferOp->getRangeMap()->getComm());
408 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
409 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
410 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
411 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
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(),
418 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
422 originalTransferOp->getDomainMap()->getComm(),
423 stridedDoFullMap->getStridedBlockId(),
424 stridedDoFullMap->getOffset());
429 originalTransferOp->getDomainMap()->lib(),
430 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
433 originalTransferOp->getDomainMap()->getComm());
437 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
438 MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
439 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
440 MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
442 Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
443 for(
size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
444 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebP[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);
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define MueLu_maxAll(rcpComm, in, out)
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
An exception safe way to call the method 'Level::SetFactoryManager()'.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Statistics0
Print statistics that do not involve significant additional computation.