53 #ifndef MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 54 #define MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 59 #include <Xpetra_BlockedCrsMatrix.hpp> 60 #include <Xpetra_CrsMatrix.hpp> 61 #include <Xpetra_CrsMatrixWrap.hpp> 62 #include <Xpetra_MapExtractor.hpp> 63 #include <Xpetra_Matrix.hpp> 64 #include <Xpetra_StridedMapFactory.hpp> 71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 RCP<ParameterList> validParamList = rcp(
new ParameterList());
76 validParamList->set<
int > (
"block row", 0,
"Block row of subblock matrix A");
77 validParamList->set<
int > (
"block col", 0,
"Block column of subblock matrix A");
79 validParamList->set< std::string > (
"Range map: Striding info",
"{}",
"Striding information for range map");
80 validParamList->set< LocalOrdinal > (
"Range map: Strided block id", -1,
"Strided block id for range map" );
81 validParamList->set< std::string > (
"Domain map: Striding info",
"{}",
"Striding information for domain map");
82 validParamList->set< LocalOrdinal > (
"Domain map: Strided block id", -1,
"Strided block id for domain map" );
84 return validParamList;
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 Input(currentLevel,
"A");
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 const ParameterList& pL = GetParameterList();
95 size_t row = Teuchos::as<size_t>(pL.get<
int>(
"block row"));
96 size_t col = Teuchos::as<size_t>(pL.get<
int>(
"block col"));
98 RCP<Matrix> Ain = Get<RCP<Matrix> >(currentLevel,
"A");
99 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
101 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
102 TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(),
Exceptions::RuntimeError,
"row [" << row <<
"] > A.Rows() [" << A->Rows() <<
"].");
103 TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(),
Exceptions::RuntimeError,
"col [" << col <<
"] > A.Cols() [" << A->Cols() <<
"].");
105 RCP<CrsMatrixWrap> Op = Teuchos::rcp(
new CrsMatrixWrap(A->getMatrix(row, col)));
108 RCP<const StridedMap> srangeMap = Teuchos::null;
109 RCP<const StridedMap> sdomainMap = Teuchos::null;
113 std::vector<size_t> rangeStridingInfo;
114 std::vector<size_t> domainStridingInfo;
115 LocalOrdinal rangeStridedBlockId = 0;
116 LocalOrdinal domainStridedBlockId = 0;
117 bool bRangeUserSpecified = CheckForUserSpecifiedBlockInfo(
true, rangeStridingInfo, rangeStridedBlockId);
118 bool bDomainUserSpecified = CheckForUserSpecifiedBlockInfo(
false, domainStridingInfo, domainStridedBlockId);
119 TEUCHOS_TEST_FOR_EXCEPTION(bRangeUserSpecified != bDomainUserSpecified,
Exceptions::RuntimeError,
"MueLu::SubBlockAFactory[" << row <<
"," << col <<
"]: the user has to specify either both domain and range map or none.");
122 RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
123 RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
125 RCP<const Map> rangeMap = rangeMapExtractor ->getMap(row);
126 RCP<const Map> domainMap = domainMapExtractor->getMap(col);
129 if(bRangeUserSpecified) srangeMap = Teuchos::rcp(
new StridedMap(rangeMap,rangeStridingInfo,rangeMap->getIndexBase(),rangeStridedBlockId,0));
130 else srangeMap = rcp_dynamic_cast<
const StridedMap>(rangeMap);
132 if(bDomainUserSpecified) sdomainMap = Teuchos::rcp(
new StridedMap(domainMap,domainStridingInfo,domainMap->getIndexBase(),domainStridedBlockId,0));
133 else sdomainMap = rcp_dynamic_cast<
const StridedMap>(domainMap);
138 if (srangeMap.is_null()) {
139 RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
140 RCP<const StridedMap> sFullRangeMap = rcp_dynamic_cast<
const StridedMap>(fullRangeMap);
141 TEUCHOS_TEST_FOR_EXCEPTION(sFullRangeMap.is_null(),
Exceptions::BadCast,
"Full rangeMap is not a strided map.");
143 std::vector<size_t> stridedData = sFullRangeMap->getStridingData();
144 if (stridedData.size() == 1 && row > 0) {
146 srangeMap = StridedMapFactory::Build(rangeMap, stridedData, 0, sFullRangeMap->getOffset());
150 srangeMap = StridedMapFactory::Build(rangeMap, stridedData, row, sFullRangeMap->getOffset());
154 if (sdomainMap.is_null()) {
155 RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
156 RCP<const StridedMap> sFullDomainMap = rcp_dynamic_cast<
const StridedMap>(fullDomainMap);
157 TEUCHOS_TEST_FOR_EXCEPTION(sFullDomainMap.is_null(),
Exceptions::BadCast,
"Full domainMap is not a strided map");
159 std::vector<size_t> stridedData = sFullDomainMap->getStridingData();
160 if (stridedData.size() == 1 && col > 0) {
162 sdomainMap = StridedMapFactory::Build(domainMap, stridedData, 0, sFullDomainMap->getOffset());
166 sdomainMap = StridedMapFactory::Build(domainMap, stridedData, col, sFullDomainMap->getOffset());
170 TEUCHOS_TEST_FOR_EXCEPTION(srangeMap.is_null(),
Exceptions::BadCast,
"rangeMap " << row <<
" is not a strided map.");
171 TEUCHOS_TEST_FOR_EXCEPTION(sdomainMap.is_null(),
Exceptions::BadCast,
"domainMap " << col <<
" is not a strided map.");
173 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") has strided maps:" 174 <<
"\n range map fixed block size = " << srangeMap ->getFixedBlockSize() <<
", strided block id = " << srangeMap ->getStridedBlockId()
175 <<
"\n domain map fixed block size = " << sdomainMap->getFixedBlockSize() <<
", strided block id = " << sdomainMap->getStridedBlockId() << std::endl;
177 if (Op->IsView(
"stridedMaps") ==
true)
178 Op->RemoveView(
"stridedMaps");
179 Op->CreateView(
"stridedMaps", srangeMap, sdomainMap);
181 TEUCHOS_TEST_FOR_EXCEPTION(Op->IsView(
"stridedMaps") ==
false,
Exceptions::RuntimeError,
"Failed to set \"stridedMaps\" view.");
183 currentLevel.Set(
"A", rcp_dynamic_cast<Matrix>(Op),
this);
186 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
188 const ParameterList & pL = GetParameterList();
191 stridedBlockId = pL.get<LocalOrdinal>(
"Range map: Strided block id");
193 stridedBlockId = pL.get<LocalOrdinal>(
"Domain map: Strided block id");
198 if(bRange ==
true) str = std::string(
"Range map: Striding info");
199 else str = std::string(
"Domain map: Striding info");
200 if(pL.isParameter(str)) {
201 std::string strStridingInfo = pL.get<std::string>(str);
202 if(strStridingInfo.empty() ==
false) {
203 Teuchos::Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
204 stridingInfo = Teuchos::createVector(arrayVal);
208 if(stridingInfo.size() > 0)
return true;
Exception indicating invalid cast attempted.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
Namespace for MueLu classes and methods.
void DeclareInput(Level ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level ¤tLevel) const
Build an object with this factory.
Class that holds all level-specific information.
RCP< const ParameterList > GetValidParameterList() const
Input.
Exception throws to report errors in the internal logical of the program.
static const RCP< const NoFactory > getRCP()
Static Get() functions.