53 #ifndef MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_ 54 #define MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 57 #include "Teuchos_ScalarTraits.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>
86 RCP<ParameterList> validParamList = rcp(
new ParameterList());
88 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
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>
97 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
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);
108 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
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>
141 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
143 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
147 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
148 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
149 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
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();
159 Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
162 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
163 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
167 RCP<const SmootherBase> Smoo = currentLevel.Get< RCP<SmootherBase> >(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
168 Inverse_.push_back(Smoo);
171 RCP<Matrix> Aii = currentLevel.Get< RCP<Matrix> >(
"A",(*it)->GetFactory(
"A").get());
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.");
189 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > residual = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
190 RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tempres = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
191 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
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);
210 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
213 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
214 bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
215 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
216 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
217 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
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);
234 RCP<BlockedMultiVector> residual = Teuchos::rcp(
new BlockedMultiVector(rangeMapExtractor_,res));
237 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
238 RCP<BlockedMultiVector> bX = Teuchos::rcp(
new BlockedMultiVector(domainMapExtractor_,rcpX));
241 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
242 RCP<const BlockedMultiVector> bB = Teuchos::rcp(
new const BlockedMultiVector(rangeMapExtractor_,rcpB));
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);
259 A_->apply(*bX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
262 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
263 bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] ==
false);
265 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(bX, i, bDomainThyraMode);
266 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
267 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode,
true);
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.
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.
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.
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.