46 #ifndef MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_
47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_
49 #include <Teuchos_Tuple.hpp>
52 #include "Xpetra_MultiVectorFactory.hpp"
57 #include <Xpetra_MapFactory.hpp>
58 #include <Xpetra_MapExtractor.hpp>
66 #include "MueLu_HierarchyUtils.hpp"
71 #include "MueLu_PerfUtils.hpp"
75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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");
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");
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
124 originalTransferOp = Get<RCP<Matrix> >(coarseLevel,
"R");
134 bool bRestrictComm =
false;
136 if (pL.
get<
bool>(
"repartition: use subcommunicators") ==
true)
137 bRestrictComm =
true;
146 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
147 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
150 std::vector<GO> fullRangeMapVector;
151 std::vector<GO> fullDomainMapVector;
152 std::vector<RCP<const Map> > subBlockRRangeMaps;
153 std::vector<RCP<const Map> > subBlockRDomainMaps;
154 subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows());
155 subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols());
157 std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
158 subBlockRebR.reserve(bOriginalTransferOp->Cols());
160 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
161 if (UseSingleSourceImporters_) {
162 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
167 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
168 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
173 if (UseSingleSourceImporters_)
174 rebalanceImporter = importers[curBlockId];
182 TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs ==
true && Rii->getRowMap()->getMinAllGlobalIndex() != 0,
184 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Rii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
185 TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs ==
true && Rii->getColMap()->getMinAllGlobalIndex() != 0,
187 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Rii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
190 if (rebalanceImporter != Teuchos::null) {
191 std::stringstream ss;
192 ss <<
"Rebalancing restriction block R(" << curBlockId <<
"," << curBlockId <<
")";
195 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
199 rebRii = MatrixFactory::Build(Rii, *rebalanceImporter, dummy, rebalanceImporter->getTargetMap());
203 params->
set(
"printLoadBalancingInfo",
true);
204 std::stringstream ss2;
205 ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
210 params->
set(
"printLoadBalancingInfo",
true);
211 std::stringstream ss2;
212 ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
217 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId), rangeMapExtractor->getThyraMode()));
219 if (orig_stridedRgMap != Teuchos::null) {
220 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
222 stridedRgMap = StridedMapFactory::Build(
223 originalTransferOp->getRangeMap()->lib(),
226 rebRii->getRangeMap()->getIndexBase(),
228 originalTransferOp->getRangeMap()->getComm(),
229 orig_stridedRgMap->getStridedBlockId(),
230 orig_stridedRgMap->getOffset());
232 stridedRgMap = Rii->getRangeMap();
234 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId), domainMapExtractor->getThyraMode()));
236 if (orig_stridedDoMap != Teuchos::null) {
237 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
239 stridedDoMap = StridedMapFactory::Build(
240 originalTransferOp->getDomainMap()->lib(),
243 rebRii->getDomainMap()->getIndexBase(),
245 originalTransferOp->getDomainMap()->getComm(),
246 orig_stridedDoMap->getStridedBlockId(),
247 orig_stridedDoMap->getOffset());
249 stridedDoMap = Rii->getDomainMap();
252 stridedRgMap->removeEmptyProcesses();
253 stridedDoMap->removeEmptyProcesses();
256 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
257 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
260 if (rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
261 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
264 subBlockRebR.push_back(rebRii);
268 subBlockRRangeMaps.push_back(rangeMapii);
271 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.
begin(), nodeRangeMapii.
end());
272 if (bThyraRangeGIDs ==
false)
273 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
277 subBlockRDomainMaps.push_back(domainMapii);
280 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.
begin(), nodeDomainMapii.
end());
281 if (bThyraDomainGIDs ==
false)
282 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
289 if (rebalanceImporter != Teuchos::null) {
290 std::stringstream ss2;
291 ss2 <<
"Rebalancing nullspace block(" << curBlockId <<
"," << curBlockId <<
")";
295 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
296 permutedNullspace->doImport(*nullspace, *rebalanceImporter,
Xpetra::INSERT);
300 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
302 coarseLevel.Set<
RCP<MultiVector> >(
"Nullspace", permutedNullspace, (*it)->GetFactory(
"Nullspace").
get());
307 coarseLevel.Set<
RCP<MultiVector> >(
"Nullspace", nullspace, (*it)->GetFactory(
"Nullspace").
get());
316 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
317 GO domainIndexBase = originalTransferOp->getDomainMap()->getIndexBase();
320 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
323 if (stridedRgFullMap != Teuchos::null) {
324 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
326 StridedMapFactory::Build(
327 originalTransferOp->getRangeMap()->lib(),
332 originalTransferOp->getRangeMap()->getComm(),
333 stridedRgFullMap->getStridedBlockId(),
334 stridedRgFullMap->getOffset());
338 originalTransferOp->getRangeMap()->lib(),
342 originalTransferOp->getRangeMap()->getComm());
345 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0, fullDomainMapVector.size());
348 if (stridedDoFullMap != Teuchos::null) {
349 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
350 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
352 StridedMapFactory::Build(
353 originalTransferOp->getDomainMap()->lib(),
358 originalTransferOp->getDomainMap()->getComm(),
359 stridedDoFullMap->getStridedBlockId(),
360 stridedDoFullMap->getOffset());
364 originalTransferOp->getDomainMap()->lib(),
368 originalTransferOp->getDomainMap()->getComm());
372 fullRangeMap->removeEmptyProcesses();
373 fullDomainMap->removeEmptyProcesses();
378 MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
380 MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
383 for (
size_t i = 0; i < subBlockRRangeMaps.size(); i++) {
385 bRebR->setMatrix(i, i, crsOpii);
388 bRebR->fillComplete();
390 Set(coarseLevel,
"R", Teuchos::rcp_dynamic_cast<Matrix>(bRebR));
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Exception indicating invalid cast attempted.
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)
void sort(View &view, const size_t &size)
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)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define SET_VALID_ENTRY(name)
Exception throws to report errors in the internal logical of the program.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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()