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>
80 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 82 #undef SET_VALID_ENTRY 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>
121 RCP<Matrix> originalAc = Get< RCP<Matrix> >(coarseLevel,
"A");
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());
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) {
163 importers[idx] = rebalanceImporter;
166 if(FactManager_.size() == 0) {
167 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
172 bool bRestrictComm =
false;
174 if (pL.
get<
bool>(
"repartition: use subcommunicators") ==
true)
175 bRestrictComm =
true;
179 XpetraList->
set(
"Restrict Communicator",
true);
181 XpetraList->
set(
"Restrict Communicator",
false);
190 for(
size_t i=0; i<bA->Rows(); i++) {
191 for(
size_t j=0; j<bA->Cols(); j++) {
195 std::stringstream ss; ss <<
"Rebalancing matrix block A(" << i <<
"," << j <<
")";
200 if( importers[i] != Teuchos::null &&
201 importers[j] != Teuchos::null &&
202 Aij != Teuchos::null) {
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();
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()));
256 if(orig_stridedRgMap != Teuchos::null) {
257 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
259 stridedRgMap = StridedMapFactory::Build(
260 bA->getRangeMap()->lib(),
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()));
271 if(orig_stridedDoMap != Teuchos::null) {
272 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
274 stridedDoMap = StridedMapFactory::Build(
275 bA->getDomainMap()->lib(),
278 rebAii->getDomainMap()->getIndexBase(),
281 orig_stridedDoMap->getStridedBlockId(),
282 orig_stridedDoMap->getOffset());
286 stridedRgMap->removeEmptyProcesses();
287 stridedDoMap->removeEmptyProcesses();
294 if(rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
295 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
297 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
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"));
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;
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());
332 if(stridedRgFullMap != Teuchos::null) {
333 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
335 StridedMapFactory::Build(
336 bA->getRangeMap()->lib(),
342 stridedRgFullMap->getStridedBlockId(),
343 stridedRgFullMap->getOffset());
347 bA->getRangeMap()->lib(),
353 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
357 if(stridedDoFullMap != Teuchos::null) {
359 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
361 StridedMapFactory::Build(
362 bA->getDomainMap()->lib(),
368 stridedDoFullMap->getStridedBlockId(),
369 stridedDoFullMap->getOffset());
374 bA->getDomainMap()->lib(),
382 fullRangeMap->removeEmptyProcesses();
383 fullDomainMap->removeEmptyProcesses();
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!");
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.
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.
Namespace for MueLu classes and methods.
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 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()