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" 63 #include "MueLu_HierarchyHelpers.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!");
108 A00->CreateView(
"stridedMaps", bA->getRangeMap(0), bA->getDomainMap(0));
109 A01->CreateView(
"stridedMaps", bA->getRangeMap(0), bA->getDomainMap(1));
110 A10->CreateView(
"stridedMaps", bA->getRangeMap(1), bA->getDomainMap(0));
112 A11->CreateView(
"stridedMaps", bA->getRangeMap(1), bA->getDomainMap(1));
114 const ParameterList& pL = GetParameterList();
115 SC omega = pL.get<Scalar>(
"omega");
118 RCP<Matrix> T = MatrixFactory2::BuildCopy(A01);
120 bool lumping = pL.get<
bool>(
"lumping");
121 bool fixing = pL.get<
bool>(
"fixing");
127 for (
size_t k = 0; k < as<size_t>(D.size()); k++)
128 if (STS::magnitude(D[k]) < 1e-4)
137 for (
size_t k = 0; k < as<size_t>(D.size()); k++)
141 RCP<Matrix> S = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A10,
false, *T,
false, GetOStream(
Statistics2));
143 if (!A11.is_null()) {
145 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*A11,
false, one, *S,
false, one, T, GetOStream(
Statistics2));
152 Set(currentLevel,
"A", S);
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Exception indicating invalid cast attempted.
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
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)
Namespace for MueLu classes and methods.
void DeclareInput(Level ¤tLevel) const
Input.
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.
void Build(Level ¤tLevel) const
Build an object with this factory.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Crs2Op(RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static const RCP< const NoFactory > getRCP()
Static Get() functions.