131 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LO, GO, NO> xdMV;
132 typedef Xpetra::BlockedMultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> xdBV;
134 bool UseSingleSource = FactManager_.size() == 0;
137 Teuchos::RCP<Matrix> nonrebCoarseA = Get< RCP<Matrix> >(coarseLevel,
"A");
139 Teuchos::RCP<Matrix> originalTransferOp = Teuchos::null;
140 originalTransferOp = Get< RCP<Matrix> >(coarseLevel,
"P");
142 RCP<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bOriginalTransferOp =
143 Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > (originalTransferOp);
144 TEUCHOS_TEST_FOR_EXCEPTION(bOriginalTransferOp==Teuchos::null,
Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: input matrix P or R is not of type BlockedCrsMatrix! error.");
146 RCP<const MapExtractor> rangeMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
147 RCP<const MapExtractor> domainMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
156 bool bThyraRangeGIDs = rangeMapExtractor->getThyraMode();
157 bool bThyraDomainGIDs = domainMapExtractor->getThyraMode();
160 std::vector<GO> fullRangeMapVector;
161 std::vector<GO> fullDomainMapVector;
162 std::vector<RCP<const Map> > subBlockPRangeMaps;
163 std::vector<RCP<const Map> > subBlockPDomainMaps;
164 subBlockPRangeMaps.reserve(bOriginalTransferOp->Rows());
165 subBlockPDomainMaps.reserve(bOriginalTransferOp->Cols());
167 std::vector<Teuchos::RCP<Matrix> > subBlockRebP;
168 subBlockRebP.reserve(bOriginalTransferOp->Rows());
171 std::vector<RCP<const Import> > importers = std::vector<RCP<const Import> >(bOriginalTransferOp->Rows(), Teuchos::null);
172 std::vector<RCP<xdMV> > newCoordinates(bOriginalTransferOp->Rows());
173 RCP<xdBV> oldCoordinates;
174 if(UseSingleSource) {
175 importers = Get<std::vector<RCP<const Import> > >(coarseLevel,
"SubImporters");
176 oldCoordinates = Get<RCP<xdBV> >(coarseLevel,
"Coordinates");
180 Teuchos::RCP<const Import> rebalanceImporter = Teuchos::null;
181 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
182 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
188 if(UseSingleSource) rebalanceImporter = importers[curBlockId];
189 else rebalanceImporter = coarseLevel.
Get<Teuchos::RCP<const Import> >(
"Importer", (*it)->GetFactory(
"Importer").get());
192 Teuchos::RCP<Matrix> Pii = bOriginalTransferOp->getMatrix(curBlockId, curBlockId);
193 Teuchos::RCP<CrsMatrixWrap> Pwii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Pii);
194 TEUCHOS_TEST_FOR_EXCEPTION(Pwii == Teuchos::null,Xpetra::Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: block " << curBlockId <<
" is not of type CrsMatrixWrap. We need an underlying CsrMatrix to replace domain map and importer!");
198 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block range " << curBlockId <<
" start with " << Pii->getRowMap()->getMinAllGlobalIndex() <<
" but should start with 0");
201 "MueLu::RebalanceBlockInterpolationFactory::Build: inconsistent Thyra GIDs. Thyra global ids for block domain " << curBlockId <<
" start with " << Pii->getColMap()->getMinAllGlobalIndex() <<
" but should start with 0");
204 if(rebalanceImporter != Teuchos::null) {
205 std::stringstream ss; ss <<
"Rebalancing prolongator block P(" << curBlockId <<
"," << curBlockId <<
")";
219 RCP<const Import> newImporter;
221 SubFactoryMonitor subM(*
this,
"Rebalancing prolongator -- fast map replacement", coarseLevel);
222 newImporter = ImportFactory::Build(rebalanceImporter->getTargetMap(), Pii->getColMap());
223 Pwii->getCrsMatrix()->replaceDomainMapAndImporter(rebalanceImporter->getTargetMap(), newImporter);
226 RCP<ParameterList> params = rcp(
new ParameterList());
227 params->set(
"printLoadBalancingInfo",
true);
228 std::stringstream ss2; ss2 <<
"P(" << curBlockId <<
"," << curBlockId <<
") rebalanced:";
232 subBlockRebP.push_back(Pii);
237 if(UseSingleSource) {
238 RCP<xdMV> localCoords = oldCoordinates->getMultiVector(curBlockId);
241 RCP<const Import> coordImporter = rebalanceImporter;
243 RCP<xdMV> permutedLocalCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), localCoords->getNumVectors());
244 permutedLocalCoords->doImport(*localCoords, *coordImporter, Xpetra::INSERT);
246 newCoordinates[curBlockId] = permutedLocalCoords;
248 else if ( (*it)->hasFactory(
"Coordinates") ==
true && coarseLevel.
IsAvailable(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()) ==
true) {
249 RCP<xdMV> coords = coarseLevel.
Get< RCP<xdMV> >(
"Coordinates",(*it)->GetFactory(
"Coordinates").get());
254 LO nodeNumElts = coords->getMap()->getLocalNumElements();
257 LO myBlkSize = 0, blkSize = 0;
259 if (nodeNumElts > 0) {
262 "MueLu::RebalanceBlockInterpolationFactory::Build: block size. " << rebalanceImporter->getSourceMap()->getLocalNumElements() <<
" not divisable by " << nodeNumElts);
263 myBlkSize = rebalanceImporter->getSourceMap()->getLocalNumElements() / nodeNumElts;
266 MueLu_maxAll(coords->getMap()->getComm(), myBlkSize, blkSize);
268 RCP<const Import> coordImporter = Teuchos::null;
270 coordImporter = rebalanceImporter;
275 RCP<const Map> origMap = coords->getMap();
276 GO indexBase = origMap->getIndexBase();
278 ArrayView<const GO> OEntries = rebalanceImporter->getTargetMap()->getLocalElementList();
279 LO numEntries = OEntries.size()/blkSize;
280 ArrayRCP<GO> Entries(numEntries);
281 for (LO i = 0; i < numEntries; i++)
282 Entries[i] = (OEntries[i*blkSize]-indexBase)/blkSize + indexBase;
284 RCP<const Map> targetMap = MapFactory::Build(origMap->lib(), origMap->getGlobalNumElements(), Entries(), indexBase, origMap->getComm());
285 coordImporter = ImportFactory::Build(origMap, targetMap);
288 RCP<xdMV> permutedCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coordImporter->getTargetMap(), coords->getNumVectors());
289 permutedCoords->doImport(*coords, *coordImporter, Xpetra::INSERT);
291 const ParameterList& pL = GetParameterList();
292 if (pL.isParameter(
"repartition: use subcommunicators") ==
true && pL.get<
bool>(
"repartition: use subcommunicators") ==
true)
293 permutedCoords->replaceMap(permutedCoords->getMap()->removeEmptyProcesses());
295 Set(coarseLevel,
"Coordinates", permutedCoords);
299 RCP<ParameterList> params = rcp(
new ParameterList());
300 params->set(
"printLoadBalancingInfo",
true);
301 std::stringstream ss; ss <<
"P(" << curBlockId <<
"," << curBlockId <<
") not rebalanced:";
304 subBlockRebP.push_back(Pii);
309 if((*it)->hasFactory(
"Coordinates") ==
true && coarseLevel.
IsAvailable(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()) ==
true) {
310 coarseLevel.
Set(
"Coordinates", coarseLevel.
Get< RCP<xdMV> >(
"Coordinates",(*it)->GetFactory(
"Coordinates").get()),
this);
316 Teuchos::RCP<const StridedMap> orig_stridedRgMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangeMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),rangeMapExtractor->getThyraMode()));
317 Teuchos::RCP<const Map> stridedRgMap = Teuchos::null;
318 if(orig_stridedRgMap != Teuchos::null) {
319 std::vector<size_t> stridingData = orig_stridedRgMap->getStridingData();
320 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
321 stridedRgMap = StridedMapFactory::Build(
322 originalTransferOp->getRangeMap()->lib(),
323 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
325 Pii->getRangeMap()->getIndexBase(),
327 originalTransferOp->getRangeMap()->getComm(),
328 orig_stridedRgMap->getStridedBlockId(),
329 orig_stridedRgMap->getOffset());
330 }
else stridedRgMap = Pii->getRangeMap();
333 Teuchos::RCP<const StridedMap> orig_stridedDoMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainMapExtractor->getMap(Teuchos::as<size_t>(curBlockId),domainMapExtractor->getThyraMode()));
335 Teuchos::RCP<const Map> stridedDoMap = Teuchos::null;
336 if(orig_stridedDoMap != Teuchos::null) {
337 std::vector<size_t> stridingData = orig_stridedDoMap->getStridingData();
338 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
339 stridedDoMap = StridedMapFactory::Build(originalTransferOp->getDomainMap()->lib(),
340 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
342 Pii->getDomainMap()->getIndexBase(),
344 originalTransferOp->getDomainMap()->getComm(),
345 orig_stridedDoMap->getStridedBlockId(),
346 orig_stridedDoMap->getOffset());
347 }
else stridedDoMap = Pii->getDomainMap();
349 TEUCHOS_TEST_FOR_EXCEPTION(stridedRgMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
350 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoMap == Teuchos::null,
Exceptions::RuntimeError,
"MueLu::RebalanceBlockInterpolationFactory::Build: failed to generate striding information. error.");
353 if(Pii->IsView(
"stridedMaps")) Pii->RemoveView(
"stridedMaps");
354 Pii->CreateView(
"stridedMaps", stridedRgMap, stridedDoMap);
357 Teuchos::RCP<const Map> rangeMapii = Pii->getRowMap(
"stridedMaps");
358 subBlockPRangeMaps.push_back(rangeMapii);
359 Teuchos::ArrayView< const GlobalOrdinal > nodeRangeMapii = Pii->getRangeMap()->getLocalElementList();
361 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMapii.begin(), nodeRangeMapii.end());
362 if(bThyraRangeGIDs ==
false)
363 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
366 Teuchos::RCP<const Map> domainMapii = Pii->getColMap(
"stridedMaps");
367 subBlockPDomainMaps.push_back(domainMapii);
368 Teuchos::ArrayView< const GlobalOrdinal > nodeDomainMapii = Pii->getDomainMap()->getLocalElementList();
370 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMapii.begin(), nodeDomainMapii.end());
371 if(bThyraDomainGIDs ==
false)
372 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
379 GO rangeIndexBase = originalTransferOp->getRangeMap()->getIndexBase();
380 GO domainIndexBase= originalTransferOp->getDomainMap()->getIndexBase();
382 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > rangePMapExtractor = bOriginalTransferOp->getRangeMapExtractor();
383 Teuchos::ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0,fullRangeMapVector.size());
384 Teuchos::RCP<const StridedMap> stridedRgFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(rangePMapExtractor->getFullMap());
385 Teuchos::RCP<const Map > fullRangeMap = Teuchos::null;
386 if(stridedRgFullMap != Teuchos::null) {
387 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
389 StridedMapFactory::Build(
390 originalTransferOp->getRangeMap()->lib(),
391 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
395 originalTransferOp->getRangeMap()->getComm(),
396 stridedRgFullMap->getStridedBlockId(),
397 stridedRgFullMap->getOffset());
401 originalTransferOp->getRangeMap()->lib(),
402 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
405 originalTransferOp->getRangeMap()->getComm());
408 RCP<const Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> > domainAMapExtractor = bOriginalTransferOp->getDomainMapExtractor();
409 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
410 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
411 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
412 if(stridedDoFullMap != Teuchos::null) {
413 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
414 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
416 StridedMapFactory::Build(
417 originalTransferOp->getDomainMap()->lib(),
418 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
422 originalTransferOp->getDomainMap()->getComm(),
423 stridedDoFullMap->getStridedBlockId(),
424 stridedDoFullMap->getOffset());
429 originalTransferOp->getDomainMap()->lib(),
430 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
433 originalTransferOp->getDomainMap()->getComm());
437 Teuchos::RCP<const MapExtractor> rebrangeMapExtractor =
438 MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bThyraRangeGIDs);
439 Teuchos::RCP<const MapExtractor> rebdomainMapExtractor =
440 MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bThyraDomainGIDs);
442 Teuchos::RCP<BlockedCrsMatrix> bRebP = Teuchos::rcp(
new BlockedCrsMatrix(rebrangeMapExtractor,rebdomainMapExtractor,10));
443 for(
size_t i = 0; i<subBlockPRangeMaps.size(); i++) {
444 Teuchos::RCP<CrsMatrixWrap> crsOpii = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(subBlockRebP[i]);
445 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null,Xpetra::Exceptions::BadCast,
"MueLu::RebalanceBlockTransferFactory::Build: block P" << i <<
" is not of type CrsMatrixWrap.");
446 bRebP->setMatrix(i,i,crsOpii);
448 bRebP->fillComplete();
450 Set(coarseLevel,
"P", Teuchos::rcp_dynamic_cast<Matrix>(bRebP));
454 if(UseSingleSource) {
455 RCP<xdBV> bcoarseCoords = rcp(
new xdBV(rebrangeMapExtractor->getBlockedMap(),newCoordinates));
456 Set(coarseLevel,
"Coordinates",bcoarseCoords);