46 #ifndef MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_DEF_HPP_ 47 #define MUELU_REBALANCEBLOCKRESTRICTIONFACTORY_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> 68 #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> >(
"R", Teuchos::null,
"Factory of the restriction operator that need to be rebalanced (only used if type=Restriction)");
82 return validParamList;
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 FactManager_.push_back(FactManager);
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 Input(coarseLevel,
"R");
94 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
95 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
99 coarseLevel.
DeclareInput(
"Importer",(*it)->GetFactory(
"Importer").get(),
this);
100 coarseLevel.
DeclareInput(
"Nullspace",(*it)->GetFactory(
"Nullspace").get(),
this);
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
109 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
111 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
112 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"R");
114 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
115 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalTransferOp);
116 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
119 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Rows() != 2,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockTransferFactory::Build: number of block rows of transfer operator is not equal 2. error.");
120 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp->Cols() != 2,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockTransferFactory::Build: number of block columns of transfer operator is not equal 2. error.");
123 std::vector<GO> fullRangeMapVector;
124 std::vector<GO> fullDomainMapVector;
125 std::vector<RCP<const Map> > subBlockRRangeMaps;
126 std::vector<RCP<const Map> > subBlockRDomainMaps;
127 subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows());
128 subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols());
130 std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
131 subBlockRebR.reserve(bOriginalTransferOp->Cols());
134 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
135 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
136 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
141 rebalanceImporter = coarseLevel.Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
144 Teuchos::RCP<CrsMatrix> Rmii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
145 Teuchos::RCP<CrsMatrixWrap> Rwii = Teuchos::rcp(
new CrsMatrixWrap(Rmii));
146 Teuchos::RCP<Matrix> Rii = Teuchos::rcp_dynamic_cast<Matrix>(Rwii);
148 Teuchos::RCP<Matrix> rebRii;
149 if(rebalanceImporter != Teuchos::null) {
150 std::stringstream ss; ss <<
"Rebalancing restriction block R(" << curBlockId <<
"," << curBlockId <<
")";
153 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
157 rebRii = MatrixFactory::Build(Rii,*rebalanceImporter,dummy,rebalanceImporter->getTargetMap());
160 RCP<ParameterList> params = rcp(
new ParameterList());
161 params->set(
"printLoadBalancingInfo",
true);
162 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
166 RCP<ParameterList> params = rcp(
new ParameterList());
167 params->set(
"printLoadBalancingInfo",
true);
168 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
173 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rgRMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
174 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rgRMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
175 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
176 if(orig_stridedRgMap != Teuchos::null) {
177 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
178 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
179 stridedRgMap = StridedMapFactory::Build(
180 originalTransferOp->getRangeMap()->lib(),
181 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
183 rebRii->getRangeMap()->getIndexBase(),
185 originalTransferOp->getRangeMap()->getComm(),
186 orig_stridedRgMap->getStridedBlockId(),
187 orig_stridedRgMap->getOffset());
189 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > doRMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
190 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(doRMapExtractor->getMap(Teuchos::as<size_t>(curBlockId)));
191 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
192 if(orig_stridedDoMap != Teuchos::null) {
193 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
194 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
195 stridedDoMap = StridedMapFactory::Build(
196 originalTransferOp->getDomainMap()->lib(),
197 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
199 rebRii->getDomainMap()->getIndexBase(),
201 originalTransferOp->getDomainMap()->getComm(),
202 orig_stridedDoMap->getStridedBlockId(),
203 orig_stridedDoMap->getOffset());
206 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
207 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
210 if(rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
211 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
214 subBlockRebR.push_back(rebRii);
217 Teuchos::RCP<const Map> rangeMapii = rebRii->getRowMap(
"stridedMaps");
218 subBlockRRangeMaps.push_back(rangeMapii);
219 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getNodeElementList();
220 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
221 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
224 Teuchos::RCP<const Map> domainMapii = rebRii->getColMap(
"stridedMaps");
225 subBlockRDomainMaps.push_back(domainMapii);
226 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getNodeElementList();
227 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
228 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
233 if(rebalanceImporter != Teuchos::null)
235 std::stringstream ss2; ss2 <<
"Rebalancing nullspace block(" << curBlockId <<
"," << curBlockId <<
")";
238 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
239 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
240 permutedNullspace->doImport(*nullspace, *rebalanceImporter, Xpetra::INSERT);
246 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", permutedNullspace, (*it)->GetFactory(
"Nullspace").get());
250 RCP<MultiVector> nullspace = coarseLevel.Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
251 coarseLevel.Set<RCP<MultiVector> >(
"Nullspace", nullspace, (*it)->GetFactory(
"Nullspace").get());
260 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
261 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
264 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangeRMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
265 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
266 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(rangeRMapExtractor->getFullMap());
267 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
268 if(stridedRgFullMap != Teuchos::null) {
269 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
271 StridedMapFactory::Build(
272 originalTransferOp->getRangeMap()->lib(),
273 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
277 originalTransferOp->getRangeMap()->getComm(),
278 stridedRgFullMap->getStridedBlockId(),
279 stridedRgFullMap->getOffset());
283 originalTransferOp->getRangeMap()->lib(),
284 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
287 originalTransferOp->getRangeMap()->getComm());
290 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
291 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
292 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
293 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
294 if(stridedDoFullMap != Teuchos::null) {
295 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
296 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
298 StridedMapFactory::Build(
299 originalTransferOp->getDomainMap()->lib(),
300 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
304 originalTransferOp->getDomainMap()->getComm(),
305 stridedDoFullMap->getStridedBlockId(),
306 stridedDoFullMap->getOffset());
311 originalTransferOp->getDomainMap()->lib(),
312 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
315 originalTransferOp->getDomainMap()->getComm());
319 Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangeMapExtractor =
320 Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullRangeMap, subBlockRRangeMaps);
321 Teuchos::RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainMapExtractor =
322 Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(fullDomainMap, subBlockRDomainMaps);
324 Teuchos::RCP<BlockedCrsMatrix> bRebR = Teuchos::rcp(
new BlockedCrsMatrix(rangeMapExtractor,domainMapExtractor,10));
325 for(
size_t i = 0; i<subBlockRRangeMaps.size(); i++) {
326 Teuchos::RCP<const CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<
const CrsMatrixWrap>(subBlockRebR[i]);
327 Teuchos::RCP<CrsMatrix> crsMatii = crsOpii->getCrsMatrix();
328 bRebR->setMatrix(i,i,crsMatii);
331 bRebR->fillComplete();
333 Set(coarseLevel,
"R", Teuchos::rcp_dynamic_cast<Matrix>(bRebR));
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Exception indicating invalid cast attempted.
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)
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.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
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()