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)) 79 #undef SET_VALID_ENTRY 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()->getNodeElementList();
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()->getNodeElementList();
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()->getNodeElementList();
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()->getNodeElementList();
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);
Exception indicating invalid cast attempted.
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 AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define SET_VALID_ENTRY(name)
Print additional debugging information.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
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 Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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()