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();
Important warning messages (one line)
Exception indicating invalid cast attempted.
RCP< SmootherPrototype > Copy() const
BraessSarazin smoother for 2x2 block matrices.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
bool IsSetup() const
Get the state of a smoother prototype.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
Namespace for MueLu classes and methods.
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)
std::string description() const
Return a simple one-line description of this object.
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)
Class that holds all level-specific information.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
virtual const Teuchos::ParameterList & GetParameterList() const
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.
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.
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()'.
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.