53 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 54 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 68 #include "MueLu_Utilities.hpp" 70 #include "MueLu_HierarchyUtils.hpp" 74 #include "MueLu_SchurComplementFactory.hpp" 75 #include "MueLu_DirectSolver.hpp" 76 #include "MueLu_SmootherFactory.hpp" 77 #include "MueLu_FactoryManager.hpp" 81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 : type_(
"Braess Sarazin"), A_(null)
94 SchurFact->SetFactory(
"A", this->
GetFactory(
"A"));
100 smoProtoSC->SetFactory(
"A", SchurFact);
105 FactManager->SetFactory(
"A", SchurFact);
106 FactManager->SetFactory(
"Smoother", SmooSCFact);
107 FactManager->SetIgnoreUserData(
true);
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 FactManager_ = FactManager;
123 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 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 " 132 "as approximation of A00 (and A00^{-1})");
133 validParamList->
set<
SC> (
"Damping factor", one,
"Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
134 validParamList->
set<
LO> (
"Sweeps", 1,
"Number of BraessSarazin sweeps (default = 1)");
136 validParamList->
set<
bool> (
"q2q1 mode",
false,
"Run in the mode matching Q2Q1 matlab code");
138 return validParamList;
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 this->Input(currentLevel,
"A");
146 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! " 147 "Introduce a FactoryManager for the SchurComplement equation.");
154 currentLevel.
DeclareInput(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
157 currentLevel.
DeclareInput(
"A", FactManager_->GetFactory(
"A").get());
166 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
171 this->GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
174 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
177 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
180 rangeMapExtractor_ = bA->getRangeMapExtractor();
181 domainMapExtractor_ = bA->getDomainMapExtractor();
184 A00_ = bA->getMatrix(0,0);
185 A01_ = bA->getMatrix(0,1);
186 A10_ = bA->getMatrix(1,0);
187 A11_ = bA->getMatrix(1,1);
190 SC omega = pL.
get<
SC>(
"Damping factor");
194 D_ = VectorFactory::Build(A00_->getRowMap());
197 if (pL.
get<
bool>(
"lumping") ==
false)
205 for (
GO row = 0; row < Ddata.
size(); row++)
206 Ddata[row] = one / (diag[row]*omega);*/
210 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
211 if (pL.
get<
bool>(
"lumping") ==
false) {
212 A00_->getLocalDiagCopy(*diagFVector);
216 diagFVector->scale(omega);
224 smoo_ = currentLevel.Get<
RCP<SmootherBase> >(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
225 S_ = currentLevel.Get<
RCP<Matrix> > (
"A", FactManager_->GetFactory(
"A").get());
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
234 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
236 #ifdef HAVE_MUELU_DEBUG 237 TEUCHOS_TEST_FOR_EXCEPTION(A_->getRangeMap()->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.");
238 TEUCHOS_TEST_FOR_EXCEPTION(A_->getDomainMap()->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.");
245 bool bA11ThyraSpecialTreatment =
false;
246 if (bA11 != Teuchos::null) {
247 if(bA11->Rows() == 1 && bA11->Cols() == 1 && rangeMapExtractor_->getThyraMode() ==
true) bA11ThyraSpecialTreatment =
true;
254 RCP<MultiVector> deltaX0 = MultiVectorFactory::Build(A00_->getRowMap(), X.getNumVectors());
255 RCP<MultiVector> deltaX1 = MultiVectorFactory::Build(A10_->getRowMap(), X.getNumVectors());
256 RCP<MultiVector> Rtmp = MultiVectorFactory::Build(A10_->getRowMap(), B.getNumVectors());
259 SC one = STS::one(), zero = STS::zero();
263 LO nSweeps = pL.
get<
LO>(
"Sweeps");
266 if (InitialGuessIsZero) {
267 R = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
268 R->update(one, B, zero);
274 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
275 S_->getLocalDiagCopy(*diagSVector);
278 for (
LO run = 0; run < nSweeps; ++run) {
289 deltaX0->putScalar(zero);
290 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
291 A10_->apply(*deltaX0, *Rtmp);
292 Rtmp->update(one, *R1, -one);
294 if (!pL.
get<
bool>(
"q2q1 mode")) {
295 deltaX1->putScalar(zero);
299 for (
GO row = 0; row < deltaX1data.
size(); row++)
300 deltaX1data[row] = 1.1*Rtmpdata[row] / Sdiag[row];
305 if(bA11ThyraSpecialTreatment ==
true) {
306 RCP<MultiVector> deltaX1_thyra = domainMapExtractor_->getVector(1, X.getNumVectors(),
true);
307 RCP<MultiVector> Rtmp_thyra = rangeMapExtractor_->getVector(1, B.getNumVectors(),
true);
309 for(
size_t k=0; k < Rtmp->getNumVectors(); k++) {
312 for(
size_t i=0; i < Rtmp->getLocalLength(); i++) {
313 thyraVecData[i] = xpetraVecData[i];
317 smoo_->Apply(*deltaX1_thyra,*Rtmp_thyra);
319 for(
size_t k=0; k < deltaX1_thyra->getNumVectors(); k++) {
322 for(
size_t i=0; i < deltaX1_thyra->getLocalLength(); i++) {
323 xpetraVecData[i] = thyraVecData[i];
329 smoo_->Apply(*deltaX1,*Rtmp);
333 deltaX0->putScalar(zero);
334 A01_->apply(*deltaX1, *deltaX0);
335 deltaX0->update(one, *R0, -one);
337 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
340 X0->update(one, *deltaX0, one);
341 X1->update(one, *deltaX1, one);
343 domainMapExtractor_->InsertVector(X0, 0, rcpX);
344 domainMapExtractor_->InsertVector(X1, 1, rcpX);
351 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
357 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
360 std::ostringstream out;
362 out <<
"{type = " << type_ <<
"}";
366 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
371 out0 <<
"Prec. type: " << type_ << std::endl;
374 if (verbLevel &
Debug) {
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()
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
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)