46 #ifndef MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ 49 #ifdef HAVE_MUELU_EXPERIMENTAL 51 #include <Teuchos_Tuple.hpp> 68 #include "MueLu_HierarchyUtils.hpp" 73 #include "MueLu_PerfUtils.hpp" 77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 RCP<ParameterList> validParamList = rcp(
new ParameterList());
81 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 83 #undef SET_VALID_ENTRY 87 validParamList->set< RCP<const FactoryBase> >(
"R", Teuchos::null,
"Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
89 return validParamList;
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 FactManager_.push_back(FactManager);
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 Input(coarseLevel,
"R");
101 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
102 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
106 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
107 coarseLevel.
DeclareInput(
"Nullspace",(*it)->GetFactory(
"Nullspace").get(),
this);
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
119 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
120 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"R");
122 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
124 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
126 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
127 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
130 bool bRestrictComm =
false;
131 const ParameterList& pL = GetParameterList();
132 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
133 bRestrictComm =
true;
142 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
143 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
146 std::vector<GO> fullRangeMapVector;
147 std::vector<GO> fullDomainMapVector;
148 std::vector<RCP<const Map> > subBlockRRangeMaps;
149 std::vector<RCP<const Map> > subBlockRDomainMaps;
150 subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows());
151 subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols());
153 std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
154 subBlockRebR.reserve(bOriginalTransferOp->Cols());
157 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
158 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
159 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
164 rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
167 Teuchos::RCP<Matrix> Rii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
170 TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs ==
true && Rii->getRowMap()->getMinAllGlobalIndex() != 0,
172 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Rii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
173 TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs ==
true && Rii->getColMap()->getMinAllGlobalIndex() != 0,
175 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Rii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
177 Teuchos::RCP<Matrix> rebRii;
178 if(rebalanceImporter != Teuchos::null) {
179 std::stringstream ss; ss <<
"Rebalancing restriction block R(" << curBlockId <<
"," << curBlockId <<
")";
182 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
186 rebRii = MatrixFactory::Build(Rii,*rebalanceImporter,dummy,rebalanceImporter->getTargetMap());
189 RCP<ParameterList> params = rcp(
new ParameterList());
190 params->set(
"printLoadBalancingInfo",
true);
191 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
195 RCP<ParameterList> params = rcp(
new ParameterList());
196 params->set(
"printLoadBalancingInfo",
true);
197 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
202 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
203 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
204 if(orig_stridedRgMap != Teuchos::null) {
205 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
206 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
207 stridedRgMap = StridedMapFactory::Build(
208 originalTransferOp->getRangeMap()->lib(),
209 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
211 rebRii->getRangeMap()->getIndexBase(),
213 originalTransferOp->getRangeMap()->getComm(),
214 orig_stridedRgMap->getStridedBlockId(),
215 orig_stridedRgMap->getOffset());
216 }
else stridedRgMap = Rii->getRangeMap();
218 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
219 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
220 if(orig_stridedDoMap != Teuchos::null) {
221 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
222 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
223 stridedDoMap = StridedMapFactory::Build(
224 originalTransferOp->getDomainMap()->lib(),
225 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
227 rebRii->getDomainMap()->getIndexBase(),
229 originalTransferOp->getDomainMap()->getComm(),
230 orig_stridedDoMap->getStridedBlockId(),
231 orig_stridedDoMap->getOffset());
232 }
else stridedDoMap = Rii->getDomainMap();
235 stridedRgMap->removeEmptyProcesses();
236 stridedDoMap->removeEmptyProcesses();
239 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
240 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
243 if(rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
244 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
247 subBlockRebR.push_back(rebRii);
250 Teuchos::RCP<const Map> rangeMapii = rebRii->getRowMap(
"stridedMaps");
251 subBlockRRangeMaps.push_back(rangeMapii);
252 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
254 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
255 if(bThyraRangeGIDs ==
false)
256 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
259 Teuchos::RCP<const Map> domainMapii = rebRii->getColMap(
"stridedMaps");
260 subBlockRDomainMaps.push_back(domainMapii);
261 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
263 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
264 if(bThyraDomainGIDs ==
false)
265 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
272 if(rebalanceImporter != Teuchos::null)
274 std::stringstream ss2; ss2 <<
"Rebalancing nullspace block(" << curBlockId <<
"," << curBlockId <<
")";
277 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
278 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
279 permutedNullspace->doImport(*nullspace, *rebalanceImporter,
Xpetra::INSERT);
283 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
285 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", permutedNullspace, (*it)->GetFactory(
"Nullspace").get());
289 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
290 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", nullspace, (*it)->GetFactory(
"Nullspace").get());
299 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
300 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
303 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
304 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getFullMap());
305 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
306 if(stridedRgFullMap != Teuchos::null) {
307 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
309 StridedMapFactory::Build(
310 originalTransferOp->getRangeMap()->lib(),
311 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
315 originalTransferOp->getRangeMap()->getComm(),
316 stridedRgFullMap->getStridedBlockId(),
317 stridedRgFullMap->getOffset());
321 originalTransferOp->getRangeMap()->lib(),
322 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
325 originalTransferOp->getRangeMap()->getComm());
328 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
329 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getFullMap());
330 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
331 if(stridedDoFullMap != Teuchos::null) {
332 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
333 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
335 StridedMapFactory::Build(
336 originalTransferOp->getDomainMap()->lib(),
337 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
341 originalTransferOp->getDomainMap()->getComm(),
342 stridedDoFullMap->getStridedBlockId(),
343 stridedDoFullMap->getOffset());
348 originalTransferOp->getDomainMap()->lib(),
349 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
352 originalTransferOp->getDomainMap()->getComm());
356 fullRangeMap->removeEmptyProcesses();
357 fullDomainMap->removeEmptyProcesses();
361 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
362 MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
363 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
364 MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
366 Teuchos::RCP<BlockedCrsMatrix> bRebR = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
367 for(
size_t i = 0; i<subBlockRRangeMaps.size(); i++) {
368 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebR[i]);
369 bRebR->setMatrix(i,i,crsOpii);
372 bRebR->fillComplete();
374 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()