89 const ParameterList& pL = GetParameterList();
90 size_t row = Teuchos::as<size_t>(pL.get<
int>(
"block row"));
91 size_t col = Teuchos::as<size_t>(pL.get<
int>(
"block col"));
93 RCP<Matrix> Ain = currentLevel.
Get< RCP<Matrix> >(
"A",this->GetFactory(
"A").get());
94 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
96 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
97 TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(),
Exceptions::RuntimeError,
"row [" << row <<
"] > A.Rows() [" << A->Rows() <<
"].");
98 TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(),
Exceptions::RuntimeError,
"col [" << col <<
"] > A.Cols() [" << A->Cols() <<
"].");
101 RCP<Matrix> Op = A->getMatrix(row, col);
108 RCP<BlockedCrsMatrix> bOp = rcp_dynamic_cast<BlockedCrsMatrix>(Op);
109 if (bOp != Teuchos::null) {
111 if (bOp->Rows() == 1 && bOp->Cols() == 1) {
113 Op = bOp->getCrsMatrix();
114 TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null,
Exceptions::BadCast,
115 "SubBlockAFactory::Build: sub block A[" << row <<
"," << col <<
"] must be a single block CrsMatrixWrap object!");
120 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") is a " << bOp->Rows() <<
"x" << bOp->Cols() <<
" block matrix" << std::endl;
121 GetOStream(
Statistics2) <<
"with altogether " << bOp->getGlobalNumRows() <<
"x" << bOp->getGlobalNumCols() <<
" rows and columns." << std::endl;
122 currentLevel.
Set(
"A", Op,
this);
138 TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null,
Exceptions::BadCast,
139 "SubBlockAFactory::Build: sub block A[" << row <<
"," << col <<
"] is NOT a BlockedCrsMatrix, but also NOT a CrsMatrixWrap object. This cannot be.");
142 RCP<const StridedMap> stridedRangeMap = Teuchos::null;
143 RCP<const StridedMap> stridedDomainMap = Teuchos::null;
146 std::vector<size_t> rangeStridingInfo;
147 std::vector<size_t> domainStridingInfo;
150 bool bRangeUserSpecified = CheckForUserSpecifiedBlockInfo(
true, rangeStridingInfo, rangeStridedBlockId);
151 bool bDomainUserSpecified = CheckForUserSpecifiedBlockInfo(
false, domainStridingInfo, domainStridedBlockId);
153 "MueLu::SubBlockAFactory[" << row <<
"," << col <<
"]: the user has to specify either both domain and range map or none.");
156 RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
157 RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
159 RCP<const Map> rangeMap = rangeMapExtractor ->getMap(row);
160 RCP<const Map> domainMap = domainMapExtractor->getMap(col);
163 if(bRangeUserSpecified) stridedRangeMap = rcp(
new StridedMap(rangeMap, rangeStridingInfo, rangeMap->getIndexBase(), rangeStridedBlockId, 0));
164 else stridedRangeMap = rcp_dynamic_cast<const StridedMap>(rangeMap);
166 if(bDomainUserSpecified) stridedDomainMap = rcp(
new StridedMap(domainMap, domainStridingInfo, domainMap->getIndexBase(), domainStridedBlockId, 0));
167 else stridedDomainMap = rcp_dynamic_cast<const StridedMap>(domainMap);
172 if (stridedRangeMap.is_null()) {
173 RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
174 RCP<const StridedMap> stridedFullRangeMap = rcp_dynamic_cast<const StridedMap>(fullRangeMap);
175 TEUCHOS_TEST_FOR_EXCEPTION(stridedFullRangeMap.is_null(),
Exceptions::BadCast,
"Full rangeMap is not a strided map.");
177 std::vector<size_t> stridedData = stridedFullRangeMap->getStridingData();
178 if (stridedData.size() == 1 && row > 0) {
180 stridedRangeMap = StridedMapFactory::Build(rangeMap, stridedData, 0, stridedFullRangeMap->getOffset());
184 stridedRangeMap = StridedMapFactory::Build(rangeMap, stridedData, row, stridedFullRangeMap->getOffset());
188 if (stridedDomainMap.is_null()) {
189 RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
190 RCP<const StridedMap> stridedFullDomainMap = rcp_dynamic_cast<const StridedMap>(fullDomainMap);
191 TEUCHOS_TEST_FOR_EXCEPTION(stridedFullDomainMap.is_null(),
Exceptions::BadCast,
"Full domainMap is not a strided map");
193 std::vector<size_t> stridedData = stridedFullDomainMap->getStridingData();
194 if (stridedData.size() == 1 && col > 0) {
196 stridedDomainMap = StridedMapFactory::Build(domainMap, stridedData, 0, stridedFullDomainMap->getOffset());
200 stridedDomainMap = StridedMapFactory::Build(domainMap, stridedData, col, stridedFullDomainMap->getOffset());
204 TEUCHOS_TEST_FOR_EXCEPTION(stridedRangeMap.is_null(),
Exceptions::BadCast,
"rangeMap " << row <<
" is not a strided map.");
205 TEUCHOS_TEST_FOR_EXCEPTION(stridedDomainMap.is_null(),
Exceptions::BadCast,
"domainMap " << col <<
" is not a strided map.");
207 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") is a single block and has strided maps:"
208 <<
"\n range map fixed block size = " << stridedRangeMap ->getFixedBlockSize() <<
", strided block id = " << stridedRangeMap ->getStridedBlockId()
209 <<
"\n domain map fixed block size = " << stridedDomainMap->getFixedBlockSize() <<
", strided block id = " << stridedDomainMap->getStridedBlockId() << std::endl;
210 GetOStream(
Statistics2) <<
"A(" << row <<
"," << col <<
") has " << Op->getGlobalNumRows() <<
"x" << Op->getGlobalNumCols() <<
" rows and columns." << std::endl;
213 if (Op->IsView(
"stridedMaps") ==
true)
214 Op->RemoveView(
"stridedMaps");
215 Op->CreateView(
"stridedMaps", stridedRangeMap, stridedDomainMap);
217 TEUCHOS_TEST_FOR_EXCEPTION(Op->IsView(
"stridedMaps") ==
false,
Exceptions::RuntimeError,
"Failed to set \"stridedMaps\" view.");
219 currentLevel.
Set(
"A", Op,
this);