47 #ifndef MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 48 #define MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 53 #include <Xpetra_BlockedCrsMatrix.hpp> 54 #include <Xpetra_MapExtractor.hpp> 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_StridedMapFactory.hpp> 63 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
65 RCP<ParameterList> validParamList = rcp(
new ParameterList());
69 validParamList->set<
int>(
"block row", 0,
"Block row of subblock matrix A");
70 validParamList->set<
int>(
"block col", 0,
"Block column of subblock matrix A");
72 validParamList->set<std::string >(
"Range map: Striding info",
"{}",
"Striding information for range map");
73 validParamList->set<
LocalOrdinal>(
"Range map: Strided block id", -1,
"Strided block id for range map");
74 validParamList->set<std::string >(
"Domain map: Striding info",
"{}",
"Striding information for domain map");
75 validParamList->set<
LocalOrdinal>(
"Domain map: Strided block id", -1,
"Strided block id for domain map");
77 return validParamList;
80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 Input(currentLevel,
"A");
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
222 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 const ParameterList & pL = GetParameterList();
228 stridedBlockId = pL.get<
LocalOrdinal>(
"Range map: Strided block id");
230 stridedBlockId = pL.get<
LocalOrdinal>(
"Domain map: Strided block id");
235 if(bRange ==
true) str = std::string(
"Range map: Striding info");
236 else str = std::string(
"Domain map: Striding info");
238 if(pL.isParameter(str)) {
239 std::string strStridingInfo = pL.get<std::string>(str);
240 if(strStridingInfo.empty() ==
false) {
241 Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
242 stridingInfo = Teuchos::createVector(arrayVal);
246 if(stridingInfo.size() > 0)
return true;
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
RCP< const ParameterList > GetValidParameterList() const override
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level ¤tLevel) const override
Build an object with this factory.
Namespace for MueLu classes and methods.
Print even more statistics.
Class that holds all level-specific information.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level ¤tLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
static const RCP< const NoFactory > getRCP()
Static Get() functions.