46 #ifndef MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_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> 66 #include "MueLu_HierarchyUtils.hpp" 71 #include "MueLu_PerfUtils.hpp" 75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
77 RCP<ParameterList> validParamList = rcp(
new ParameterList());
79 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 81 #undef SET_VALID_ENTRY 85 validParamList->set< RCP<const FactoryBase> >(
"R", Teuchos::null,
"Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
86 validParamList->set<RCP<const FactoryBase> >(
"Importer", Teuchos::null,
"Generating factory of the matrix Importer for rebalancing");
87 validParamList->set<RCP<const FactoryBase> >(
"SubImporters", Teuchos::null,
"Generating factory of the matrix sub-Importers for rebalancing");
88 validParamList->set<RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the Nullspace operator");
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>
100 Input(coarseLevel,
"R");
102 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
103 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
107 if(!UseSingleSourceImporters_) coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
108 coarseLevel.
DeclareInput(
"Nullspace",(*it)->GetFactory(
"Nullspace").get(),
this);
112 if(UseSingleSourceImporters_) {
113 Input(coarseLevel,
"SubImporters");
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
127 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
128 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"R");
130 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
131 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalTransferOp);
132 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
134 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
135 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
138 bool bRestrictComm =
false;
139 const ParameterList& pL = GetParameterList();
140 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
141 bRestrictComm =
true;
150 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
151 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
154 std::vector<GO> fullRangeMapVector;
155 std::vector<GO> fullDomainMapVector;
156 std::vector<RCP<const Map> > subBlockRRangeMaps;
157 std::vector<RCP<const Map> > subBlockRDomainMaps;
158 subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows());
159 subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols());
161 std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
162 subBlockRebR.reserve(bOriginalTransferOp->Cols());
164 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
165 if(UseSingleSourceImporters_) {
166 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
170 Teuchos::RCP<const Import> rebalanceImporter;
171 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
172 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
177 if(UseSingleSourceImporters_) rebalanceImporter = importers[curBlockId];
178 else rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
181 Teuchos::RCP<Matrix> Rii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
184 TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs ==
true && Rii->getRowMap()->getMinAllGlobalIndex() != 0,
186 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Rii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
187 TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs ==
true && Rii->getColMap()->getMinAllGlobalIndex() != 0,
189 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Rii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
191 Teuchos::RCP<Matrix> rebRii;
192 if(rebalanceImporter != Teuchos::null) {
193 std::stringstream ss; ss <<
"Rebalancing restriction block R(" << curBlockId <<
"," << curBlockId <<
")";
196 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
200 rebRii = MatrixFactory::Build(Rii,*rebalanceImporter,dummy,rebalanceImporter->getTargetMap());
203 RCP<ParameterList> params = rcp(
new ParameterList());
204 params->set(
"printLoadBalancingInfo",
true);
205 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
209 RCP<ParameterList> params = rcp(
new ParameterList());
210 params->set(
"printLoadBalancingInfo",
true);
211 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
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 = rebRii->getRangeMap()->getNodeElementList();
221 stridedRgMap = StridedMapFactory::Build(
222 originalTransferOp->getRangeMap()->lib(),
223 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
225 rebRii->getRangeMap()->getIndexBase(),
227 originalTransferOp->getRangeMap()->getComm(),
228 orig_stridedRgMap->getStridedBlockId(),
229 orig_stridedRgMap->getOffset());
230 }
else stridedRgMap = Rii->getRangeMap();
232 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
233 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
234 if(orig_stridedDoMap != Teuchos::null) {
235 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
236 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
237 stridedDoMap = StridedMapFactory::Build(
238 originalTransferOp->getDomainMap()->lib(),
239 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
241 rebRii->getDomainMap()->getIndexBase(),
243 originalTransferOp->getDomainMap()->getComm(),
244 orig_stridedDoMap->getStridedBlockId(),
245 orig_stridedDoMap->getOffset());
246 }
else stridedDoMap = Rii->getDomainMap();
249 stridedRgMap->removeEmptyProcesses();
250 stridedDoMap->removeEmptyProcesses();
253 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
254 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
257 if(rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
258 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
261 subBlockRebR.push_back(rebRii);
264 Teuchos::RCP<const Map> rangeMapii = rebRii->getRowMap(
"stridedMaps");
265 subBlockRRangeMaps.push_back(rangeMapii);
266 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
268 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
269 if(bThyraRangeGIDs ==
false)
270 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
273 Teuchos::RCP<const Map> domainMapii = rebRii->getColMap(
"stridedMaps");
274 subBlockRDomainMaps.push_back(domainMapii);
275 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
277 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
278 if(bThyraDomainGIDs ==
false)
279 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
286 if(rebalanceImporter != Teuchos::null)
288 std::stringstream ss2; ss2 <<
"Rebalancing nullspace block(" << curBlockId <<
"," << curBlockId <<
")";
291 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
292 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
293 permutedNullspace->doImport(*nullspace, *rebalanceImporter, Xpetra::INSERT);
297 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
299 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", permutedNullspace, (*it)->GetFactory(
"Nullspace").get());
303 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
304 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", nullspace, (*it)->GetFactory(
"Nullspace").get());
313 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
314 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
317 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
318 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getFullMap());
319 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
320 if(stridedRgFullMap != Teuchos::null) {
321 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
323 StridedMapFactory::Build(
324 originalTransferOp->getRangeMap()->lib(),
325 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
329 originalTransferOp->getRangeMap()->getComm(),
330 stridedRgFullMap->getStridedBlockId(),
331 stridedRgFullMap->getOffset());
335 originalTransferOp->getRangeMap()->lib(),
336 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
339 originalTransferOp->getRangeMap()->getComm());
342 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
343 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getFullMap());
344 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
345 if(stridedDoFullMap != Teuchos::null) {
346 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
347 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
349 StridedMapFactory::Build(
350 originalTransferOp->getDomainMap()->lib(),
351 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
355 originalTransferOp->getDomainMap()->getComm(),
356 stridedDoFullMap->getStridedBlockId(),
357 stridedDoFullMap->getOffset());
362 originalTransferOp->getDomainMap()->lib(),
363 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
366 originalTransferOp->getDomainMap()->getComm());
370 fullRangeMap->removeEmptyProcesses();
371 fullDomainMap->removeEmptyProcesses();
375 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
376 MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
377 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
378 MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
380 Teuchos::RCP<BlockedCrsMatrix> bRebR = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
381 for(
size_t i = 0; i<subBlockRRangeMaps.size(); i++) {
382 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebR[i]);
383 bRebR->setMatrix(i,i,crsOpii);
386 bRebR->fillComplete();
388 Set(coarseLevel,
"R", Teuchos::rcp_dynamic_cast<Matrix>(bRebR));
Exception indicating invalid cast attempted.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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)
#define SET_VALID_ENTRY(name)
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()