47 #ifndef MUELU_ZEROSUBBLOCKAFACTORY_DEF_HPP_ 48 #define MUELU_ZEROSUBBLOCKAFACTORY_DEF_HPP_ 53 #include <Xpetra_BlockedCrsMatrix.hpp> 54 #include <Xpetra_Matrix.hpp> 55 #include <Xpetra_MatrixFactory.hpp> 56 #include <Xpetra_Vector.hpp> 57 #include <Xpetra_VectorFactory.hpp> 64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 RCP<ParameterList> validParamList = rcp(
new ParameterList());
70 validParamList->set<
int>(
"block row", 0,
"Block row of subblock in block matrix A");
71 validParamList->set<
int>(
"block col", 0,
"Block column of subblock in block matrix A");
73 return validParamList;
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 Input(currentLevel,
"A");
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 const ParameterList& pL = GetParameterList();
86 const size_t row = Teuchos::as<size_t>(pL.get<
int>(
"block row"));
87 const size_t col = Teuchos::as<size_t>(pL.get<
int>(
"block col"));
90 RCP<Matrix> Ain = currentLevel.
Get<RCP<Matrix>>(
"A", this->GetFactory(
"A").get());
91 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
94 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
95 TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(),
Exceptions::RuntimeError,
"row [" << row <<
"] > A.Rows() [" << A->Rows() <<
"].");
96 TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(),
Exceptions::RuntimeError,
"col [" << col <<
"] > A.Cols() [" << A->Cols() <<
"].");
99 RCP<const StridedMap> stridedRangeMap = Teuchos::null;
100 RCP<const StridedMap> stridedDomainMap = Teuchos::null;
103 RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
104 RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
109 if (stridedRangeMap.is_null()) {
111 "Range map extractor contains non-strided maps in block row " << row <<
". This should not be.");
112 stridedRangeMap = rcp_dynamic_cast<
const StridedMap>(rangeMapExtractor->getMap(row));
115 if (stridedDomainMap.is_null()) {
117 "Domain map extractor contains non-strided maps in block row " << row <<
". This should not be.");
118 stridedDomainMap = rcp_dynamic_cast<
const StridedMap>(domainMapExtractor->getMap(col));
121 TEUCHOS_TEST_FOR_EXCEPTION(stridedRangeMap.is_null(),
Exceptions::BadCast,
"rangeMap " << row <<
" is not a strided map.");
122 TEUCHOS_TEST_FOR_EXCEPTION(stridedDomainMap.is_null(),
Exceptions::BadCast,
"domainMap " << col <<
" is not a strided map.");
124 RCP<Matrix> Op = MatrixFactory::Build(stridedRangeMap, stridedDomainMap, static_cast<size_t>(0));
125 TEUCHOS_ASSERT(!Op.is_null());
127 Op->fillComplete(stridedDomainMap, stridedRangeMap);
128 TEUCHOS_ASSERT(Op->isFillComplete());
130 GetOStream(
Statistics1) <<
"A(" << row <<
"," << col <<
") is a single block and has strided maps:" 131 <<
"\n range map fixed block size = " << stridedRangeMap ->getFixedBlockSize() <<
", strided block id = " << stridedRangeMap ->getStridedBlockId()
132 <<
"\n domain map fixed block size = " << stridedDomainMap->getFixedBlockSize() <<
", strided block id = " << stridedDomainMap->getStridedBlockId() << std::endl;
133 GetOStream(
Statistics2) <<
"A(" << row <<
"," << col <<
") has " << Op->getGlobalNumRows() <<
"x" << Op->getGlobalNumCols() <<
" rows and columns." << std::endl;
135 if (Op->IsView(
"stridedMaps") ==
true)
136 Op->RemoveView(
"stridedMaps");
137 Op->CreateView(
"stridedMaps", stridedRangeMap, stridedDomainMap);
139 currentLevel.
Set(
"A", Op,
this);
Exception indicating invalid cast attempted.
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.
void DeclareInput(Level ¤tLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level ¤tLevel) const override
Build a zero sub-block object with this factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print even more statistics.
Class that holds all level-specific information.
RCP< const ParameterList > GetValidParameterList() const override
Input.
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.