47#ifndef MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_
48#define MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_
50#include <Xpetra_Matrix.hpp>
51#include <Xpetra_CrsMatrix.hpp>
52#include <Xpetra_CrsMatrixWrap.hpp>
53#include <Xpetra_MatrixFactory.hpp>
54#include <Xpetra_MapExtractor.hpp>
55#include <Xpetra_MapExtractorFactory.hpp>
56#include <Xpetra_StridedMap.hpp>
57#include <Xpetra_StridedMapFactory.hpp>
58#include <Xpetra_BlockedCrsMatrix.hpp>
60#include <Xpetra_VectorFactory.hpp>
65#include "MueLu_HierarchyUtils.hpp"
68#include "MueLu_PerfUtils.hpp"
69#include "MueLu_RAPFactory.hpp"
73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 RCP<ParameterList> validParamList = rcp(
new ParameterList());
77#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
81 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A for rebalancing");
82 validParamList->set<RCP<const FactoryBase> >(
"Importer", Teuchos::null,
"Generating factory of the matrix Importer for rebalancing");
83 validParamList->set<RCP<const FactoryBase> >(
"SubImporters", Teuchos::null,
"Generating factory of the matrix sub-Importers for rebalancing");
85 return validParamList;
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 FactManager_.push_back(FactManager);
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 Input(coarseLevel,
"A");
97 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
98 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
102 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
106 if(FactManager_.size() == 0) {
107 Input(coarseLevel,
"SubImporters");
112 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
118 RCP<Matrix> originalAc = Get< RCP<Matrix> >(coarseLevel,
"A");
120 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalAc);
121 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockAcFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
122 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).");
125 std::vector<GO> fullRangeMapVector;
126 std::vector<GO> fullDomainMapVector;
127 std::vector<RCP<const Map> > subBlockARangeMaps;
128 std::vector<RCP<const Map> > subBlockADomainMaps;
129 subBlockARangeMaps.reserve(bA->Rows());
130 subBlockADomainMaps.reserve(bA->Cols());
133 RCP<const MapExtractor> rangeMapExtractor = bA->getRangeMapExtractor();
134 RCP<const MapExtractor> domainMapExtractor = bA->getDomainMapExtractor();
143 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
144 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
147 std::vector<RCP<Matrix> > subBlockRebA =
148 std::vector<RCP<Matrix> >(bA->Cols() * bA->Rows(), Teuchos::null);
152 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bA->Rows(), Teuchos::null);
153 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
155 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
159 RCP<const Import> rebalanceImporter = coarseLevel.
Get<RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
160 importers[idx] = rebalanceImporter;
163 if(FactManager_.size() == 0) {
164 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
169 bool bRestrictComm =
false;
170 const ParameterList& pL = GetParameterList();
171 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
172 bRestrictComm =
true;
174 RCP<ParameterList> XpetraList = Teuchos::rcp(
new ParameterList());
176 XpetraList->set(
"Restrict Communicator",
true);
178 XpetraList->set(
"Restrict Communicator",
false);
182 RCP<const Teuchos::Comm<int> > rebalancedComm = Teuchos::null;
187 for(
size_t i=0; i<bA->Rows(); i++) {
188 for(
size_t j=0; j<bA->Cols(); j++) {
190 RCP<Matrix> Aij = bA->getMatrix(i, j);
192 std::stringstream ss; ss <<
"Rebalancing matrix block A(" << i <<
"," << j <<
")";
195 RCP<Matrix> rebAij = Teuchos::null;
197 if( importers[i] != Teuchos::null &&
198 importers[j] != Teuchos::null &&
199 Aij != Teuchos::null) {
200 RCP<const Map> targetRangeMap = importers[i]->getTargetMap();
201 RCP<const Map> targetDomainMap = importers[j]->getTargetMap();
208 RCP<Matrix> cAij = Aij;
211 Teuchos::RCP<const Import> rebAijImport = ImportFactory::Build(importers[j]->getTargetMap(),cAij->getColMap());
212 TEUCHOS_TEST_FOR_EXCEPTION(rebAijImport.is_null() ==
true,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Importer associated with block " << j <<
" is null.");
214 Teuchos::RCP<const CrsMatrixWrap> cAwij = Teuchos::rcp_dynamic_cast<const CrsMatrixWrap>(cAij);
215 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.");
216 Teuchos::RCP<CrsMatrix> cAmij = cAwij->getCrsMatrix();
223 rebAij = MatrixFactory::Build(cAij, *(importers[i]), *(importers[j]), targetDomainMap, targetRangeMap, XpetraList);
230 subBlockRebA[i*bA->Cols() + j] = rebAij;
232 if (!rebAij.is_null()) {
234 if(rebalancedComm.is_null()) rebalancedComm = rebAij->getRowMap()->getComm();
237 RCP<ParameterList> params = rcp(
new ParameterList());
238 params->set(
"printLoadBalancingInfo",
true);
239 std::stringstream ss2; ss2 <<
"A(" << i <<
"," << j <<
") rebalanced:";
249 if ( subBlockRebA[i*bA->Cols() + i].is_null() ==
false ) {
250 RCP<Matrix> rebAii = subBlockRebA[i*bA->Cols() + i];
251 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(i,rangeMapExtractor->getThyraMode()));
252 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
253 if(orig_stridedRgMap != Teuchos::null) {
254 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
255 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebAii->getRangeMap()->getLocalElementList();
256 stridedRgMap = StridedMapFactory::Build(
257 bA->getRangeMap()->lib(),
258 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
260 rebAii->getRangeMap()->getIndexBase(),
263 orig_stridedRgMap->getStridedBlockId(),
264 orig_stridedRgMap->getOffset());
266 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(i,domainMapExtractor->getThyraMode()));
267 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
268 if(orig_stridedDoMap != Teuchos::null) {
269 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
270 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebAii->getDomainMap()->getLocalElementList();
271 stridedDoMap = StridedMapFactory::Build(
272 bA->getDomainMap()->lib(),
273 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
275 rebAii->getDomainMap()->getIndexBase(),
278 orig_stridedDoMap->getStridedBlockId(),
279 orig_stridedDoMap->getOffset());
283 stridedRgMap->removeEmptyProcesses();
284 stridedDoMap->removeEmptyProcesses();
287 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
288 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
291 if(rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
292 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
294 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
295 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMap = rebAii->getRangeMap()->getLocalElementList();
297 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
298 if(bThyraRangeGIDs ==
false)
299 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
301 subBlockADomainMaps.push_back(rebAii->getColMap(
"stridedMaps"));
302 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMap = rebAii->getDomainMap()->getLocalElementList();
304 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
305 if(bThyraDomainGIDs ==
false)
306 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
313 if (rebalancedComm == Teuchos::null) {
314 GetOStream(
Debug,-1) <<
"RebalanceBlockedAc: deactivate proc " << originalAc->getRowMap()->getComm()->getRank() << std::endl;
316 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::null;
317 coarseLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
323 GO rangeIndexBase = bA->getRangeMap()->getIndexBase();
324 GO domainIndexBase = bA->getDomainMap()->getIndexBase();
326 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
327 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getFullMap());
328 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
329 if(stridedRgFullMap != Teuchos::null) {
330 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
332 StridedMapFactory::Build(
333 bA->getRangeMap()->lib(),
334 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
339 stridedRgFullMap->getStridedBlockId(),
340 stridedRgFullMap->getOffset());
344 bA->getRangeMap()->lib(),
345 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
350 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
352 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getFullMap());
353 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
354 if(stridedDoFullMap != Teuchos::null) {
355 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: full map in domain map extractor has no striding information! error.");
356 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
358 StridedMapFactory::Build(
359 bA->getDomainMap()->lib(),
360 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
365 stridedDoFullMap->getStridedBlockId(),
366 stridedDoFullMap->getOffset());
371 bA->getDomainMap()->lib(),
372 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
379 fullRangeMap->removeEmptyProcesses();
380 fullDomainMap->removeEmptyProcesses();
384 RCP<const MapExtractor> rebRangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockARangeMaps, bThyraRangeGIDs);
385 RCP<const MapExtractor> rebDomainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockADomainMaps, bThyraDomainGIDs);
387 TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->NumMaps() != rebRangeMapExtractor->NumMaps(),
Exceptions::RuntimeError,
388 "MueLu::RebalanceBlockedAc::Build: Rebalanced RangeMapExtractor has " << rebRangeMapExtractor->NumMaps()
389 <<
" sub maps. Original RangeMapExtractor has " << rangeMapExtractor->NumMaps() <<
". They must match!");
390 TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->NumMaps() != rebDomainMapExtractor->NumMaps(),
Exceptions::RuntimeError,
391 "MueLu::RebalanceBlockedAc::Build: Rebalanced DomainMapExtractor has " << rebDomainMapExtractor->NumMaps()
392 <<
" sub maps. Original DomainMapExtractor has " << domainMapExtractor->NumMaps() <<
". They must match!");
394 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::rcp(
new BlockedCrsMatrix(rebRangeMapExtractor,rebDomainMapExtractor,10));
395 for(
size_t i=0; i<bA->Rows(); i++) {
396 for(
size_t j=0; j<bA->Cols(); j++) {
397 reb_bA->setMatrix(i,j,subBlockRebA[i*bA->Cols() + j]);
401 reb_bA->fillComplete();
403 coarseLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
408 if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
412 for (std::vector<RCP<const FactoryBase> >::const_iterator it2 = rebalanceFacts_.begin(); it2 != rebalanceFacts_.end(); ++it2) {
413 GetOStream(
Runtime0) <<
"RebalanceBlockedAc: call rebalance factory " << (*it2).get() <<
": " << (*it2)->description() << std::endl;
414 (*it2)->CallBuild(coarseLevel);
419 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
421 rebalanceFacts_.push_back(factory);
#define SET_VALID_ENTRY(name)
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void AddRebalanceFactory(const RCP< const FactoryBase > &factory)
Add rebalancing factory in the end of list of rebalancing factories in RebalanceAcFactory.
An exception safe way to call the method 'Level::SetFactoryManager()'.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Namespace for MueLu classes and methods.
@ Debug
Print additional debugging information.
@ Runtime0
One-liner description of what is happening.
@ Statistics0
Print statistics that do not involve significant additional computation.