47 #ifndef MUELU_BLOCKEDJACOBISMOOTHER_DEF_HPP_
48 #define MUELU_BLOCKEDJACOBISMOOTHER_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 Jacobi")
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 Jacobi");
86 validParamList->
set<
LocalOrdinal>(
"Sweeps", 1,
"Number of 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 Jacobi Smoother", currentLevel);
139 A_ = Factory::Get<RCP<Matrix> >(currentLevel,
"A");
144 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedJacobiSmoother::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::BlockedJacobiSmoother::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) {
202 bool bCopyResultX =
false;
203 bool bReorderX =
false;
204 bool bReorderB =
false;
214 if (rbA.is_null() ==
false) {
219 if (bX.is_null() ==
true) {
221 rcpX.
swap(vectorWithBlockedMap);
227 if (bB.is_null() ==
true) {
229 rcpB.
swap(vectorWithBlockedMap);
234 if (bX.is_null() ==
true) {
236 rcpX.
swap(vectorWithBlockedMap);
240 if (bB.is_null() ==
true) {
242 rcpB.
swap(vectorWithBlockedMap);
247 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
248 bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
251 if (rbA.is_null() ==
false) {
255 if (bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps() || bReorderX) {
258 rcpX.
swap(reorderedBlockedVector);
261 if (bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps() || bReorderB) {
264 rcpB.
swap(reorderedBlockedVector);
271 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
280 residual->update(1.0, *rcpB, 0.0);
281 if (InitialGuessIsZero ==
false || run > 0) {
285 for (
size_t i = 0; i < Inverse_.size(); i++) {
287 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
288 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
294 Inverse_.at(i)->Apply(*tXi, *ri,
false);
297 if (InitialGuessIsZero && run == 0) {
298 Xi->update(omega, *tXi, 0.0);
300 Xi->update(omega, *tXi, 1.0);
305 if (bCopyResultX ==
true) {
307 X.update(one, *Xmerged, zero);
311 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
316 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
318 std::ostringstream out;
320 out <<
"{type = " << type_ <<
"}";
324 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
334 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
337 if (verbLevel &
Debug) {
342 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.
MueLu::DefaultLocalOrdinal LocalOrdinal
RCP< const ParameterList > GetValidParameterList() const
Input.
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)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
std::string description() const
Return a simple one-line description of this object.
void swap(RCP< T > &r_ptr)
Print additional debugging information.
void Setup(Level ¤tLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Class that holds all level-specific information.
void DeclareInput(Level ¤tLevel) const
Input.
bool IsSetup() const
Get the state of a smoother prototype.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
block Jacobi method for blocked matrices
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver.
BlockedJacobiSmoother()
Constructor.
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 DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
virtual ~BlockedJacobiSmoother()
Destructor.
RCP< SmootherPrototype > Copy() const
virtual std::string description() const
Return a simple one-line description of this object.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.
std::string toString(const T &t)