46 #ifndef MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 47 #define MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 49 #include <Xpetra_BlockedCrsMatrix.hpp> 50 #include <Xpetra_MultiVectorFactory.hpp> 51 #include <Xpetra_VectorFactory.hpp> 52 #include <Xpetra_MatrixFactory.hpp> 53 #include <Xpetra_Matrix.hpp> 54 #include <Xpetra_MatrixMatrix.hpp> 55 #include <Xpetra_CrsMatrixWrap.hpp> 56 #include <Xpetra_BlockedCrsMatrix.hpp> 57 #include <Xpetra_CrsMatrix.hpp> 60 #include "MueLu_Utilities.hpp" 63 #include "MueLu_SchurComplementFactory.hpp" 67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 RCP<ParameterList> validParamList = rcp(
new ParameterList());
71 SC one = Teuchos::ScalarTraits<SC>::one();
73 validParamList->set<RCP<const FactoryBase> >(
"A",
NoFactory::getRCP(),
"Generating factory of the matrix A used for building Schur complement\n" 74 "(must be a 2x2 blocked operator)");
75 validParamList->set<SC> (
"omega", one,
"Scaling parameter in S = A(1,1) - 1/omega A(1,0) diag{A(0,0)}^{-1} A(0,1)");
76 validParamList->set<
bool> (
"lumping",
false,
"Use lumping to construct diag(A(0,0)), i.e. use row sum of the abs values on the diagonal " 77 "as approximation of A00 (and A00^{-1})");
78 validParamList->set<
bool> (
"fixing",
false,
"Fix diagonal by replacing small entries with 1.0");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 Input(currentLevel,
"A");
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 typedef Teuchos::ScalarTraits<SC> STS;
93 SC zero = STS::zero(), one = STS::one();
95 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
96 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
98 "MueLu::SchurComplementFactory::Build: input matrix A is not of type BlockedCrsMatrix!");
101 "MueLu::SchurComplementFactory::Build: input matrix A is a " << bA->Rows() <<
"x" << bA->Cols() <<
" block matrix. We expect a 2x2 blocked operator.");
103 RCP<Matrix> A00 = bA->getMatrix(0,0);
104 RCP<Matrix> A01 = bA->getMatrix(0,1);
105 RCP<Matrix> A10 = bA->getMatrix(1,0);
106 RCP<Matrix> A11 = bA->getMatrix(1,1);
108 RCP<BlockedCrsMatrix> bA01 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A01);
109 bool bIsBlocked = (bA01 == Teuchos::null ? false :
true);
111 const ParameterList& pL = GetParameterList();
112 SC omega = pL.get<Scalar>(
"omega");
115 "MueLu::SchurComplementFactory::Build: Scaling parameter omega must not be zero to avoid division by zero.");
117 RCP<Matrix> S = Teuchos::null;
119 if(A01.is_null() ==
false && A10.is_null() ==
false) {
120 bool lumping = pL.get<
bool>(
"lumping");
121 bool fixing = pL.get<
bool>(
"fixing");
122 RCP<Vector> diag = Teuchos::null;
124 diag = VectorFactory::Build(A00->getRangeMap(),
true);
125 A00->getLocalDiagCopy(*diag);
127 RCP<BlockedCrsMatrix> bA00 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A00);
128 TEUCHOS_TEST_FOR_EXCEPTION(bA00.is_null()==
false,
MueLu::Exceptions::RuntimeError,
"MueLu::SchurComplementFactory::Build: Mass lumping not implemented. Implement a mass lumping kernel!");
134 D->scale(Teuchos::as<Scalar>(-one/omega));
137 RCP<Matrix> T = MatrixFactory::BuildCopy(A01);
143 "MueLu::SchurComplementFactory::Build: RangeMap of A01 and domain map of A10 are not the same.");
144 RCP<ParameterList> myparams = rcp(
new ParameterList);
145 S = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A10,
false, *T,
false, GetOStream(
Statistics2),
true,
true,std::string(
"SchurComplementFactory"),myparams);
149 RCP<BlockedCrsMatrix> bA10 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A10);
150 RCP<BlockedCrsMatrix> bT = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(T);
153 "MueLu::SchurComplementFactory::Build: Block rows and cols of A01 and A10 are not compatible.");
155 "MueLu::SchurComplementFactory::Build: The scaled A01 operator has " << bT->Rows() <<
"x" << bT->Cols() <<
" blocks, " 156 "but should have " << bA01->Rows() <<
"x" << bA01->Cols() <<
" blocks.");
158 "MueLu::SchurComplementFactory::Build: Block rows and cols of A01 and A10 are not compatible.");
160 S = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(*bA10,
false, *bT,
false, GetOStream(
Statistics2));
163 if (!A11.is_null()) {
165 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*A11,
false, one, *S,
false, one, T, GetOStream(
Statistics2));
170 "MueLu::SchurComplementFactory::Build: RangeMap of A11 and S are not the same.");
172 "MueLu::SchurComplementFactory::Build: DomainMap of A11 and S are not the same.");
177 if (!A11.is_null()) {
178 S = MatrixFactory::BuildCopy(A11);
180 S = MatrixFactory::Build(A11->getRowMap(), 10 );
181 S->fillComplete(A11->getDomainMap(),A11->getRangeMap());
192 RCP<BlockedCrsMatrix> bS = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(S);
194 if (bS != Teuchos::null && bS->Rows() == 1 && bS->Cols() == 1) {
195 RCP<Matrix> temp = bS->getCrsMatrix();
201 Set(currentLevel,
"A", S);
Exception indicating invalid cast attempted.
Timer to be used in factories. Similar to Monitor but with additional timers.
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
void DeclareInput(Level ¤tLevel) const
Input.
Namespace for MueLu classes and methods.
Print even more statistics.
Class that holds all level-specific information.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetInverse(Teuchos::RCP< const Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())
void Build(Level ¤tLevel) const
Build an object with this factory.
Exception throws to report errors in the internal logical of the program.
static const RCP< const NoFactory > getRCP()
Static Get() functions.