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;
MueLu::DefaultLocalOrdinal LocalOrdinal
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const RCP< const NoFactory > getRCP()
Static Get() functions.
void DeclareInput(Level ¤tLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
void Build(Level ¤tLevel) const override
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const override
Input.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.