47 #ifndef MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_
48 #define MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_
52 #include <Xpetra_CrsMatrixWrap.hpp>
54 #include <Xpetra_MapExtractor.hpp>
56 #include <Xpetra_StridedMap.hpp>
57 #include <Xpetra_StridedMapFactory.hpp>
65 #include "MueLu_HierarchyUtils.hpp"
68 #include "MueLu_PerfUtils.hpp"
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
78 #undef SET_VALID_ENTRY
81 validParamList->
set<
RCP<const FactoryBase> >(
"Importer", Teuchos::null,
"Generating factory of the matrix Importer for rebalancing");
82 validParamList->
set<
RCP<const FactoryBase> >(
"SubImporters", Teuchos::null,
"Generating factory of the matrix sub-Importers for rebalancing");
84 return validParamList;
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 FactManager_.push_back(FactManager);
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 Input(coarseLevel,
"A");
96 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
97 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
101 coarseLevel.
DeclareInput(
"Importer", (*it)->GetFactory(
"Importer").get(),
this);
105 if (FactManager_.size() == 0) {
106 Input(coarseLevel,
"SubImporters");
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 RCP<Matrix> originalAc = Get<RCP<Matrix> >(coarseLevel,
"A");
120 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != bA->Cols(),
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Blocked operator has " << bA->Rows() <<
" and " << bA->Cols() <<
". We only support square matrices (with same number of blocks and columns).");
123 std::vector<GO> fullRangeMapVector;
124 std::vector<GO> fullDomainMapVector;
125 std::vector<RCP<const Map> > subBlockARangeMaps;
126 std::vector<RCP<const Map> > subBlockADomainMaps;
127 subBlockARangeMaps.reserve(bA->Rows());
128 subBlockADomainMaps.reserve(bA->Cols());
141 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
142 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
145 std::vector<RCP<Matrix> > subBlockRebA =
146 std::vector<RCP<Matrix> >(bA->Cols() * bA->Rows(), Teuchos::null);
150 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bA->Rows(), Teuchos::null);
151 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
153 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
158 importers[idx] = rebalanceImporter;
161 if (FactManager_.size() == 0) {
162 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
166 bool bRestrictComm =
false;
168 if (pL.
get<
bool>(
"repartition: use subcommunicators") ==
true)
169 bRestrictComm =
true;
173 XpetraList->
set(
"Restrict Communicator",
true);
175 XpetraList->
set(
"Restrict Communicator",
false);
184 for (
size_t i = 0; i < bA->Rows(); i++) {
185 for (
size_t j = 0; j < bA->Cols(); j++) {
189 std::stringstream ss;
190 ss <<
"Rebalancing matrix block A(" << i <<
"," << j <<
")";
195 if (importers[i] != Teuchos::null &&
196 importers[j] != Teuchos::null &&
197 Aij != Teuchos::null) {
210 TEUCHOS_TEST_FOR_EXCEPTION(rebAijImport.
is_null() ==
true,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Importer associated with block " << j <<
" is null.");
213 TEUCHOS_TEST_FOR_EXCEPTION(cAwij.is_null() ==
true,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Block (" << i <<
"," << j <<
") is not of type CrsMatrix. We cannot rebalanced (nested) operators.");
221 rebAij = MatrixFactory::Build(cAij, *(importers[i]), *(importers[j]), targetDomainMap, targetRangeMap, XpetraList);
228 subBlockRebA[i * bA->Cols() + j] = rebAij;
230 if (!rebAij.is_null()) {
232 if (rebalancedComm.
is_null()) rebalancedComm = rebAij->getRowMap()->getComm();
236 params->
set(
"printLoadBalancingInfo",
true);
237 std::stringstream ss2;
238 ss2 <<
"A(" << i <<
"," << j <<
") rebalanced:";
248 if (subBlockRebA[i * bA->Cols() + i].is_null() ==
false) {
249 RCP<Matrix> rebAii = subBlockRebA[i * bA->Cols() + i];
250 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(i, rangeMapExtractor->getThyraMode()));
252 if (orig_stridedRgMap != Teuchos::null) {
253 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
255 stridedRgMap = StridedMapFactory::Build(
256 bA->getRangeMap()->lib(),
259 rebAii->getRangeMap()->getIndexBase(),
262 orig_stridedRgMap->getStridedBlockId(),
263 orig_stridedRgMap->getOffset());
265 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(i, domainMapExtractor->getThyraMode()));
267 if (orig_stridedDoMap != Teuchos::null) {
268 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
270 stridedDoMap = StridedMapFactory::Build(
271 bA->getDomainMap()->lib(),
274 rebAii->getDomainMap()->getIndexBase(),
277 orig_stridedDoMap->getStridedBlockId(),
278 orig_stridedDoMap->getOffset());
282 stridedRgMap->removeEmptyProcesses();
283 stridedDoMap->removeEmptyProcesses();
286 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
287 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
290 if (rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
291 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
293 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
296 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.
begin(), nodeRangeMap.
end());
297 if (bThyraRangeGIDs ==
false)
298 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
300 subBlockADomainMaps.push_back(rebAii->getColMap(
"stridedMaps"));
303 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.
begin(), nodeDomainMap.
end());
304 if (bThyraDomainGIDs ==
false)
305 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
312 if (rebalancedComm == Teuchos::null) {
313 GetOStream(
Debug, -1) <<
"RebalanceBlockedAc: deactivate proc " << originalAc->getRowMap()->getComm()->getRank() << std::endl;
316 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
322 GO rangeIndexBase = bA->getRangeMap()->getIndexBase();
323 GO domainIndexBase = bA->getDomainMap()->getIndexBase();
325 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
328 if (stridedRgFullMap != Teuchos::null) {
329 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
331 StridedMapFactory::Build(
332 bA->getRangeMap()->lib(),
338 stridedRgFullMap->getStridedBlockId(),
339 stridedRgFullMap->getOffset());
343 bA->getRangeMap()->lib(),
349 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0, fullDomainMapVector.size());
353 if (stridedDoFullMap != Teuchos::null) {
354 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: full map in domain map extractor has no striding information! error.");
355 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
357 StridedMapFactory::Build(
358 bA->getDomainMap()->lib(),
364 stridedDoFullMap->getStridedBlockId(),
365 stridedDoFullMap->getOffset());
369 bA->getDomainMap()->lib(),
377 fullRangeMap->removeEmptyProcesses();
378 fullDomainMap->removeEmptyProcesses();
382 RCP<const MapExtractor> rebRangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockARangeMaps, bThyraRangeGIDs);
383 RCP<const MapExtractor> rebDomainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockADomainMaps, bThyraDomainGIDs);
385 TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->NumMaps() != rebRangeMapExtractor->NumMaps(),
Exceptions::RuntimeError,
386 "MueLu::RebalanceBlockedAc::Build: Rebalanced RangeMapExtractor has " << rebRangeMapExtractor->NumMaps()
387 <<
" sub maps. Original RangeMapExtractor has " << rangeMapExtractor->NumMaps() <<
". They must match!");
388 TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->NumMaps() != rebDomainMapExtractor->NumMaps(),
Exceptions::RuntimeError,
389 "MueLu::RebalanceBlockedAc::Build: Rebalanced DomainMapExtractor has " << rebDomainMapExtractor->NumMaps()
390 <<
" sub maps. Original DomainMapExtractor has " << domainMapExtractor->NumMaps() <<
". They must match!");
393 for (
size_t i = 0; i < bA->Rows(); i++) {
394 for (
size_t j = 0; j < bA->Cols(); j++) {
395 reb_bA->setMatrix(i, j, subBlockRebA[i * bA->Cols() + j]);
399 reb_bA->fillComplete();
401 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
406 if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
410 for (std::vector<
RCP<const FactoryBase> >::const_iterator it2 = rebalanceFacts_.begin(); it2 != rebalanceFacts_.end(); ++it2) {
411 GetOStream(
Runtime0) <<
"RebalanceBlockedAc: call rebalance factory " << (*it2).get() <<
": " << (*it2)->description() << std::endl;
412 (*it2)->CallBuild(coarseLevel);
417 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
419 rebalanceFacts_.push_back(factory);
Exception indicating invalid cast attempted.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
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)
#define SET_VALID_ENTRY(name)
Print additional debugging information.
One-liner description of what is happening.
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)
void AddRebalanceFactory(const RCP< const FactoryBase > &factory)
Add rebalancing factory in the end of list of rebalancing factories in RebalanceAcFactory.
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)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()