47#ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
48#define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
50#include <Teuchos_ArrayViewDecl.hpp>
51#include <Teuchos_ScalarTraits.hpp>
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MultiVectorFactory.hpp>
58#include <Xpetra_VectorFactory.hpp>
59#include <Xpetra_ReorderedBlockedCrsMatrix.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>
74 TEUCHOS_TEST_FOR_EXCEPTION(pos != 0,
Exceptions::RuntimeError,
"MueLu::BraessSarazinSmoother::AddFactoryManager: parameter \'pos\' must be zero! error.");
75 FactManager_ = FactManager;
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 RCP<ParameterList> validParamList = rcp(
new ParameterList());
82 SC one = Teuchos::ScalarTraits<SC>::one();
84 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
86 validParamList->set<
bool>(
"lumping",
false,
"Use lumping to construct diag(A(0,0)), "
87 "i.e. use row sum of the abs values on the diagonal as approximation of A00 (and A00^{-1})");
88 validParamList->set<SC>(
"Damping factor", one,
"Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
89 validParamList->set<LO>(
"Sweeps", 1,
"Number of BraessSarazin sweeps (default = 1)");
90 validParamList->set<
bool>(
"q2q1 mode",
false,
"Run in the mode matching Q2Q1 matlab code");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
97 this->Input(currentLevel,
"A");
100 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! "
101 "Introduce a FactoryManager for the SchurComplement equation.");
108 currentLevel.
DeclareInput(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
111 currentLevel.
DeclareInput(
"A", FactManager_->GetFactory(
"A").get());
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 this->GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
123 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
124 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
126 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
129 rangeMapExtractor_ = bA->getRangeMapExtractor();
130 domainMapExtractor_ = bA->getDomainMapExtractor();
133 A00_ = bA->getMatrix(0,0);
134 A01_ = bA->getMatrix(0,1);
135 A10_ = bA->getMatrix(1,0);
136 A11_ = bA->getMatrix(1,1);
139 SC omega = pL.get<SC>(
"Damping factor");
143 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
144 if (pL.get<
bool>(
"lumping") ==
false) {
145 A00_->getLocalDiagCopy(*diagFVector);
149 diagFVector->scale(omega);
153 RCP<BlockedVector> bD = Teuchos::rcp_dynamic_cast<BlockedVector>(D_);
154 if(bD.is_null() ==
false && bD->getBlockedMap()->getNumMaps() == 1) {
155 RCP<Vector> nestedVec = bD->getMultiVector(0,bD->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
163 smoo_ = currentLevel.
Get<RCP<SmootherBase> >(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
164 S_ = currentLevel.
Get<RCP<Matrix> > (
"A", FactManager_->GetFactory(
"A").get());
170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
173 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
175 RCP<MultiVector> rcpX = rcpFromRef(X);
176 RCP<const MultiVector> rcpB = rcpFromRef(B);
179 bool bCopyResultX =
false;
180 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
182 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
183 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
185 if(bX.is_null() ==
true) {
186 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
191 if(bB.is_null() ==
true) {
192 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
197 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
198 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
201 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
202 if(rbA.is_null() ==
false) {
204 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
207 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
209 Teuchos::RCP<MultiVector> test =
210 buildReorderedBlockedMultiVector(brm, bX);
213 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
215 Teuchos::RCP<const MultiVector> test =
216 buildReorderedBlockedMultiVector(brm, bB);
224 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
225 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
227 RCP<MultiVector> deltaX = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
228 RCP<BlockedMultiVector> bdeltaX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(deltaX);
229 RCP<MultiVector> deltaX0 = bdeltaX->getMultiVector(0,bDomainThyraMode);
230 RCP<MultiVector> deltaX1 = bdeltaX->getMultiVector(1,bDomainThyraMode);
232 RCP<MultiVector> Rtmp = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
235 typedef Teuchos::ScalarTraits<SC> STS;
236 SC one = STS::one(), zero = STS::zero();
240 LO nSweeps = pL.get<LO>(
"Sweeps");
243 if (InitialGuessIsZero) {
244 R = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
245 R->update(one, *rcpB, zero);
251 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
252 S_->getLocalDiagCopy(*diagSVector);
254 for (LO run = 0; run < nSweeps; ++run) {
258 RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0, bRangeThyraMode);
259 RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1, bRangeThyraMode);
262 deltaX0->putScalar(zero);
263 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
264 A10_->apply(*deltaX0, *Rtmp);
265 Rtmp->update(one, *R1, -one);
267 if (!pL.get<
bool>(
"q2q1 mode")) {
268 deltaX1->putScalar(zero);
271 if(Teuchos::rcp_dynamic_cast<BlockedVector>(diagSVector) == Teuchos::null) {
272 ArrayRCP<SC> Sdiag = diagSVector->getDataNonConst(0);
273 ArrayRCP<SC> deltaX1data = deltaX1->getDataNonConst(0);
274 ArrayRCP<SC> Rtmpdata = Rtmp->getDataNonConst(0);
275 for (GO row = 0; row < deltaX1data.size(); row++)
276 deltaX1data[row] = Teuchos::as<SC>(1.1)*Rtmpdata[row] / Sdiag[row];
282 smoo_->Apply(*deltaX1,*Rtmp);
285 deltaX0->putScalar(zero);
286 A01_->apply(*deltaX1, *deltaX0);
287 deltaX0->update(one, *R0, -one);
289 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
291 RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
292 RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
295 X0->update(one, *deltaX0, one);
296 X1->update(one, *deltaX1, one);
298 domainMapExtractor_->InsertVector(X0, 0, rcpX, bDomainThyraMode);
299 domainMapExtractor_->InsertVector(X1, 1, rcpX, bDomainThyraMode);
301 if (run < nSweeps-1) {
307 if (bCopyResultX ==
true) {
308 RCP<MultiVector> Xmerged = bX->Merge();
309 X.update(one, *Xmerged, zero);
314 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
315 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
320 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
323 std::ostringstream out;
325 out <<
"{type = " << type_ <<
"}";
329 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
334 out0 <<
"Prec. type: " << type_ << std::endl;
337 if (verbLevel &
Debug) {
342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
345 return Teuchos::OrdinalTraits<size_t>::invalid();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
BraessSarazin smoother for 2x2 block matrices.
RCP< const ParameterList > GetValidParameterList() const
Input.
std::string description() const
Return a simple one-line description of this object.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
void DeclareInput(Level ¤tLevel) const
Input.
RCP< SmootherPrototype > Copy() const
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void Setup(Level ¤tLevel)
Setup routine.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual const Teuchos::ParameterList & GetParameterList() const =0
An exception safe way to call the method 'Level::SetFactoryManager()'.
bool IsSetup() const
Get the state of a smoother prototype.
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
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())
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.