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>
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
80 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 82 #undef SET_VALID_ENTRY 84 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A for rebalancing");
85 validParamList->set<RCP<const FactoryBase> >(
"Importer", Teuchos::null,
"Generating factory of the matrix Importer for rebalancing");
86 validParamList->set<RCP<const FactoryBase> >(
"SubImporters", Teuchos::null,
"Generating factory of the matrix sub-Importers for rebalancing");
88 return validParamList;
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 FactManager_.push_back(FactManager);
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 Input(coarseLevel,
"A");
100 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
101 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
105 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
109 if(FactManager_.size() == 0) {
110 Input(coarseLevel,
"SubImporters");
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
121 RCP<Matrix> originalAc = Get< RCP<Matrix> >(coarseLevel,
"A");
123 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalAc);
124 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockAcFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
125 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).");
128 std::vector<GO> fullRangeMapVector;
129 std::vector<GO> fullDomainMapVector;
130 std::vector<RCP<const Map> > subBlockARangeMaps;
131 std::vector<RCP<const Map> > subBlockADomainMaps;
132 subBlockARangeMaps.reserve(bA->Rows());
133 subBlockADomainMaps.reserve(bA->Cols());
136 Teuchos::RCP<const MapExtractorClass> rangeMapExtractor = bA->getRangeMapExtractor();
137 Teuchos::RCP<const MapExtractorClass> domainMapExtractor = bA->getDomainMapExtractor();
146 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
147 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
150 std::vector<RCP<Matrix> > subBlockRebA =
151 std::vector<RCP<Matrix> >(bA->Cols() * bA->Rows(), Teuchos::null);
155 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bA->Rows(), Teuchos::null);
156 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
158 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
162 RCP<const Import> rebalanceImporter = coarseLevel.Get<RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
163 importers[idx] = rebalanceImporter;
166 if(FactManager_.size() == 0) {
167 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
172 bool bRestrictComm =
false;
173 const ParameterList& pL = GetParameterList();
174 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
175 bRestrictComm =
true;
177 RCP<ParameterList> XpetraList = Teuchos::rcp(
new ParameterList());
179 XpetraList->set(
"Restrict Communicator",
true);
181 XpetraList->set(
"Restrict Communicator",
false);
185 RCP<const Teuchos::Comm<int> > rebalancedComm = Teuchos::null;
190 for(
size_t i=0; i<bA->Rows(); i++) {
191 for(
size_t j=0; j<bA->Cols(); j++) {
193 RCP<Matrix> Aij = bA->getMatrix(i, j);
195 std::stringstream ss; ss <<
"Rebalancing matrix block A(" << i <<
"," << j <<
")";
198 RCP<Matrix> rebAij = Teuchos::null;
200 if( importers[i] != Teuchos::null &&
201 importers[j] != Teuchos::null &&
202 Aij != Teuchos::null) {
203 RCP<const Map> targetRangeMap = importers[i]->getTargetMap();
204 RCP<const Map> targetDomainMap = importers[j]->getTargetMap();
211 RCP<Matrix> cAij = Aij;
214 Teuchos::RCP<const Import> rebAijImport = ImportFactory::Build(importers[j]->getTargetMap(),cAij->getColMap());
215 TEUCHOS_TEST_FOR_EXCEPTION(rebAijImport.is_null() ==
true,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Importer associated with block " << j <<
" is null.");
217 Teuchos::RCP<const CrsMatrixWrap> cAwij = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(cAij);
218 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.");
219 Teuchos::RCP<CrsMatrix> cAmij = cAwij->getCrsMatrix();
226 rebAij = MatrixFactory::Build(cAij, *(importers[i]), *(importers[j]), targetDomainMap, targetRangeMap, XpetraList);
233 subBlockRebA[i*bA->Cols() + j] = rebAij;
235 if (!rebAij.is_null()) {
237 if(rebalancedComm.is_null()) rebalancedComm = rebAij->getRowMap()->getComm();
240 RCP<ParameterList> params = rcp(
new ParameterList());
241 params->set(
"printLoadBalancingInfo",
true);
242 std::stringstream ss2; ss2 <<
"A(" << i <<
"," << j <<
") rebalanced:";
252 if ( subBlockRebA[i*bA->Cols() + i].is_null() == false ) {
253 RCP<Matrix> rebAii = subBlockRebA[i*bA->Cols() + i];
254 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(i,rangeMapExtractor->getThyraMode()));
255 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
256 if(orig_stridedRgMap != Teuchos::null) {
257 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
258 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebAii->getRangeMap()->getNodeElementList();
259 stridedRgMap = StridedMapFactory::Build(
260 bA->getRangeMap()->lib(),
261 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
263 rebAii->getRangeMap()->getIndexBase(),
266 orig_stridedRgMap->getStridedBlockId(),
267 orig_stridedRgMap->getOffset());
269 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(i,domainMapExtractor->getThyraMode()));
270 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
271 if(orig_stridedDoMap != Teuchos::null) {
272 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
273 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebAii->getDomainMap()->getNodeElementList();
274 stridedDoMap = StridedMapFactory::Build(
275 bA->getDomainMap()->lib(),
276 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
278 rebAii->getDomainMap()->getIndexBase(),
281 orig_stridedDoMap->getStridedBlockId(),
282 orig_stridedDoMap->getOffset());
286 stridedRgMap->removeEmptyProcesses();
287 stridedDoMap->removeEmptyProcesses();
290 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
291 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
294 if(rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
295 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
297 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
298 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMap = rebAii->getRangeMap()->getNodeElementList();
300 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
301 if(bThyraRangeGIDs ==
false)
302 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
304 subBlockADomainMaps.push_back(rebAii->getColMap(
"stridedMaps"));
305 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMap = rebAii->getDomainMap()->getNodeElementList();
307 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
308 if(bThyraDomainGIDs ==
false)
309 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
316 if (rebalancedComm == Teuchos::null) {
317 GetOStream(
Debug,-1) <<
"RebalanceBlockedAc: deactivate proc " << originalAc->getRowMap()->getComm()->getRank() << std::endl;
319 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::null;
320 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
326 GO rangeIndexBase = bA->getRangeMap()->getIndexBase();
327 GO domainIndexBase = bA->getDomainMap()->getIndexBase();
329 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
330 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getFullMap());
331 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
332 if(stridedRgFullMap != Teuchos::null) {
333 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
335 StridedMapFactory::Build(
336 bA->getRangeMap()->lib(),
337 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
342 stridedRgFullMap->getStridedBlockId(),
343 stridedRgFullMap->getOffset());
347 bA->getRangeMap()->lib(),
348 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
353 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
355 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getFullMap());
356 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
357 if(stridedDoFullMap != Teuchos::null) {
358 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: full map in domain map extractor has no striding information! error.");
359 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
361 StridedMapFactory::Build(
362 bA->getDomainMap()->lib(),
363 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
368 stridedDoFullMap->getStridedBlockId(),
369 stridedDoFullMap->getOffset());
374 bA->getDomainMap()->lib(),
375 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
382 fullRangeMap->removeEmptyProcesses();
383 fullDomainMap->removeEmptyProcesses();
387 Teuchos::RCP<const MapExtractorClass> rebRangeMapExtractor = MapExtractorFactoryClass::Build(fullRangeMap, subBlockARangeMaps, bThyraRangeGIDs);
388 Teuchos::RCP<const MapExtractorClass> rebDomainMapExtractor = MapExtractorFactoryClass::Build(fullDomainMap, subBlockADomainMaps, bThyraDomainGIDs);
390 TEUCHOS_TEST_FOR_EXCEPTION(rangeMapExtractor->NumMaps() != rebRangeMapExtractor->NumMaps(),
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: Rebalanced RangeMapExtractor has " << rebRangeMapExtractor <<
" sub maps. Original RangeMapExtractor has " << rangeMapExtractor->NumMaps() <<
". They must match!");
391 TEUCHOS_TEST_FOR_EXCEPTION(domainMapExtractor->NumMaps() != rebDomainMapExtractor->NumMaps(),
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: Rebalanced DomainMapExtractor has " << rebDomainMapExtractor <<
" sub maps. Original DomainMapExtractor has " << domainMapExtractor->NumMaps() <<
". They must match!");
393 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::rcp(
new BlockedCrsMatrix(rebRangeMapExtractor,rebDomainMapExtractor,10));
394 for(
size_t i=0; i<bA->Rows(); i++) {
395 for(
size_t j=0; j<bA->Cols(); j++) {
397 reb_bA->setMatrix(i,j,subBlockRebA[i*bA->Cols() + j]);
401 reb_bA->fillComplete();
404 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
409 if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
413 for (std::vector<RCP<const FactoryBase> >::const_iterator it2 = rebalanceFacts_.begin(); it2 != rebalanceFacts_.end(); ++it2) {
414 GetOStream(
Runtime0) <<
"RebalanceBlockedAc: call rebalance factory " << (*it2).get() <<
": " << (*it2)->description() << std::endl;
415 (*it2)->CallBuild(coarseLevel);
420 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
427 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()
RebalanceBlockAcFactory()