46 #ifndef MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ 49 #ifdef HAVE_MUELU_EXPERIMENTAL 51 #include <Teuchos_Tuple.hpp> 53 #include "Xpetra_Vector.hpp" 54 #include "Xpetra_VectorFactory.hpp" 55 #include "Xpetra_MultiVector.hpp" 56 #include "Xpetra_MultiVectorFactory.hpp" 57 #include <Xpetra_Matrix.hpp> 58 #include <Xpetra_BlockedCrsMatrix.hpp> 59 #include <Xpetra_MapFactory.hpp> 60 #include <Xpetra_MapExtractor.hpp> 61 #include <Xpetra_MapExtractorFactory.hpp> 62 #include <Xpetra_MatrixFactory.hpp> 63 #include <Xpetra_Import.hpp> 64 #include <Xpetra_ImportFactory.hpp> 69 #include "MueLu_HierarchyUtils.hpp" 72 #include "MueLu_PerfUtils.hpp" 76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
80 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
81 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.");
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 84 #undef SET_VALID_ENTRY 90 return validParamList;
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 FactManager_.push_back(FactManager);
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
101 Input(coarseLevel,
"P");
102 Input(coarseLevel,
"A");
104 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
105 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
108 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
112 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
118 Teuchos::RCP<Matrix> nonrebCoarseA = Get< RCP<Matrix> >(coarseLevel,
"A");
120 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
121 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"P");
123 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
124 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
125 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
127 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
128 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
137 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
138 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
141 std::vector<GO> fullRangeMapVector;
142 std::vector<GO> fullDomainMapVector;
143 std::vector<RCP<const Map> > subBlockPRangeMaps;
144 std::vector<RCP<const Map> > subBlockPDomainMaps;
145 subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows());
146 subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols());
148 std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
149 subBlockRebP.reserve(bOriginalTransferOp->Rows());
152 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
153 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
154 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
159 rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
162 Teuchos::RCP<Matrix> Pii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
163 Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Pii);
164 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!");
167 TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs ==
true && Pii->getRowMap()->getMinAllGlobalIndex() != 0,
169 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Pii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
170 TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs ==
true && Pii->getColMap()->getMinAllGlobalIndex() != 0,
172 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Pii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
175 if(rebalanceImporter != Teuchos::null) {
176 std::stringstream ss; ss <<
"Rebalancing prolongator block P(" << curBlockId <<
"," << curBlockId <<
")";
190 RCP<const Import> newImporter;
192 SubFactoryMonitor subM(*
this,
"Rebalancing prolongator -- fast map replacement", coarseLevel);
193 newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
194 Pwii->getCrsMatrix()->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
197 RCP<ParameterList> params = rcp(
new ParameterList());
198 params->set(
"printLoadBalancingInfo",
true);
199 std::stringstream ss2; ss2 <<
"P(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
203 subBlockRebP.push_back(Pii);
206 RCP<ParameterList> params = rcp(
new ParameterList());
207 params->set(
"printLoadBalancingInfo",
true);
208 std::stringstream ss; ss <<
"P(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
211 subBlockRebP.push_back(Pii);
216 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
217 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
218 if(orig_stridedRgMap != Teuchos::null) {
219 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
220 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
221 stridedRgMap = StridedMapFactory::Build(
222 originalTransferOp->getRangeMap()->lib(),
223 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
225 Pii->getRangeMap()->getIndexBase(),
227 originalTransferOp->getRangeMap()->getComm(),
228 orig_stridedRgMap->getStridedBlockId(),
229 orig_stridedRgMap->getOffset());
230 }
else stridedRgMap = Pii->getRangeMap();
233 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
235 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
236 if(orig_stridedDoMap != Teuchos::null) {
237 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
238 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
239 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
240 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
242 Pii->getDomainMap()->getIndexBase(),
244 originalTransferOp->getDomainMap()->getComm(),
245 orig_stridedDoMap->getStridedBlockId(),
246 orig_stridedDoMap->getOffset());
247 }
else stridedDoMap = Pii->getDomainMap();
249 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
250 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
253 if(Pii->IsView(
"stridedMaps")) Pii->RemoveView(
"stridedMaps");
254 Pii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
257 Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap(
"stridedMaps");
258 subBlockPRangeMaps.push_back(rangeMapii);
259 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
261 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
262 if(bThyraRangeGIDs ==
false)
263 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
266 Teuchos::RCP<const Map> domainMapii = Pii->getColMap(
"stridedMaps");
267 subBlockPDomainMaps.push_back(domainMapii);
268 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
270 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
271 if(bThyraDomainGIDs ==
false)
272 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
279 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
280 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
282 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
283 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
284 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangePMapExtractor->getFullMap());
285 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
286 if(stridedRgFullMap != Teuchos::null) {
287 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
289 StridedMapFactory::Build(
290 originalTransferOp->getRangeMap()->lib(),
291 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
295 originalTransferOp->getRangeMap()->getComm(),
296 stridedRgFullMap->getStridedBlockId(),
297 stridedRgFullMap->getOffset());
301 originalTransferOp->getRangeMap()->lib(),
302 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
305 originalTransferOp->getRangeMap()->getComm());
308 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
309 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
310 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
311 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
312 if(stridedDoFullMap != Teuchos::null) {
313 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
314 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
316 StridedMapFactory::Build(
317 originalTransferOp->getDomainMap()->lib(),
318 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
322 originalTransferOp->getDomainMap()->getComm(),
323 stridedDoFullMap->getStridedBlockId(),
324 stridedDoFullMap->getOffset());
329 originalTransferOp->getDomainMap()->lib(),
330 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
333 originalTransferOp->getDomainMap()->getComm());
337 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
338 MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
339 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
340 MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
342 Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
343 for(
size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
344 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebP[i]);
345 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null,Xpetra::Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: block P" << i <<
" is not of type CrsMatrixWrap.");
346 bRebP->setMatrix(i,i,crsOpii);
348 bRebP->fillComplete();
350 Set(coarseLevel,
"P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
Exception indicating invalid cast attempted.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
Exception throws to report errors in the internal logical of the program.
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()