47 #ifndef MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_ 48 #define MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_ 50 #ifdef HAVE_MUELU_EXPERIMENTAL 52 #include <Xpetra_Matrix.hpp> 53 #include <Xpetra_CrsMatrix.hpp> 54 #include <Xpetra_CrsMatrixWrap.hpp> 55 #include <Xpetra_MatrixFactory.hpp> 56 #include <Xpetra_MapExtractor.hpp> 57 #include <Xpetra_MapExtractorFactory.hpp> 58 #include <Xpetra_StridedMap.hpp> 59 #include <Xpetra_StridedMapFactory.hpp> 60 #include <Xpetra_BlockedCrsMatrix.hpp> 62 #include <Xpetra_VectorFactory.hpp> 67 #include "MueLu_HierarchyUtils.hpp" 70 #include "MueLu_PerfUtils.hpp" 71 #include "MueLu_RAPFactory.hpp" 75 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 RCP<ParameterList> validParamList = rcp(
new ParameterList());
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 84 #undef SET_VALID_ENTRY 86 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A 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 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
115 RCP<Matrix> originalAc = Get< RCP<Matrix> >(coarseLevel,
"A");
117 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalAc);
118 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockAcFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
119 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).");
122 std::vector<GO> fullRangeMapVector;
123 std::vector<GO> fullDomainMapVector;
124 std::vector<RCP<const Map> > subBlockARangeMaps;
125 std::vector<RCP<const Map> > subBlockADomainMaps;
126 subBlockARangeMaps.reserve(bA->Rows());
127 subBlockADomainMaps.reserve(bA->Cols());
130 Teuchos::RCP<const MapExtractorClass> rangeMapExtractor = bA->getRangeMapExtractor();
131 Teuchos::RCP<const MapExtractorClass> domainMapExtractor = bA->getDomainMapExtractor();
140 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
141 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
144 std::vector<RCP<Matrix> > subBlockRebA =
145 std::vector<RCP<Matrix> >(bA->Cols() * bA->Rows(), Teuchos::null);
149 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bA->Rows(), Teuchos::null);
150 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
152 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
156 RCP<const Import> rebalanceImporter = coarseLevel.Get<RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
157 importers[idx] = rebalanceImporter;
162 bool bRestrictComm =
false;
163 const ParameterList& pL = GetParameterList();
164 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
165 bRestrictComm =
true;
167 RCP<ParameterList> XpetraList = Teuchos::rcp(
new ParameterList());
169 XpetraList->set(
"Restrict Communicator",
true);
171 XpetraList->set(
"Restrict Communicator",
false);
175 RCP<const Teuchos::Comm<int> > rebalancedComm = Teuchos::null;
180 for(
size_t i=0; i<bA->Rows(); i++) {
181 for(
size_t j=0; j<bA->Cols(); j++) {
183 RCP<Matrix> Aij = bA->getMatrix(i, j);
185 std::stringstream ss; ss <<
"Rebalancing matrix block A(" << i <<
"," << j <<
")";
188 RCP<Matrix> rebAij = Teuchos::null;
190 if( importers[i] != Teuchos::null &&
191 importers[j] != Teuchos::null &&
192 Aij != Teuchos::null) {
193 RCP<const Map> targetRangeMap = importers[i]->getTargetMap();
194 RCP<const Map> targetDomainMap = importers[j]->getTargetMap();
201 RCP<Matrix> cAij = Aij;
204 Teuchos::RCP<const Import> rebAijImport = ImportFactory::Build(importers[j]->getTargetMap(),cAij->getColMap());
205 TEUCHOS_TEST_FOR_EXCEPTION(rebAijImport.is_null() ==
true,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: Importer associated with block " << j <<
" is null.");
207 Teuchos::RCP<const CrsMatrixWrap> cAwij = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(cAij);
208 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.");
209 Teuchos::RCP<CrsMatrix> cAmij = cAwij->getCrsMatrix();
216 rebAij = MatrixFactory::Build(cAij, *(importers[i]), *(importers[j]), targetDomainMap, targetRangeMap, XpetraList);
223 subBlockRebA[i*bA->Cols() + j] = rebAij;
225 if (!rebAij.is_null()) {
227 if(rebalancedComm.is_null()) rebalancedComm = rebAij->getRowMap()->getComm();
230 RCP<ParameterList> params = rcp(
new ParameterList());
231 params->set(
"printLoadBalancingInfo",
true);
232 std::stringstream ss2; ss2 <<
"A(" << i <<
"," << j <<
") rebalanced:";
242 if ( subBlockRebA[i*bA->Cols() + i].is_null() == false ) {
243 RCP<Matrix> rebAii = subBlockRebA[i*bA->Cols() + i];
244 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(i,rangeMapExtractor->getThyraMode()));
245 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
246 if(orig_stridedRgMap != Teuchos::null) {
247 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
248 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebAii->getRangeMap()->getNodeElementList();
249 stridedRgMap = StridedMapFactory::Build(
250 bA->getRangeMap()->lib(),
251 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
253 rebAii->getRangeMap()->getIndexBase(),
256 orig_stridedRgMap->getStridedBlockId(),
257 orig_stridedRgMap->getOffset());
259 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(i,domainMapExtractor->getThyraMode()));
260 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
261 if(orig_stridedDoMap != Teuchos::null) {
262 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
263 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebAii->getDomainMap()->getNodeElementList();
264 stridedDoMap = StridedMapFactory::Build(
265 bA->getDomainMap()->lib(),
266 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
268 rebAii->getDomainMap()->getIndexBase(),
271 orig_stridedDoMap->getStridedBlockId(),
272 orig_stridedDoMap->getOffset());
276 stridedRgMap->removeEmptyProcesses();
277 stridedDoMap->removeEmptyProcesses();
280 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
281 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
284 if(rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
285 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
287 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
288 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMap = rebAii->getRangeMap()->getNodeElementList();
290 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
291 if(bThyraRangeGIDs ==
false)
292 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
294 subBlockADomainMaps.push_back(rebAii->getColMap(
"stridedMaps"));
295 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMap = rebAii->getDomainMap()->getNodeElementList();
297 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
298 if(bThyraDomainGIDs ==
false)
299 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
306 if (rebalancedComm == Teuchos::null) {
307 GetOStream(
Debug,-1) <<
"RebalanceBlockedAc: deactivate proc " << originalAc->getRowMap()->getComm()->getRank() << std::endl;
309 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::null;
310 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
316 GO rangeIndexBase = bA->getRangeMap()->getIndexBase();
317 GO domainIndexBase = bA->getDomainMap()->getIndexBase();
319 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
320 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getFullMap());
321 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
322 if(stridedRgFullMap != Teuchos::null) {
323 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
325 StridedMapFactory::Build(
326 bA->getRangeMap()->lib(),
327 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
332 stridedRgFullMap->getStridedBlockId(),
333 stridedRgFullMap->getOffset());
337 bA->getRangeMap()->lib(),
338 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
343 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
345 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getFullMap());
346 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
347 if(stridedDoFullMap != Teuchos::null) {
348 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockedAc::Build: full map in domain map extractor has no striding information! error.");
349 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
351 StridedMapFactory::Build(
352 bA->getDomainMap()->lib(),
353 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
358 stridedDoFullMap->getStridedBlockId(),
359 stridedDoFullMap->getOffset());
364 bA->getDomainMap()->lib(),
365 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
372 fullRangeMap->removeEmptyProcesses();
373 fullDomainMap->removeEmptyProcesses();
377 Teuchos::RCP<const MapExtractorClass> rebRangeMapExtractor = MapExtractorFactoryClass::Build(fullRangeMap, subBlockARangeMaps, bThyraRangeGIDs);
378 Teuchos::RCP<const MapExtractorClass> rebDomainMapExtractor = MapExtractorFactoryClass::Build(fullDomainMap, subBlockADomainMaps, bThyraDomainGIDs);
380 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!");
381 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!");
383 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::rcp(
new BlockedCrsMatrix(rebRangeMapExtractor,rebDomainMapExtractor,10));
384 for(
size_t i=0; i<bA->Rows(); i++) {
385 for(
size_t j=0; j<bA->Cols(); j++) {
387 reb_bA->setMatrix(i,j,subBlockRebA[i*bA->Cols() + j]);
391 reb_bA->fillComplete();
394 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
399 if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
403 for (std::vector<RCP<const FactoryBase> >::const_iterator it2 = rebalanceFacts_.begin(); it2 != rebalanceFacts_.end(); ++it2) {
404 GetOStream(
Runtime0) <<
"RebalanceBlockedAc: call rebalance factory " << (*it2).get() <<
": " << (*it2)->description() << std::endl;
405 (*it2)->CallBuild(coarseLevel);
410 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
417 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()