53 #ifndef MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 54 #define MUELU_SUBBLOCKAFACTORY_DEF_HPP_ 59 #include <Xpetra_BlockedCrsMatrix.hpp> 60 #include <Xpetra_MapExtractor.hpp> 61 #include <Xpetra_Matrix.hpp> 62 #include <Xpetra_StridedMapFactory.hpp> 69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 RCP<ParameterList> validParamList = rcp(
new ParameterList());
74 validParamList->set<
int > (
"block row", 0,
"Block row of subblock matrix A");
75 validParamList->set<
int > (
"block col", 0,
"Block column of subblock matrix A");
77 validParamList->set< std::string > (
"Range map: Striding info",
"{}",
"Striding information for range map");
78 validParamList->set<
LocalOrdinal > (
"Range map: Strided block id", -1,
"Strided block id for range map" );
79 validParamList->set< std::string > (
"Domain map: Striding info",
"{}",
"Striding information for domain map");
80 validParamList->set<
LocalOrdinal > (
"Domain map: Strided block id", -1,
"Strided block id for domain map" );
82 return validParamList;
85 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 Input(currentLevel,
"A");
90 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 = currentLevel.
Get< RCP<Matrix> >(
"A",this->GetFactory(
"A").get());
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() <<
"].");
106 RCP<Matrix> Op = A->getMatrix(row, col);
113 RCP<BlockedCrsMatrix> bOp = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Op);
114 if (bOp != Teuchos::null) {
116 if (bOp->Rows() == 1 && bOp->Cols() == 1) {
118 Op = bOp->getCrsMatrix();
119 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null,
Exceptions::BadCast,
"SubBlockAFactory::Build: sub block A[" << row <<
"," << col <<
"] must be a single block CrsMatrixWrap object!");
124 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") is a " << bOp->Rows() <<
"x" << bOp->Cols() <<
" block matrix" << std::endl;
125 GetOStream(
Statistics2) <<
"with altogether " << bOp->getGlobalNumRows() <<
"x" << bOp->getGlobalNumCols() <<
" rows and columns." << std::endl;
126 currentLevel.
Set(
"A", Op,
this);
142 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null,
Exceptions::BadCast,
"SubBlockAFactory::Build: sub block A[" << row <<
"," << col <<
"] is NOT a BlockedCrsMatrix but also NOT a CrsMatrixWrap object? This cannot be.");
145 RCP<const StridedMap> srangeMap = Teuchos::null;
146 RCP<const StridedMap> sdomainMap = Teuchos::null;
150 std::vector<size_t> rangeStridingInfo;
151 std::vector<size_t> domainStridingInfo;
154 bool bRangeUserSpecified = CheckForUserSpecifiedBlockInfo(
true, rangeStridingInfo, rangeStridedBlockId);
155 bool bDomainUserSpecified = CheckForUserSpecifiedBlockInfo(
false, domainStridingInfo, domainStridedBlockId);
156 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.");
159 RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
160 RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
162 RCP<const Map> rangeMap = rangeMapExtractor ->getMap(row);
163 RCP<const Map> domainMap = domainMapExtractor->getMap(col);
166 if(bRangeUserSpecified) srangeMap = Teuchos::rcp(
new StridedMap(rangeMap,rangeStridingInfo,rangeMap->getIndexBase(),rangeStridedBlockId,0));
167 else srangeMap = rcp_dynamic_cast<
const StridedMap>(rangeMap);
169 if(bDomainUserSpecified) sdomainMap = Teuchos::rcp(
new StridedMap(domainMap,domainStridingInfo,domainMap->getIndexBase(),domainStridedBlockId,0));
170 else sdomainMap = rcp_dynamic_cast<
const StridedMap>(domainMap);
175 if (srangeMap.is_null()) {
176 RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
177 RCP<const StridedMap> sFullRangeMap = rcp_dynamic_cast<
const StridedMap>(fullRangeMap);
178 TEUCHOS_TEST_FOR_EXCEPTION(sFullRangeMap.is_null(),
Exceptions::BadCast,
"Full rangeMap is not a strided map.");
180 std::vector<size_t> stridedData = sFullRangeMap->getStridingData();
181 if (stridedData.size() == 1 && row > 0) {
183 srangeMap = StridedMapFactory::Build(rangeMap, stridedData, 0, sFullRangeMap->getOffset());
187 srangeMap = StridedMapFactory::Build(rangeMap, stridedData, row, sFullRangeMap->getOffset());
191 if (sdomainMap.is_null()) {
192 RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
193 RCP<const StridedMap> sFullDomainMap = rcp_dynamic_cast<
const StridedMap>(fullDomainMap);
194 TEUCHOS_TEST_FOR_EXCEPTION(sFullDomainMap.is_null(),
Exceptions::BadCast,
"Full domainMap is not a strided map");
196 std::vector<size_t> stridedData = sFullDomainMap->getStridingData();
197 if (stridedData.size() == 1 && col > 0) {
199 sdomainMap = StridedMapFactory::Build(domainMap, stridedData, 0, sFullDomainMap->getOffset());
203 sdomainMap = StridedMapFactory::Build(domainMap, stridedData, col, sFullDomainMap->getOffset());
207 TEUCHOS_TEST_FOR_EXCEPTION(srangeMap.is_null(),
Exceptions::BadCast,
"rangeMap " << row <<
" is not a strided map.");
208 TEUCHOS_TEST_FOR_EXCEPTION(sdomainMap.is_null(),
Exceptions::BadCast,
"domainMap " << col <<
" is not a strided map.");
210 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") is a single block and has strided maps:" 211 <<
"\n range map fixed block size = " << srangeMap ->getFixedBlockSize() <<
", strided block id = " << srangeMap ->getStridedBlockId()
212 <<
"\n domain map fixed block size = " << sdomainMap->getFixedBlockSize() <<
", strided block id = " << sdomainMap->getStridedBlockId() << std::endl;
213 GetOStream(
Statistics2) <<
"A(" << row <<
"," << col <<
") has " << Op->getGlobalNumRows() <<
"x" << Op->getGlobalNumCols() <<
" rows and columns." << std::endl;
216 if (Op->IsView(
"stridedMaps") ==
true)
217 Op->RemoveView(
"stridedMaps");
218 Op->CreateView(
"stridedMaps", srangeMap, sdomainMap);
220 TEUCHOS_TEST_FOR_EXCEPTION(Op->IsView(
"stridedMaps") ==
false,
Exceptions::RuntimeError,
"Failed to set \"stridedMaps\" view.");
222 currentLevel.
Set(
"A", Op,
this);
225 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
227 const ParameterList & pL = GetParameterList();
230 stridedBlockId = pL.get<
LocalOrdinal>(
"Range map: Strided block id");
232 stridedBlockId = pL.get<
LocalOrdinal>(
"Domain map: Strided block id");
237 if(bRange ==
true) str = std::string(
"Range map: Striding info");
238 else str = std::string(
"Domain map: Striding info");
239 if(pL.isParameter(str)) {
240 std::string strStridingInfo = pL.get<std::string>(str);
241 if(strStridingInfo.empty() ==
false) {
242 Teuchos::Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
243 stridingInfo = Teuchos::createVector(arrayVal);
247 if(stridingInfo.size() > 0)
return true;
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
RCP< const ParameterList > GetValidParameterList() const
Input.
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
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print even more statistics.
void Build(Level ¤tLevel) const
Build an object with this factory.
Class that holds all level-specific information.
void DeclareInput(Level ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
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.
static const RCP< const NoFactory > getRCP()
Static Get() functions.