111 RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel,
"A");
113 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
114 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
116 const int numFactManagers = FactManager_.size();
120 "Number of block rows [" << A->Rows() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
122 "Number of block cols [" << A->Cols() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
125 std::vector<RCP<Matrix> > subBlockP (numFactManagers);
126 std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
127 std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
129 std::vector<GO> fullRangeMapVector;
130 std::vector<GO> fullDomainMapVector;
132 RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
133 RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
135 const ParameterList& pL = GetParameterList();
136 const bool backwards = pL.get<
bool>(
"backwards");
141 for (
int k = 0; k < numFactManagers; k++) {
142 int i = (backwards ? numFactManagers-1 - k : k);
143 if (restrictionMode_) GetOStream(
Runtime1) <<
"Generating R for block " << k <<
"/"<<numFactManagers <<std::endl;
144 else GetOStream(
Runtime1) <<
"Generating P for block " << k <<
"/"<<numFactManagers <<std::endl;
146 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
151 if (!restrictionMode_) subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"P", factManager->GetFactory(
"P").get());
152 else subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"R", factManager->GetFactory(
"R").get());
155 TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView(
"stridedMaps") ==
false,
Exceptions::BadCast,
156 "subBlock P operator [" << i <<
"] has no strided map information.");
161 RCP<const StridedMap> strPartialMap = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getRowMap(
"stridedMaps"));
162 std::vector<size_t> stridedRgData = strPartialMap->getStridingData();
163 subBlockPRangeMaps[i] = StridedMapFactory::Build(
164 subBlockP[i]->getRangeMap(),
166 strPartialMap->getStridedBlockId(),
167 strPartialMap->getOffset());
171 ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getLocalElementList();
172 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
175 RCP<const StridedMap> strPartialMap2 = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getColMap(
"stridedMaps"));
176 std::vector<size_t> stridedRgData2 = strPartialMap2->getStridingData();
177 subBlockPDomainMaps[i] = StridedMapFactory::Build(
178 subBlockP[i]->getDomainMap(),
180 strPartialMap2->getStridedBlockId(),
181 strPartialMap2->getOffset());
185 ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getLocalElementList();
186 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
196 bool bDoSortRangeGIDs = rangeAMapExtractor->getThyraMode() ? false :
true;
197 bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false :
true;
199 if (bDoSortRangeGIDs) {
200 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
201 fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
203 if (bDoSortDomainGIDs) {
204 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
205 fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
209 GO rangeIndexBase = 0;
210 GO domainIndexBase = 0;
211 if (!restrictionMode_) {
213 rangeIndexBase = A->getRangeMap() ->getIndexBase();
214 domainIndexBase = A->getDomainMap()->getIndexBase();
218 rangeIndexBase = A->getDomainMap()->getIndexBase();
219 domainIndexBase = A->getRangeMap()->getIndexBase();
225 RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<const StridedMap>(rangeAMapExtractor->getFullMap());
226 RCP<const Map > fullRangeMap = Teuchos::null;
228 ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
229 if (stridedRgFullMap != Teuchos::null) {
230 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
231 fullRangeMap = StridedMapFactory::Build(
232 A->getRangeMap()->lib(),
233 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
237 A->getRangeMap()->getComm(),
239 stridedRgFullMap->getOffset());
241 fullRangeMap = MapFactory::Build(
242 A->getRangeMap()->lib(),
243 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
246 A->getRangeMap()->getComm());
249 ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
250 RCP<const StridedMap> stridedDoFullMap = rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
251 RCP<const Map > fullDomainMap = null;
252 if (stridedDoFullMap != null) {
254 "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information!");
256 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
257 fullDomainMap = StridedMapFactory::Build(
258 A->getDomainMap()->lib(),
259 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
263 A->getDomainMap()->getComm(),
265 stridedDoFullMap->getOffset());
267 fullDomainMap = MapFactory::Build(
268 A->getDomainMap()->lib(),
269 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
272 A->getDomainMap()->getComm());
276 RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
277 RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());
279 RCP<BlockedCrsMatrix> P = rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
280 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
281 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
283 RCP<CrsMatrixWrap> crsOpii = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
284 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null, Xpetra::Exceptions::BadCast,
285 "Block [" << i <<
","<< j <<
"] must be of type CrsMatrixWrap.");
286 P->setMatrix(i, i, crsOpii);
288 P->setMatrix(i, j, Teuchos::null);
294 if (!restrictionMode_) {
296 coarseLevel.
Set(
"P", rcp_dynamic_cast<Matrix>(P),
this);
298 coarseLevel.
Set(
"CoarseMap",P->getBlockedDomainMap(),
this);
303 coarseLevel.
Set(
"R", Teuchos::rcp_dynamic_cast<Matrix>(P),
this);