46 #ifndef MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 47 #define MUELU_SCHURCOMPLEMENTFACTORY_DEF_HPP_ 49 #ifdef HAVE_MUELU_EXPERIMENTAL 51 #include <Xpetra_BlockedCrsMatrix.hpp> 52 #include <Xpetra_MultiVectorFactory.hpp> 53 #include <Xpetra_VectorFactory.hpp> 54 #include <Xpetra_MatrixFactory.hpp> 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_MatrixMatrix.hpp> 57 #include <Xpetra_CrsMatrixWrap.hpp> 58 #include <Xpetra_BlockedCrsMatrix.hpp> 59 #include <Xpetra_CrsMatrix.hpp> 62 #include "MueLu_Utilities.hpp" 65 #include "MueLu_SchurComplementFactory.hpp" 69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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>
95 SC zero = STS::zero(), one = STS::one();
97 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"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.");
111 bool bIsBlocked = (bA01 == Teuchos::null ? false :
true);
114 SC omega = pL.
get<Scalar>(
"omega");
117 "MueLu::SchurComplementFactory::Build: Scaling parameter omega must not be zero to avoid division by zero.");
122 bool lumping = pL.
get<
bool>(
"lumping");
123 bool fixing = pL.
get<
bool>(
"fixing");
127 diag = VectorFactory::Build(A00->getRangeMap(),
true);
128 A00->getLocalDiagCopy(*diag);
138 D->scale(Teuchos::as<Scalar>(-STS::one()/omega));
148 "MueLu::SchurComplementFactory::Build: RangeMap of A01 and domain map of A10 are not the same.");
150 S = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A10,
false, *T,
false, GetOStream(
Statistics2));
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.");
164 S = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixMultiplyBlock(*bA10,
false, *bT,
false, GetOStream(
Statistics2));
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.");
182 S = MatrixFactory::BuildCopy(A11);
184 S = MatrixFactory::Build(A11->getRowMap(), 10 );
185 S->fillComplete(A11->getDomainMap(),A11->getRangeMap());
198 if (bS != Teuchos::null && bS->Rows() == 1 && bS->Cols() == 1) {
205 Set(currentLevel,
"A", S);
Exception indicating invalid cast attempted.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void swap(RCP< T > &r_ptr)
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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.