46 #ifndef MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_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> 68 #include "MueLu_HierarchyUtils.hpp" 73 #include "MueLu_PerfUtils.hpp" 77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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>
120 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"R");
123 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalTransferOp);
130 bool bRestrictComm =
false;
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());
158 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
159 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
172 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Rii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 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");
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());
190 params->
set(
"printLoadBalancingInfo",
true);
191 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
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()));
204 if(orig_stridedRgMap != Teuchos::null) {
205 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
207 stridedRgMap = StridedMapFactory::Build(
208 originalTransferOp->getRangeMap()->lib(),
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()));
220 if(orig_stridedDoMap != Teuchos::null) {
221 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
223 stridedDoMap = StridedMapFactory::Build(
224 originalTransferOp->getDomainMap()->lib(),
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();
243 if(rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
244 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
247 subBlockRebR.push_back(rebRii);
251 subBlockRRangeMaps.push_back(rangeMapii);
254 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.
begin(), nodeRangeMapii.
end());
255 if(bThyraRangeGIDs ==
false)
256 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
260 subBlockRDomainMaps.push_back(domainMapii);
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 <<
")";
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());
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());
306 if(stridedRgFullMap != Teuchos::null) {
307 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
309 StridedMapFactory::Build(
310 originalTransferOp->getRangeMap()->lib(),
315 originalTransferOp->getRangeMap()->getComm(),
316 stridedRgFullMap->getStridedBlockId(),
317 stridedRgFullMap->getOffset());
321 originalTransferOp->getRangeMap()->lib(),
325 originalTransferOp->getRangeMap()->getComm());
328 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
331 if(stridedDoFullMap != Teuchos::null) {
333 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
335 StridedMapFactory::Build(
336 originalTransferOp->getDomainMap()->lib(),
341 originalTransferOp->getDomainMap()->getComm(),
342 stridedDoFullMap->getStridedBlockId(),
343 stridedDoFullMap->getOffset());
348 originalTransferOp->getDomainMap()->lib(),
352 originalTransferOp->getDomainMap()->getComm());
356 fullRangeMap->removeEmptyProcesses();
357 fullDomainMap->removeEmptyProcesses();
362 MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
364 MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
367 for(
size_t i = 0; i<subBlockRRangeMaps.size(); 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.
T & get(const std::string &name, T def_value)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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)
Namespace for MueLu classes and methods.
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.
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()