53 #ifndef MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_ 54 #define MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 61 #include <Xpetra_Matrix.hpp> 62 #include <Xpetra_BlockedCrsMatrix.hpp> 63 #include <Xpetra_MultiVectorFactory.hpp> 67 #include "MueLu_Utilities.hpp" 69 #include "MueLu_HierarchyUtils.hpp" 74 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
76 : type_(
"blocked GaussSeidel"), A_(
Teuchos::null)
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 validParamList->
set< Scalar > (
"Damping factor", 1.0,
"Damping/Scaling factor in BGS");
90 validParamList->
set< LocalOrdinal > (
"Sweeps", 1,
"Number of BGS sweeps (default = 1)");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 size_t myPos = Teuchos::as<size_t>(pos);
101 if (myPos < FactManager_.size()) {
103 FactManager_.at(myPos) = FactManager;
104 }
else if( myPos == FactManager_.size()) {
106 FactManager_.push_back(FactManager);
109 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
112 FactManager_.push_back(FactManager);
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").
get());
123 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
124 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
128 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
131 currentLevel.
DeclareInput(
"A",(*it)->GetFactory(
"A").get());
137 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
143 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
147 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
152 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.");
153 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.");
156 rangeMapExtractor_ = bA->getRangeMapExtractor();
157 domainMapExtractor_ = bA->getDomainMapExtractor();
162 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
163 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
168 Inverse_.push_back(Smoo);
172 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
178 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
183 #ifdef HAVE_MUELU_DEBUG 184 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.");
185 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.");
197 LocalOrdinal nSweeps = pL.
get<LocalOrdinal>(
"Sweeps");
198 Scalar omega = pL.
get<Scalar>(
"Damping factor");
201 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
204 for(
size_t i = 0; i<Inverse_.size(); i++) {
208 residual->update(1.0,B,0.0);
213 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
214 bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
220 Inverse_.at(i)->Apply(*tXi, *ri,
false);
223 Xi->update(omega,*tXi,1.0);
226 domainMapExtractor_->InsertVector(Xi, i, rcpX, bDomainThyraMode);
233 RCP<MultiVector> res = MultiVectorFactory::Build(B.getMap(), B.getNumVectors(),
true);
246 LocalOrdinal nSweeps = pL.
get<LocalOrdinal>(
"Sweeps");
247 Scalar omega = pL.
get<Scalar>(
"Damping factor");
250 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
253 for(
size_t i = 0; i<Inverse_.size(); i++) {
257 residual->update(1.0,*bB,0.0);
262 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
263 bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
270 Inverse_.at(i)->Apply(*tXi, *ri,
false);
273 Xi->update(omega,*tXi,1.0);
277 domainMapExtractor_->InsertVector(Xi, i, bX, bDomainThyraMode);
282 for(
size_t r = 0; r < domainMapExtractor_->NumMaps(); ++r) {
283 bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[r] ==
false);
284 domainMapExtractor_->InsertVector(bX->getMultiVector(r,bDomainThyraMode),r,rcpX,bDomainThyraMode);
289 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
294 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
296 std::ostringstream out;
298 out <<
"{type = " << type_ <<
"}";
302 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
308 LocalOrdinal nSweeps = pL.
get<LocalOrdinal>(
"Sweeps");
309 Scalar omega = pL.
get<Scalar>(
"Damping factor");
312 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
315 if (verbLevel &
Debug) {
Important warning messages (one line)
Exception indicating invalid cast attempted.
virtual ~BlockedGaussSeidelSmoother()
Destructor.
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)
block Gauss-Seidel method for blocked matrices
Print additional debugging information.
Namespace for MueLu classes and methods.
std::string description() const
Return a simple one-line description of this object.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void DeclareInput(Level ¤tLevel) const
Input.
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...
RCP< SmootherPrototype > Copy() const
RCP< const ParameterList > GetValidParameterList() const
Input.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
BlockedGaussSeidelSmoother()
Constructor.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual const Teuchos::ParameterList & GetParameterList() const
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 print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
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)