47 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
48 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
50 #include <Teuchos_ArrayViewDecl.hpp>
57 #include <Xpetra_MultiVectorFactory.hpp>
63 #include "MueLu_Utilities.hpp"
65 #include "MueLu_HierarchyUtils.hpp"
68 #include "MueLu_FactoryManager.hpp"
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 FactManager_ = FactManager;
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 validParamList->
set<
bool>(
"lumping",
false,
87 "Use lumping to construct diag(A(0,0)), "
88 "i.e. use row sum of the abs values on the diagonal as approximation of A00 (and A00^{-1})");
89 validParamList->
set<
SC>(
"Damping factor", one,
"Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
90 validParamList->
set<
LO>(
"Sweeps", 1,
"Number of BraessSarazin sweeps (default = 1)");
91 validParamList->
set<
bool>(
"q2q1 mode",
false,
"Run in the mode matching Q2Q1 matlab code");
93 return validParamList;
96 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 this->Input(currentLevel,
"A");
101 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! "
102 "Introduce a FactoryManager for the SchurComplement equation.");
109 currentLevel.
DeclareInput(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
112 currentLevel.
DeclareInput(
"A", FactManager_->GetFactory(
"A").get());
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 this->GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
124 A_ = Factory::Get<RCP<Matrix> >(currentLevel,
"A");
127 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
130 rangeMapExtractor_ = bA->getRangeMapExtractor();
131 domainMapExtractor_ = bA->getDomainMapExtractor();
134 A00_ = bA->getMatrix(0, 0);
135 A01_ = bA->getMatrix(0, 1);
136 A10_ = bA->getMatrix(1, 0);
137 A11_ = bA->getMatrix(1, 1);
140 SC omega = pL.
get<
SC>(
"Damping factor");
144 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
145 if (pL.
get<
bool>(
"lumping") ==
false) {
146 A00_->getLocalDiagCopy(*diagFVector);
150 diagFVector->scale(omega);
155 if (bD.
is_null() ==
false && bD->getBlockedMap()->getNumMaps() == 1) {
156 RCP<Vector> nestedVec = bD->getMultiVector(0, bD->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
164 smoo_ = currentLevel.Get<
RCP<SmootherBase> >(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
165 S_ = currentLevel.Get<
RCP<Matrix> >(
"A", FactManager_->GetFactory(
"A").get());
171 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
174 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
180 bool bCopyResultX =
false;
181 bool bReorderX =
false;
182 bool bReorderB =
false;
192 if (rbA.is_null() ==
false) {
197 if (bX.is_null() ==
true) {
199 rcpX.
swap(vectorWithBlockedMap);
205 if (bB.is_null() ==
true) {
207 rcpB.
swap(vectorWithBlockedMap);
212 if (bX.is_null() ==
true) {
214 rcpX.
swap(vectorWithBlockedMap);
218 if (bB.is_null() ==
true) {
220 rcpB.
swap(vectorWithBlockedMap);
225 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
226 bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
229 if (rbA.is_null() ==
false) {
233 if (bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) {
236 rcpX.
swap(reorderedBlockedVector);
239 if (bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) {
242 rcpB.
swap(reorderedBlockedVector);
249 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
250 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
252 RCP<MultiVector> deltaX = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
257 RCP<MultiVector> Rtmp = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
260 SC one = STS::one(), zero = STS::zero();
264 LO nSweeps = pL.
get<
LO>(
"Sweeps");
267 if (InitialGuessIsZero) {
268 R = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
269 R->update(one, *rcpB, zero);
275 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
276 S_->getLocalDiagCopy(*diagSVector);
278 for (
LO run = 0; run < nSweeps; ++run) {
282 RCP<MultiVector> R0 = rangeMapExtractor_->ExtractVector(R, 0, bRangeThyraMode);
283 RCP<MultiVector> R1 = rangeMapExtractor_->ExtractVector(R, 1, bRangeThyraMode);
286 deltaX0->putScalar(zero);
287 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
288 A10_->apply(*deltaX0, *Rtmp);
289 Rtmp->update(one, *R1, -one);
291 if (!pL.
get<
bool>(
"q2q1 mode")) {
292 deltaX1->putScalar(zero);
295 if (Teuchos::rcp_dynamic_cast<BlockedVector>(diagSVector) == Teuchos::null) {
299 for (
GO row = 0; row < deltaX1data.
size(); row++)
300 deltaX1data[row] = Teuchos::as<SC>(1.1) * Rtmpdata[row] / Sdiag[row];
306 smoo_->Apply(*deltaX1, *Rtmp);
309 deltaX0->putScalar(zero);
310 A01_->apply(*deltaX1, *deltaX0);
311 deltaX0->update(one, *R0, -one);
313 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
315 RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
316 RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
319 X0->update(one, *deltaX0, one);
320 X1->update(one, *deltaX1, one);
322 domainMapExtractor_->InsertVector(X0, 0, rcpX, bDomainThyraMode);
323 domainMapExtractor_->InsertVector(X1, 1, rcpX, bDomainThyraMode);
325 if (run < nSweeps - 1) {
330 if (bCopyResultX ==
true) {
332 X.update(one, *Xmerged, zero);
336 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
345 std::ostringstream out;
347 out <<
"{type = " << type_ <<
"}";
351 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
356 out0 <<
"Prec. type: " << type_ << std::endl;
359 if (verbLevel &
Debug) {
364 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
Important warning messages (one line)
virtual const Teuchos::ParameterList & GetParameterList() const
Exception indicating invalid cast attempted.
BraessSarazin smoother for 2x2 block matrices.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
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)
Print additional debugging information.
RCP< SmootherPrototype > Copy() const
void DeclareInput(Level ¤tLevel) const
Input.
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Class that holds all level-specific information.
bool IsSetup() const
Get the state of a smoother prototype.
RCP< const ParameterList > GetValidParameterList() const
Input.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
void Setup(Level ¤tLevel)
Setup routine.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
Exception throws to report errors in the internal logical of the program.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
An exception safe way to call the method 'Level::SetFactoryManager()'.
Teuchos::RCP< const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > buildReorderedBlockedMultiVector(Teuchos::RCP< const Xpetra::BlockReorderManager > brm, Teuchos::RCP< const Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > bvec)
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
std::string description() const
Return a simple one-line description of this object.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)