53#ifndef MUELU_INDEFBLOCKEDDIAGONALSMOOTHER_DEF_HPP_
54#define MUELU_INDEFBLOCKEDDIAGONALSMOOTHER_DEF_HPP_
56#include "Teuchos_ArrayViewDecl.hpp"
57#include "Teuchos_ScalarTraits.hpp"
61#include <Xpetra_Matrix.hpp>
62#include <Xpetra_CrsMatrixWrap.hpp>
63#include <Xpetra_BlockedCrsMatrix.hpp>
64#include <Xpetra_MultiVectorFactory.hpp>
65#include <Xpetra_VectorFactory.hpp>
66#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
70#include "MueLu_Utilities.hpp"
72#include "MueLu_HierarchyUtils.hpp"
74#include "MueLu_SubBlockAFactory.hpp"
77#include "MueLu_SchurComplementFactory.hpp"
78#include "MueLu_DirectSolver.hpp"
79#include "MueLu_SmootherFactory.hpp"
80#include "MueLu_FactoryManager.hpp"
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 : type_(
"IndefiniteBlockDiagonalSmoother"), A_(
Teuchos::null)
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 RCP<ParameterList> validParamList = rcp(
new ParameterList());
97 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A (must be a 2x2 block matrix)");
98 validParamList->set<
Scalar > (
"Damping factor", 1.0,
"Damping/Scaling factor");
99 validParamList->set<
LocalOrdinal > (
"Sweeps", 1,
"Number of SIMPLE sweeps (default = 1)");
102 return validParamList;
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
108 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::IndefBlockedDiagonalSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
110 size_t myPos = Teuchos::as<size_t>(pos);
112 if (myPos < FactManager_.size()) {
114 FactManager_.at(myPos) = FactManager;
115 }
else if( myPos == FactManager_.size()) {
117 FactManager_.push_back(FactManager);
119 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
120 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
123 FactManager_.push_back(FactManager);
127 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").get());
131 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.size() != 2,
Exceptions::RuntimeError,
"MueLu::IndefBlockedDiagonalSmoother::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\"!");
134 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
135 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
139 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
144 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
146 FactoryMonitor m(*
this,
"Setup for indefinite blocked diagonal smoother", currentLevel);
149 this->GetOStream(
Warnings0) <<
"MueLu::IndefBlockedDiagonalSmoother::Setup(): Setup() has already been called";
152 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
154 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
155 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null,
Exceptions::BadCast,
"MueLu::IndefBlockedDiagonalSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
158 rangeMapExtractor_ = bA->getRangeMapExtractor();
159 domainMapExtractor_ = bA->getDomainMapExtractor();
162 Teuchos::RCP<Matrix> A00 = bA->getMatrix(0, 0);
163 Teuchos::RCP<Matrix> A01 = bA->getMatrix(0, 1);
164 Teuchos::RCP<Matrix> A10 = bA->getMatrix(1, 0);
165 Teuchos::RCP<Matrix> A11 = bA->getMatrix(1, 1);
199 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
201 velPredictSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",velpredictFactManager->GetFactory(
"Smoother").get());
204 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
206 schurCompSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother", schurFactManager->GetFactory(
"Smoother").get());
211 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
216 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
218 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
225 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
226 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
229 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
230 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
233 bool bCopyResultX =
false;
234 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
236 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
237 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
239 if(bX.is_null() ==
true) {
240 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
245 if(bB.is_null() ==
true) {
246 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
251 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
252 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
255 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
256 if(rbA.is_null() ==
false) {
258 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
261 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
263 Teuchos::RCP<MultiVector> test =
264 buildReorderedBlockedMultiVector(brm, bX);
267 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
269 Teuchos::RCP<const MultiVector> test =
270 buildReorderedBlockedMultiVector(brm, bB);
279 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
280 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
281 Teuchos::RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
282 Teuchos::RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
285 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
286 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
287 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0,bDomainThyraMode);
288 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1,bDomainThyraMode);
293 residual->update(one,*rcpB,zero);
294 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
304 bxtilde->putScalar(zero);
305 velPredictSmoo_->Apply(*xtilde1,*r1);
309 schurCompSmoo_->Apply(*xtilde2,*r2);
312 rcpX->update(omega,*bxtilde,one);
314 Teuchos::RCP<MultiVector> x1 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
315 Teuchos::RCP<MultiVector> x2 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
319 x1->update(omega,*xtilde1,one);
320 x2->update(omega,*xtilde2,one);
323 domainMapExtractor_->InsertVector(x1, 0, rcpX, bDomainThyraMode);
324 domainMapExtractor_->InsertVector(x2, 1, rcpX, bDomainThyraMode);
328 if (bCopyResultX ==
true) {
329 RCP<MultiVector> Xmerged = bX->Merge();
330 X.update(one, *Xmerged, zero);
334 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
335 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
340 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
342 std::ostringstream out;
344 out <<
"{type = " << type_ <<
"}";
348 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
353 out0 <<
"Prec. type: " << type_ << std::endl;
356 if (verbLevel &
Debug) {
360 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
363 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.
Cheap Blocked diagonal smoother for indefinite 2x2 block matrices.
IndefBlockedDiagonalSmoother()
Constructor.
virtual ~IndefBlockedDiagonalSmoother()
Destructor.
std::string description() const
Return a simple one-line description of this object.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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 DeclareInput(Level ¤tLevel) const
Input.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
RCP< SmootherPrototype > Copy() const
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
RCP< const ParameterList > GetValidParameterList() const
Input.
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.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.