46 #ifndef MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 47 #define MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 49 #ifdef HAVE_MUELU_EXPERIMENTAL 62 #include "MueLu_Utilities.hpp" 65 #include "MueLu_SchurComplementFactory.hpp" 69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
71 RCP<ParameterList> validParamList = rcp(
new ParameterList());
73 SC one = Teuchos::ScalarTraits<SC>::one();
75 validParamList->set<RCP<const FactoryBase> >(
"A",
NoFactory::getRCP(),
"Generating factory of the matrix A used for building Schur complement\n" 76 "(must be a 2x2 blocked operator)");
77 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)");
78 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 " 79 "as approximation of A00 (and A00^{-1})");
80 validParamList->set<
bool> (
"fixing",
false,
"Fix diagonal by replacing small entries with 1.0");
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 typedef Teuchos::ScalarTraits<SC> STS;
95 SC zero = STS::zero(), one = STS::one();
97 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
98 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A);
100 "MueLu::SchurComplementFactory::Build: input matrix A is not of type BlockedCrsMatrix!");
103 "MueLu::SchurComplementFactory::Build: input matrix A is a " << bA->Rows() <<
"x" << bA->Cols() <<
" block matrix. We expect a 2x2 blocked operator.");
105 RCP<Matrix> A00 = bA->getMatrix(0,0);
106 RCP<Matrix> A01 = bA->getMatrix(0,1);
107 RCP<Matrix> A10 = bA->getMatrix(1,0);
108 RCP<Matrix> A11 = bA->getMatrix(1,1);
110 RCP<BlockedCrsMatrix> bA01 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A01);
111 bool bIsBlocked = (bA01 == Teuchos::null ? false :
true);
113 const ParameterList& pL = GetParameterList();
114 SC omega = pL.get<Scalar>(
"omega");
117 "MueLu::SchurComplementFactory::Build: Scaling parameter omega must not be zero to avoid division by zero.");
119 RCP<Matrix> S = Teuchos::null;
121 if(A01.is_null() ==
false && A10.is_null() ==
false) {
122 bool lumping = pL.get<
bool>(
"lumping");
123 bool fixing = pL.get<
bool>(
"fixing");
125 RCP<Vector> diag = Teuchos::null;
127 diag = VectorFactory::Build(A00->getRangeMap(),
true);
128 A00->getLocalDiagCopy(*diag);
138 D->scale(Teuchos::as<Scalar>(-STS::one()/omega));
142 RCP<Matrix> T = MatrixFactory::BuildCopy(A01);
148 "MueLu::SchurComplementFactory::Build: RangeMap of A01 and domain map of A10 are not the same.");
153 RCP<BlockedCrsMatrix> bA10 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A10);
154 RCP<BlockedCrsMatrix> bT = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(T);
157 "MueLu::SchurComplementFactory::Build: Block rows and cols of A01 and A10 are not compatible.");
159 "MueLu::SchurComplementFactory::Build: The scaled A01 operator has " << bT->Rows() <<
"x" << bT->Cols() <<
" blocks, " 160 "but should have " << bA01->Rows() <<
"x" << bA01->Cols() <<
" blocks.");
162 "MueLu::SchurComplementFactory::Build: Block rows and cols of A01 and A10 are not compatible.");
167 if (!A11.is_null()) {
169 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*A11,
false, one, *S,
false, one, T, GetOStream(
Statistics2));
174 "MueLu::SchurComplementFactory::Build: RangeMap of A11 and S are not the same.");
176 "MueLu::SchurComplementFactory::Build: DomainMap of A11 and S are not the same.");
181 if (!A11.is_null()) {
182 S = MatrixFactory::BuildCopy(A11);
184 S = MatrixFactory::Build(A11->getRowMap(), 10 );
185 S->fillComplete(A11->getDomainMap(),A11->getRangeMap());
196 RCP<BlockedCrsMatrix> bS = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(S);
198 if (bS != Teuchos::null && bS->Rows() == 1 && bS->Cols() == 1) {
199 RCP<Matrix> temp = bS->getCrsMatrix();
205 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 RCP< BlockedCrsMatrix > TwoMatrixMultiplyBlock(const BlockedCrsMatrix &A, bool transposeA, const BlockedCrsMatrix &B, bool transposeB, Teuchos::FancyOStream &fos, bool doFillComplete=true, bool doOptimizeStorage=true)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
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())
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
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.