53 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 54 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 69 #include "MueLu_Utilities.hpp" 71 #include "MueLu_HierarchyUtils.hpp" 75 #include "MueLu_SchurComplementFactory.hpp" 76 #include "MueLu_DirectSolver.hpp" 77 #include "MueLu_SmootherFactory.hpp" 78 #include "MueLu_FactoryManager.hpp" 82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 : type_(
"Braess Sarazin"), A_(null)
95 SchurFact->SetFactory(
"A", this->
GetFactory(
"A"));
101 smoProtoSC->SetFactory(
"A", SchurFact);
106 FactManager->SetFactory(
"A", SchurFact);
107 FactManager->SetFactory(
"Smoother", SmooSCFact);
108 FactManager->SetIgnoreUserData(
true);
114 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
118 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 FactManager_ = FactManager;
124 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
132 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 " 133 "as approximation of A00 (and A00^{-1})");
134 validParamList->
set<
SC> (
"Damping factor", one,
"Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
135 validParamList->
set<
LO> (
"Sweeps", 1,
"Number of BraessSarazin sweeps (default = 1)");
137 validParamList->
set<
bool> (
"q2q1 mode",
false,
"Run in the mode matching Q2Q1 matlab code");
139 return validParamList;
142 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 this->Input(currentLevel,
"A");
147 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! " 148 "Introduce a FactoryManager for the SchurComplement equation.");
155 currentLevel.
DeclareInput(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
158 currentLevel.
DeclareInput(
"A", FactManager_->GetFactory(
"A").get());
167 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
172 this->GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
175 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
178 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
181 rangeMapExtractor_ = bA->getRangeMapExtractor();
182 domainMapExtractor_ = bA->getDomainMapExtractor();
185 A00_ = bA->getMatrix(0,0);
186 A01_ = bA->getMatrix(0,1);
187 A10_ = bA->getMatrix(1,0);
188 A11_ = bA->getMatrix(1,1);
191 SC omega = pL.
get<
SC>(
"Damping factor");
195 D_ = VectorFactory::Build(A00_->getRowMap());
198 if (pL.
get<
bool>(
"lumping") ==
false)
206 for (
GO row = 0; row < Ddata.
size(); row++)
207 Ddata[row] = one / (diag[row]*omega);*/
211 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
212 if (pL.
get<
bool>(
"lumping") ==
false) {
213 A00_->getLocalDiagCopy(*diagFVector);
217 diagFVector->scale(omega);
222 if(bD.
is_null() ==
false && bD->getBlockedMap()->getNumMaps() == 1) {
223 RCP<Vector> nestedVec = bD->getMultiVector(0,bD->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
232 smoo_ = currentLevel.Get<
RCP<SmootherBase> >(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
233 S_ = currentLevel.Get<
RCP<Matrix> > (
"A", FactManager_->GetFactory(
"A").get());
239 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
242 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
244 #if 0 //def HAVE_MUELU_DEBUG 251 if(rcpBDebugB.is_null() ==
false) {
254 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullRangeMap()->isSameAs(*(B.getMap())) ==
false,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Apply(): The map of RHS vector B is not the same as range map of the blocked operator A. Please check the map of B and A.");
256 if(rcpBDebugX.
is_null() ==
false) {
259 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullDomainMap()->isSameAs(*(X.getMap())) ==
false,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Apply(): The map of the solution vector X is not the same as domain map of the blocked operator A. Please check the map of X and A.");
268 bool bCopyResultX =
false;
274 if(bX.is_null() ==
true) {
280 if(bB.is_null() ==
true) {
286 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
287 bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
291 if(rbA.is_null() ==
false) {
296 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
302 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
313 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
314 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
316 RCP<MultiVector> deltaX = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
321 RCP<MultiVector> Rtmp = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
325 SC one = STS::one(), zero = STS::zero();
329 LO nSweeps = pL.
get<
LO>(
"Sweeps");
332 if (InitialGuessIsZero) {
333 R = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
334 R->update(one, *rcpB, zero);
340 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
341 S_->getLocalDiagCopy(*diagSVector);
343 for (
LO run = 0; run < nSweeps; ++run) {
347 RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0, bRangeThyraMode);
348 RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1, bRangeThyraMode);
351 deltaX0->putScalar(zero);
352 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
353 A10_->apply(*deltaX0, *Rtmp);
354 Rtmp->update(one, *R1, -one);
356 if (!pL.
get<
bool>(
"q2q1 mode")) {
357 deltaX1->putScalar(zero);
360 if(Teuchos::rcp_dynamic_cast<BlockedVector>(diagSVector) == Teuchos::null) {
364 for (
GO row = 0; row < deltaX1data.
size(); row++)
365 deltaX1data[row] = 1.1*Rtmpdata[row] / Sdiag[row];
371 smoo_->Apply(*deltaX1,*Rtmp);
374 deltaX0->putScalar(zero);
375 A01_->apply(*deltaX1, *deltaX0);
376 deltaX0->update(one, *R0, -one);
378 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
380 RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
381 RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
384 X0->update(one, *deltaX0, one);
385 X1->update(one, *deltaX1, one);
387 domainMapExtractor_->InsertVector(X0, 0, rcpX, bDomainThyraMode);
388 domainMapExtractor_->InsertVector(X1, 1, rcpX, bDomainThyraMode);
390 if (run < nSweeps-1) {
396 if (bCopyResultX ==
true) {
398 X.update(one, *Xmerged, zero);
403 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
409 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
412 std::ostringstream out;
414 out <<
"{type = " << type_ <<
"}";
418 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
423 out0 <<
"Prec. type: " << type_ << std::endl;
426 if (verbLevel &
Debug) {
431 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Exception indicating invalid cast attempted.
This class specifies the default factory that should generate some data on a Level if the data does n...
RCP< SmootherPrototype > Copy() const
BraessSarazin smoother for 2x2 block matrices.
bool IsSetup() const
Get the state of a smoother prototype.
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.
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Namespace for MueLu classes and methods.
std::string description() const
Return a simple one-line description of this object.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
virtual ~BraessSarazinSmoother()
Destructor.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual const Teuchos::ParameterList & GetParameterList() const
BraessSarazinSmoother()
Constructor.
Factory for building the Schur Complement for a 2x2 block matrix.
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())
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
void DeclareInput(Level ¤tLevel) const
Input.
void Setup(Level ¤tLevel)
Setup routine.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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)
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.
RCP< const ParameterList > GetValidParameterList() const
Input.
std::string toString(const T &t)