47#ifndef MUELU_UZAWASMOOTHER_DEF_HPP_
48#define MUELU_UZAWASMOOTHER_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 RCP<ParameterList> validParamList = rcp(
new ParameterList());
76 validParamList->set<RCP<const FactoryBase>>(
"A", Teuchos::null,
"Generating factory of the matrix A");
77 validParamList->set<
Scalar> (
"Damping factor", 1.0,
"Damping/Scaling factor in SIMPLE");
78 validParamList->set<
LocalOrdinal>(
"Sweeps", 1,
"Number of SIMPLE sweeps (default = 1)");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::UzawaSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
87 size_t myPos = Teuchos::as<size_t>(pos);
89 if (myPos < FactManager_.size()) {
91 FactManager_.at(myPos) = FactManager;
92 }
else if( myPos == FactManager_.size()) {
94 FactManager_.push_back(FactManager);
96 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
97 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
100 FactManager_.push_back(FactManager);
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
108 AddFactoryManager(FactManager, 0);
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 AddFactoryManager(FactManager, 1);
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").get());
120 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.size() != 2,
Exceptions::RuntimeError,
"MueLu::UzawaSmoother::DeclareInput: You have to declare two FactoryManagers with a \"Smoother\" object: One for predicting the primary variable and one for the SchurComplement system. The smoother for the SchurComplement system needs a SchurComplementFactory as input for variable \"A\". make sure that you use the same proper damping factors for omega both in the SchurComplementFactory and in the SIMPLE smoother!");
123 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
124 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
128 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 FactoryMonitor m(*
this,
"Setup blocked Uzawa Smoother", currentLevel);
138 this->GetOStream(
Warnings0) <<
"MueLu::UzawaSmoother::Setup(): Setup() has already been called";
141 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
143 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
144 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null,
Exceptions::BadCast,
"MueLu::UzawaSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
147 rangeMapExtractor_ = bA->getRangeMapExtractor();
148 domainMapExtractor_ = bA->getDomainMapExtractor();
151 F_ = bA->getMatrix(0, 0);
152 G_ = bA->getMatrix(0, 1);
153 D_ = bA->getMatrix(1, 0);
154 Z_ = bA->getMatrix(1, 1);
159 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
161 velPredictSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",velpredictFactManager->GetFactory(
"Smoother").get());
164 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
166 schurCompSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother", schurFactManager->GetFactory(
"Smoother").get());
172 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
179 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
181 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
182 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
190 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
191 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
194 bool bCopyResultX =
false;
195 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
197 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
198 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
200 if(bX.is_null() ==
true) {
201 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
206 if(bB.is_null() ==
true) {
207 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
212 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
213 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
216 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
217 if(rbA.is_null() ==
false) {
219 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
222 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
224 Teuchos::RCP<MultiVector> test =
225 buildReorderedBlockedMultiVector(brm, bX);
228 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
230 Teuchos::RCP<const MultiVector> test =
231 buildReorderedBlockedMultiVector(brm, bB);
240 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
241 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
242 Teuchos::RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
243 Teuchos::RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
246 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
247 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
248 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0,bDomainThyraMode);
249 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1,bDomainThyraMode);
254 residual->update(one,*rcpB,zero);
255 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
259 bxtilde->putScalar(zero);
260 velPredictSmoo_->Apply(*xtilde1,*r1);
264 RCP<MultiVector> schurCompRHS = MultiVectorFactory::Build(r2->getMap(),rcpB->getNumVectors());
265 D_->apply(*xtilde1,*schurCompRHS);
266 schurCompRHS->update(one,*r2,-one);
270 schurCompSmoo_->Apply(*xtilde2,*schurCompRHS);
272 rcpX->update(omega,*bxtilde,one);
275 if (bCopyResultX ==
true) {
276 RCP<MultiVector> Xmerged = bX->Merge();
277 X.update(one, *Xmerged, zero);
281 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
282 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
287 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
289 std::ostringstream out;
291 out <<
"{type = " << type_ <<
"}";
295 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
300 out0 <<
"Prec. type: " << type_ << std::endl;
303 if (verbLevel &
Debug) {
308 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
311 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)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
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.
Block triangle Uzawa smoother for 2x2 block matrices.
RCP< SmootherPrototype > Copy() const
void Setup(Level ¤tLevel)
Setup routine.
RCP< const ParameterList > GetValidParameterList() const
Input.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager at a specific position.
void Apply(MultiVector &X, MultiVector const &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
std::string description() const
Return a simple one-line description of this object.
void SetSchurCompFactoryManager(RCP< FactoryManager > FactManager)
Set factory manager for internal SchurComplement handling.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void DeclareInput(Level ¤tLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void SetVelocityPredictionFactoryManager(RCP< FactoryManager > FactManager)
Set factory manager for internal velocity prediction.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.