46 #ifndef MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKINTERPOLATIONFACTORY_DEF_HPP_ 49 #ifdef HAVE_MUELU_EXPERIMENTAL 51 #include <Teuchos_Tuple.hpp> 53 #include "Xpetra_Vector.hpp" 54 #include "Xpetra_VectorFactory.hpp" 55 #include "Xpetra_MultiVector.hpp" 56 #include "Xpetra_MultiVectorFactory.hpp" 57 #include <Xpetra_Matrix.hpp> 58 #include <Xpetra_BlockedCrsMatrix.hpp> 59 #include <Xpetra_MapFactory.hpp> 60 #include <Xpetra_MapExtractor.hpp> 61 #include <Xpetra_MapExtractorFactory.hpp> 62 #include <Xpetra_MatrixFactory.hpp> 63 #include <Xpetra_Import.hpp> 64 #include <Xpetra_ImportFactory.hpp> 69 #include "MueLu_HierarchyHelpers.hpp" 72 #include "MueLu_PerfUtils.hpp" 76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
80 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Factory of the prolongation operator that need to be rebalanced (only used if type=Interpolation)");
85 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 87 #undef SET_VALID_ENTRY 95 return validParamList;
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 FactManager_.push_back(FactManager);
103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 Input(coarseLevel,
"P");
108 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
109 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
112 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
122 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
123 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"P");
125 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
126 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
127 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
130 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Rows() != 2,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockTransferFactory::Build: number of block rows of transfer operator is not equal 2. error.");
131 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Cols() != 2,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockTransferFactory::Build: number of block columns of transfer operator is not equal 2. error.");
134 std::vector<GO> fullRangeMapVector;
135 std::vector<GO> fullDomainMapVector;
136 std::vector<RCP<const Map> > subBlockPRangeMaps;
137 std::vector<RCP<const Map> > subBlockPDomainMaps;
138 subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows());
139 subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols());
141 std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
142 subBlockRebP.reserve(bOriginalTransferOp->Rows());
145 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
146 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
147 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
152 rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
155 Teuchos::RCP<CrsMatrix> Pmii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
156 Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp(
new CrsMatrixWrap(Pmii));
157 Teuchos::RCP<Matrix> Pii = Teuchos::rcp_dynamic_cast<Matrix>(Pwii);
160 if(rebalanceImporter != Teuchos::null) {
161 std::stringstream ss; ss <<
"Rebalancing prolongator block P(" << curBlockId <<
"," << curBlockId <<
")";
176 RCP<const Import> newImporter;
178 SubFactoryMonitor subM(*
this,
"Rebalancing prolongator -- fast map replacement", coarseLevel);
179 newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
180 Pmii->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
183 RCP<ParameterList> params = rcp(
new ParameterList());
184 params->set(
"printLoadBalancingInfo",
true);
185 std::stringstream ss2; ss2 <<
"P(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
189 subBlockRebP.push_back(Pii);
192 RCP<ParameterList> params = rcp(
new ParameterList());
193 params->set(
"printLoadBalancingInfo",
true);
194 std::stringstream ss; ss <<
"P(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
197 subBlockRebP.push_back(Pii);
201 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgPMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
202 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rgPMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
203 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
204 if(orig_stridedRgMap != Teuchos::null) {
205 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
206 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
207 stridedRgMap = StridedMapFactory::Build(
208 originalTransferOp->getRangeMap()->lib(),
209 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
211 Pii->getRangeMap()->getIndexBase(),
213 originalTransferOp->getRangeMap()->getComm(),
214 orig_stridedRgMap->getStridedBlockId(),
215 orig_stridedRgMap->getOffset());
216 }
else stridedRgMap = Pii->getRangeMap();
217 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doPMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
218 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(doPMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
219 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
220 if(orig_stridedDoMap != Teuchos::null) {
221 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
222 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
223 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
224 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
226 Pii->getDomainMap()->getIndexBase(),
228 originalTransferOp->getDomainMap()->getComm(),
229 orig_stridedDoMap->getStridedBlockId(),
230 orig_stridedDoMap->getOffset());
231 }
else stridedDoMap = Pii->getDomainMap();
233 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
234 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
237 if(Pii->IsView(
"stridedMaps")) Pii->RemoveView(
"stridedMaps");
238 Pii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
241 Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap(
"stridedMaps");
242 subBlockPRangeMaps.push_back(rangeMapii);
243 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getNodeElementList();
244 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
245 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
248 Teuchos::RCP<const Map> domainMapii = Pii->getColMap(
"stridedMaps");
249 subBlockPDomainMaps.push_back(domainMapii);
250 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getNodeElementList();
251 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
252 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
259 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
260 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
262 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
263 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
264 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangePMapExtractor->getFullMap());
265 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
266 if(stridedRgFullMap != Teuchos::null) {
267 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
269 StridedMapFactory::Build(
270 originalTransferOp->getRangeMap()->lib(),
271 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
275 originalTransferOp->getRangeMap()->getComm(),
276 stridedRgFullMap->getStridedBlockId(),
277 stridedRgFullMap->getOffset());
281 originalTransferOp->getRangeMap()->lib(),
282 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
285 originalTransferOp->getRangeMap()->getComm());
288 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
289 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
290 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
291 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
292 if(stridedDoFullMap != Teuchos::null) {
293 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
294 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
296 StridedMapFactory::Build(
297 originalTransferOp->getDomainMap()->lib(),
298 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
302 originalTransferOp->getDomainMap()->getComm(),
303 stridedDoFullMap->getStridedBlockId(),
304 stridedDoFullMap->getOffset());
309 originalTransferOp->getDomainMap()->lib(),
310 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
313 originalTransferOp->getDomainMap()->getComm());
317 Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangeMapExtractor =
318 Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullRangeMap, subBlockPRangeMaps);
319 Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainMapExtractor =
320 Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullDomainMap, subBlockPDomainMaps);
322 Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor,domainMapExtractor,10));
323 for(
size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
324 Teuchos::RCP<const CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(subBlockRebP[i]);
325 Teuchos::RCP<CrsMatrix> crsMatii = crsOpii->getCrsMatrix();
326 bRebP->setMatrix(i,i,crsMatii);
328 bRebP->fillComplete();
330 Set(coarseLevel,
"P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
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
Specifies the data that this class needs, and the factories that generate that data.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
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 AddFactoryManager(RCP< const FactoryManagerBase > FactManager)
Add a factory manager.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()