125 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
127 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
128 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"R");
130 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
131 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(originalTransferOp);
132 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
134 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
135 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
138 bool bRestrictComm =
false;
139 const ParameterList& pL = GetParameterList();
140 if (pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
141 bRestrictComm =
true;
150 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
151 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
154 std::vector<GO> fullRangeMapVector;
155 std::vector<GO> fullDomainMapVector;
156 std::vector<RCP<const Map> > subBlockRRangeMaps;
157 std::vector<RCP<const Map> > subBlockRDomainMaps;
158 subBlockRRangeMaps.reserve(bOriginalTransferOp->Rows());
159 subBlockRDomainMaps.reserve(bOriginalTransferOp->Cols());
161 std::vector<Teuchos::RCP<Matrix> > subBlockRebR;
162 subBlockRebR.reserve(bOriginalTransferOp->Cols());
164 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
165 if(UseSingleSourceImporters_) {
166 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
170 Teuchos::RCP<const Import> rebalanceImporter;
171 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
172 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
177 if(UseSingleSourceImporters_) rebalanceImporter = importers[curBlockId];
178 else rebalanceImporter = coarseLevel.
Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
181 Teuchos::RCP<Matrix> Rii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
184 TEUCHOS_TEST_FOR_EXCEPTION(bThyraRangeGIDs ==
true && Rii->getRowMap()->getMinAllGlobalIndex() != 0,
186 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Rii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
187 TEUCHOS_TEST_FOR_EXCEPTION(bThyraDomainGIDs ==
true && Rii->getColMap()->getMinAllGlobalIndex() != 0,
189 "MueLu::RebalanceBlockRestrictionFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Rii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
191 Teuchos::RCP<Matrix> rebRii;
192 if(rebalanceImporter != Teuchos::null) {
193 std::stringstream ss; ss <<
"Rebalancing restriction block R(" << curBlockId <<
"," << curBlockId <<
")";
196 SubFactoryMonitor subM(*
this,
"Rebalancing restriction -- fusedImport", coarseLevel);
200 rebRii = MatrixFactory::Build(Rii,*rebalanceImporter,dummy,rebalanceImporter->getTargetMap());
203 RCP<ParameterList> params = rcp(
new ParameterList());
204 params->set(
"printLoadBalancingInfo",
true);
205 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
209 RCP<ParameterList> params = rcp(
new ParameterList());
210 params->set(
"printLoadBalancingInfo",
true);
211 std::stringstream ss2; ss2 <<
"R(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
216 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
217 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
218 if(orig_stridedRgMap != Teuchos::null) {
219 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
220 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getLocalElementList();
221 stridedRgMap = StridedMapFactory::Build(
222 originalTransferOp->getRangeMap()->lib(),
223 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
225 rebRii->getRangeMap()->getIndexBase(),
227 originalTransferOp->getRangeMap()->getComm(),
228 orig_stridedRgMap->getStridedBlockId(),
229 orig_stridedRgMap->getOffset());
230 }
else stridedRgMap = Rii->getRangeMap();
232 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
233 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
234 if(orig_stridedDoMap != Teuchos::null) {
235 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
236 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getLocalElementList();
237 stridedDoMap = StridedMapFactory::Build(
238 originalTransferOp->getDomainMap()->lib(),
239 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
241 rebRii->getDomainMap()->getIndexBase(),
243 originalTransferOp->getDomainMap()->getComm(),
244 orig_stridedDoMap->getStridedBlockId(),
245 orig_stridedDoMap->getOffset());
246 }
else stridedDoMap = Rii->getDomainMap();
249 stridedRgMap->removeEmptyProcesses();
250 stridedDoMap->removeEmptyProcesses();
253 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
254 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockRestrictionFactory::Build: failed to generate striding information. error.");
257 if(rebRii->IsView(
"stridedMaps")) rebRii->RemoveView(
"stridedMaps");
258 rebRii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
261 subBlockRebR.push_back(rebRii);
264 Teuchos::RCP<const Map> rangeMapii = rebRii->getRowMap(
"stridedMaps");
265 subBlockRRangeMaps.push_back(rangeMapii);
266 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = rebRii->getRangeMap()->getLocalElementList();
268 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
269 if(bThyraRangeGIDs ==
false)
270 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
273 Teuchos::RCP<const Map> domainMapii = rebRii->getColMap(
"stridedMaps");
274 subBlockRDomainMaps.push_back(domainMapii);
275 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = rebRii->getDomainMap()->getLocalElementList();
277 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
278 if(bThyraDomainGIDs ==
false)
279 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
286 if(rebalanceImporter != Teuchos::null)
288 std::stringstream ss2; ss2 <<
"Rebalancing nullspace block(" << curBlockId <<
"," << curBlockId <<
")";
291 RCP<MultiVector> nullspace = coarseLevel.
Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
292 RCP<MultiVector> permutedNullspace = MultiVectorFactory::Build(rebalanceImporter->getTargetMap(), nullspace->getNumVectors());
293 permutedNullspace->doImport(*nullspace, *rebalanceImporter, Xpetra::INSERT);
297 permutedNullspace->replaceMap(permutedNullspace->getMap()->removeEmptyProcesses());
299 coarseLevel.
Set<RCP<MultiVector> >(
"Nullspace", permutedNullspace, (*it)->GetFactory(
"Nullspace").get());
303 RCP<MultiVector> nullspace = coarseLevel.
Get<RCP<MultiVector> >(
"Nullspace", (*it)->GetFactory(
"Nullspace").get());
304 coarseLevel.
Set<RCP<MultiVector> >(
"Nullspace", nullspace, (*it)->GetFactory(
"Nullspace").get());
313 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
314 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
317 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
318 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getFullMap());
319 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
320 if(stridedRgFullMap != Teuchos::null) {
321 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
323 StridedMapFactory::Build(
324 originalTransferOp->getRangeMap()->lib(),
325 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
329 originalTransferOp->getRangeMap()->getComm(),
330 stridedRgFullMap->getStridedBlockId(),
331 stridedRgFullMap->getOffset());
335 originalTransferOp->getRangeMap()->lib(),
336 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
339 originalTransferOp->getRangeMap()->getComm());
342 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
343 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getFullMap());
344 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
345 if(stridedDoFullMap != Teuchos::null) {
346 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
347 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
349 StridedMapFactory::Build(
350 originalTransferOp->getDomainMap()->lib(),
351 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
355 originalTransferOp->getDomainMap()->getComm(),
356 stridedDoFullMap->getStridedBlockId(),
357 stridedDoFullMap->getOffset());
362 originalTransferOp->getDomainMap()->lib(),
363 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
366 originalTransferOp->getDomainMap()->getComm());
370 fullRangeMap->removeEmptyProcesses();
371 fullDomainMap->removeEmptyProcesses();
375 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
376 MapExtractorFactory::Build(fullRangeMap, subBlockRRangeMaps, bThyraRangeGIDs);
377 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
378 MapExtractorFactory::Build(fullDomainMap, subBlockRDomainMaps, bThyraDomainGIDs);
380 Teuchos::RCP<BlockedCrsMatrix> bRebR = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
381 for(
size_t i = 0; i<subBlockRRangeMaps.size(); i++) {
382 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebR[i]);
383 bRebR->setMatrix(i,i,crsOpii);
386 bRebR->fillComplete();
388 Set(coarseLevel,
"R", Teuchos::rcp_dynamic_cast<Matrix>(bRebR));