47 #ifndef MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_
48 #define MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_
50 #include "Teuchos_ArrayViewDecl.hpp"
60 #include <Xpetra_MultiVectorFactory.hpp>
65 #include "MueLu_HierarchyUtils.hpp"
70 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
72 : type_(
"blocked GaussSeidel")
77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 validParamList->
set<
Scalar>(
"Damping factor", 1.0,
"Damping/Scaling factor in BGS");
86 validParamList->
set<
LocalOrdinal>(
"Sweeps", 1,
"Number of BGS sweeps (default = 1)");
88 return validParamList;
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
95 size_t myPos = Teuchos::as<size_t>(pos);
97 if (myPos < FactManager_.size()) {
99 FactManager_.at(myPos) = FactManager;
100 }
else if (myPos == FactManager_.size()) {
102 FactManager_.push_back(FactManager);
105 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
108 FactManager_.push_back(FactManager);
112 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 currentLevel.
DeclareInput(
"A", this->GetFactory(
"A").
get());
119 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
120 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
124 currentLevel.
DeclareInput(
"PreSmoother", (*it)->GetFactory(
"Smoother").get());
127 currentLevel.
DeclareInput(
"A", (*it)->GetFactory(
"A").get());
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
139 A_ = Factory::Get<RCP<Matrix> >(currentLevel,
"A");
144 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Setup: number of block rows of A is " << bA->Rows() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
145 TEUCHOS_TEST_FOR_EXCEPTION(bA->Cols() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Setup: number of block cols of A is " << bA->Cols() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
148 rangeMapExtractor_ = bA->getRangeMapExtractor();
149 domainMapExtractor_ = bA->getDomainMapExtractor();
152 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
153 for (it = FactManager_.begin(); it != FactManager_.end(); ++it) {
158 Inverse_.push_back(Smoo);
162 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii) != Teuchos::null);
168 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
172 #if 0 // def HAVE_MUELU_DEBUG
179 if(rcpBDebugB.is_null() ==
false) {
186 if(rcpBDebugX.
is_null() ==
false) {
201 bool bCopyResultX =
false;
202 bool bReorderX =
false;
203 bool bReorderB =
false;
213 if (rbA.is_null() ==
false) {
218 if (bX.is_null() ==
true) {
220 rcpX.
swap(vectorWithBlockedMap);
226 if (bB.is_null() ==
true) {
228 rcpB.
swap(vectorWithBlockedMap);
233 if (bX.is_null() ==
true) {
235 rcpX.
swap(vectorWithBlockedMap);
239 if (bB.is_null() ==
true) {
241 rcpB.
swap(vectorWithBlockedMap);
246 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
247 bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
250 if (rbA.is_null() ==
false) {
254 if (bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) {
257 rcpX.
swap(reorderedBlockedVector);
260 if (bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) {
263 rcpB.
swap(reorderedBlockedVector);
270 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
278 if (InitialGuessIsZero ==
true)
285 for (
size_t i = 0; i < Inverse_.size(); i++) {
288 residual->update(1.0, *rcpB, 0.0);
289 if (InitialGuessIsZero ==
false || i > 0 || run > 0)
293 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
294 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
300 Inverse_.at(i)->Apply(*tXi, *ri,
false);
303 Xi->update(omega, *tXi, 1.0);
307 if (bCopyResultX ==
true) {
309 X.update(one, *Xmerged, zero);
313 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
318 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
320 std::ostringstream out;
322 out <<
"{type = " << type_ <<
"}";
326 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
336 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
339 if (verbLevel &
Debug) {
343 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
Important warning messages (one line)
virtual const Teuchos::ParameterList & GetParameterList() const
Exception indicating invalid cast attempted.
virtual ~BlockedGaussSeidelSmoother()
Destructor.
MueLu::DefaultLocalOrdinal LocalOrdinal
T & get(const std::string &name, T def_value)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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)
block Gauss-Seidel method for blocked matrices
void swap(RCP< T > &r_ptr)
Print additional debugging information.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level verbLevel to an FancyOStream object out.
std::string description() const
Return a simple one-line description of this object.
void DeclareInput(Level ¤tLevel) const
Input.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
void Setup(Level ¤tLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
bool IsSetup() const
Get the state of a smoother prototype.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
BlockedGaussSeidelSmoother()
Constructor.
RCP< SmootherPrototype > Copy() const
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
RCP< const ParameterList > GetValidParameterList() const
Input.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver. Solves the linear system AX=B using the constructed solver.
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
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 AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.
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.
std::string toString(const T &t)