53 #ifndef MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_ 54 #define MUELU_REBALANCEBLOCKACFACTORY_DEF_HPP_ 56 #ifdef HAVE_MUELU_EXPERIMENTAL 58 #include <Xpetra_Matrix.hpp> 59 #include <Xpetra_CrsMatrix.hpp> 60 #include <Xpetra_CrsMatrixWrap.hpp> 61 #include <Xpetra_MatrixFactory.hpp> 62 #include <Xpetra_MapExtractor.hpp> 63 #include <Xpetra_MapExtractorFactory.hpp> 64 #include <Xpetra_StridedMap.hpp> 65 #include <Xpetra_StridedMapFactory.hpp> 66 #include <Xpetra_BlockedCrsMatrix.hpp> 71 #include "MueLu_HierarchyHelpers.hpp" 74 #include "MueLu_PerfUtils.hpp" 75 #include "MueLu_RAPFactory.hpp" 79 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 RCP<ParameterList> validParamList = rcp(
new ParameterList());
86 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 88 #undef SET_VALID_ENTRY 90 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A for rebalancing");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 FactManager_.push_back(FactManager);
100 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 Input(coarseLevel,
"A");
104 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
105 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
109 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
119 RCP<Matrix> originalAc = Get< RCP<Matrix> >(coarseLevel,
"A");
121 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bA = Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalAc);
122 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
125 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != 2,
Exceptions::RuntimeError,
"MueLu::BlockedPFactory::Build: number of block rows of A is not equal 2. error.");
126 TEUCHOS_TEST_FOR_EXCEPTION(bA->Cols() != 2,
Exceptions::RuntimeError,
"MueLu::BlockedPFactory::Build: number of block columns of A is not equal 2. error.");
129 Teuchos::RCP<const MapExtractorClass> rangeMapExtractor = bA->getRangeMapExtractor();
130 Teuchos::RCP<const MapExtractorClass> domainMapExtractor = bA->getDomainMapExtractor();
132 std::vector<GO> fullRangeMapVector;
133 std::vector<GO> fullDomainMapVector;
135 std::vector<RCP<const Map> > subBlockARangeMaps;
136 std::vector<RCP<const Map> > subBlockADomainMaps;
137 subBlockARangeMaps.reserve(bA->Rows());
138 subBlockADomainMaps.reserve(bA->Cols());
140 std::vector<Teuchos::RCP<Matrix> > subBlockRebA;
141 subBlockRebA.reserve(bA->Cols() * bA->Rows());
143 for(
size_t i=0; i<bA->Rows(); i++) {
144 for(
size_t j=0; j<bA->Cols(); j++) {
146 Teuchos::RCP<CrsMatrix> Amij = bA->getMatrix(i, j);
147 Teuchos::RCP<CrsMatrixWrap> Awij = Teuchos::rcp(
new CrsMatrixWrap(Amij));
148 Teuchos::RCP<Matrix> Aij = Teuchos::rcp_dynamic_cast<Matrix>(Awij);
150 subBlockRebA.push_back(Aij);
154 size_t curBlockId = 0;
155 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
156 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
160 Teuchos::RCP<const Import> rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
165 Teuchos::RCP<Matrix> Aii = subBlockRebA[curBlockId*bA->Cols() + curBlockId];
167 Teuchos::RCP<Matrix> rebAii;
168 if(rebalanceImporter != Teuchos::null) {
169 std::stringstream ss; ss <<
"Rebalancing matrix block A(" << curBlockId <<
"," << curBlockId <<
")";
171 RCP<const Map> targetMap = rebalanceImporter->getTargetMap();
175 ParameterList XpetraList;
178 XpetraList.set(
"Restrict Communicator",
false );
181 rebAii = MatrixFactory::Build(Aii, *rebalanceImporter, targetMap, targetMap, rcp(&XpetraList,
false));
186 if (!rebAii.is_null()) {
187 RCP<ParameterList> params = rcp(
new ParameterList());
188 params->set(
"printLoadBalancingInfo",
true);
189 std::stringstream ss2; ss2 <<
"A(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
205 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgAMapExtractor = bA->getRangeMapExtractor();
206 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rgAMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
207 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
208 if(orig_stridedRgMap != Teuchos::null) {
209 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
210 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebAii->getRangeMap()->getNodeElementList();
211 stridedRgMap = StridedMapFactory::Build(
212 bA->getRangeMap()->lib(),
213 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
215 rebAii->getRangeMap()->getIndexBase(),
217 bA->getRangeMap()->getComm(),
218 orig_stridedRgMap->getStridedBlockId(),
219 orig_stridedRgMap->getOffset());
221 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doAMapExtractor = bA->getDomainMapExtractor();
222 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(doAMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
223 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
224 if(orig_stridedDoMap != Teuchos::null) {
225 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
226 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebAii->getDomainMap()->getNodeElementList();
227 stridedDoMap = StridedMapFactory::Build(
228 bA->getDomainMap()->lib(),
229 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
231 rebAii->getDomainMap()->getIndexBase(),
233 bA->getDomainMap()->getComm(),
234 orig_stridedDoMap->getStridedBlockId(),
235 orig_stridedDoMap->getOffset());
238 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
239 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockAcFactory::Build: failed to generate striding information. error.");
242 if(rebAii->IsView(
"stridedMaps")) rebAii->RemoveView(
"stridedMaps");
243 rebAii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
245 subBlockRebA[curBlockId*bA->Cols() + curBlockId] = rebAii;
248 for(
size_t j=0; j<bA->Cols(); j++) {
249 if(j==curBlockId)
continue;
252 Teuchos::RCP<Matrix> Aij = subBlockRebA[curBlockId*bA->Cols() + j];
254 Teuchos::RCP<Matrix> rebAij;
255 if(rebalanceImporter!=Teuchos::null) {
256 std::stringstream ss3; ss3 <<
"Rebalancing matrix block A(" << curBlockId <<
"," << j <<
")";
258 RCP<const Map> targetMap = rebalanceImporter->getTargetMap();
260 Teuchos::RCP<Map> dummy;
261 rebAij = MatrixFactory::Build(Aij, *rebalanceImporter, dummy, targetMap);
264 if (!rebAij.is_null() && Aij->IsView(
"stridedMaps"))
265 rebAij->CreateView(
"stridedMaps", Aij);
267 if (!rebAij.is_null()) {
268 RCP<ParameterList> params = rcp(
new ParameterList());
269 params->set(
"printLoadBalancingInfo",
true);
270 std::stringstream ss4; ss4 <<
"A(" << curBlockId <<
"," << j <<
") rebalanced:";
282 subBlockRebA[curBlockId*bA->Cols() + j] = rebAij;
286 for(
size_t i=0; i<bA->Rows(); i++) {
287 if(i==curBlockId)
continue;
290 Teuchos::RCP<Matrix> Aij = subBlockRebA[i*bA->Cols() + curBlockId];
292 Teuchos::RCP<Matrix> rebAij;
293 if(rebalanceImporter!=Teuchos::null) {
294 std::stringstream ss; ss <<
"Rebalancing matrix block (" << i <<
"," << curBlockId <<
")";
296 RCP<const Map> targetMap = rebalanceImporter->getTargetMap();
299 Teuchos::RCP<const CrsMatrixWrap> rebAwij = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(rebAij);
300 Teuchos::RCP<CrsMatrix> rebAmij = rebAwij->getCrsMatrix();
301 Teuchos::RCP<const Import> rebAijImport = ImportFactory::Build(targetMap,Aij->getColMap());
302 rebAmij->replaceDomainMapAndImporter(targetMap,rebAijImport);
305 if (!rebAij.is_null() && Aij->IsView(
"stridedMaps"))
306 rebAij->CreateView(
"stridedMaps", Aij);
308 if (!rebAij.is_null()) {
309 RCP<ParameterList> params = rcp(
new ParameterList());
310 params->set(
"printLoadBalancingInfo",
true);
311 std::stringstream ss2; ss2 <<
"A(" << i <<
"," << curBlockId <<
") rebalanced:";
323 subBlockRebA[i*bA->Cols() + curBlockId] = rebAij;
328 subBlockARangeMaps.push_back(rebAii->getRowMap(
"stridedMaps"));
329 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMap = rebAii->getRangeMap()->getNodeElementList();
330 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
331 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
333 subBlockADomainMaps.push_back(rebAii->getColMap(
"stridedMaps"));
334 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMap = rebAii->getDomainMap()->getNodeElementList();
335 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
336 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
344 GO rangeIndexBase = 0;
345 GO domainIndexBase = 0;
346 rangeIndexBase = bA->getRangeMap()->getIndexBase();
347 domainIndexBase= bA->getDomainMap()->getIndexBase();
349 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
350 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getFullMap());
351 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
352 if(stridedRgFullMap != Teuchos::null) {
353 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
355 StridedMapFactory::Build(
356 bA->getRangeMap()->lib(),
357 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
361 bA->getRangeMap()->getComm(),
362 stridedRgFullMap->getStridedBlockId(),
363 stridedRgFullMap->getOffset());
367 bA->getRangeMap()->lib(),
368 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
371 bA->getRangeMap()->getComm());
374 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
376 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getFullMap());
377 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
378 if(stridedDoFullMap != Teuchos::null) {
379 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
380 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
382 StridedMapFactory::Build(
383 bA->getDomainMap()->lib(),
384 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
388 bA->getDomainMap()->getComm(),
389 stridedDoFullMap->getStridedBlockId(),
390 stridedDoFullMap->getOffset());
395 bA->getDomainMap()->lib(),
396 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
399 bA->getDomainMap()->getComm());
403 Teuchos::RCP<const MapExtractorClass> rebRangeMapExtractor = MapExtractorFactoryClass::Build(fullRangeMap, subBlockARangeMaps);
404 Teuchos::RCP<const MapExtractorClass> rebDomainMapExtractor = MapExtractorFactoryClass::Build(fullDomainMap, subBlockADomainMaps);
406 Teuchos::RCP<BlockedCrsMatrix> reb_bA = Teuchos::rcp(
new BlockedCrsMatrix(rebRangeMapExtractor,rebDomainMapExtractor,10));
407 for(
size_t i=0; i<bA->Rows(); i++) {
408 for(
size_t j=0; j<bA->Cols(); j++) {
409 Teuchos::RCP<const CrsMatrixWrap> crsOpij = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(subBlockRebA[i*bA->Cols() + j]);
410 Teuchos::RCP<CrsMatrix> crsMatij = crsOpij->getCrsMatrix();
411 reb_bA->setMatrix(i,j,crsMatij);
414 reb_bA->fillComplete();
416 coarseLevel.Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(reb_bA),
this);
422 if (rebalanceFacts_.begin() != rebalanceFacts_.end()) {
426 for (std::vector<RCP<const FactoryBase> >::const_iterator it2 = rebalanceFacts_.begin(); it2 != rebalanceFacts_.end(); ++it2) {
427 GetOStream(
Runtime0) <<
"RebalanceBlockedAc: call rebalance factory " << (*it2).get() <<
": " << (*it2)->description() << std::endl;
428 (*it2)->CallBuild(coarseLevel);
433 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
440 rebalanceFacts_.push_back(factory);
Exception indicating invalid cast attempted.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
Timer to be used in factories. Similar to Monitor but with additional timers.
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 DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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()