47#ifndef MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_
48#define MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_
50#include "Teuchos_ArrayViewDecl.hpp"
51#include "Teuchos_ScalarTraits.hpp"
55#include <Xpetra_BlockReorderManager.hpp>
56#include <Xpetra_Matrix.hpp>
57#include <Xpetra_BlockedCrsMatrix.hpp>
58#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
59#include <Xpetra_ReorderedBlockedMultiVector.hpp>
60#include <Xpetra_MultiVectorFactory.hpp>
64#include "MueLu_Utilities.hpp"
66#include "MueLu_HierarchyUtils.hpp"
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 : type_(
"blocked GaussSeidel"), A_(
Teuchos::null)
78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 RCP<ParameterList> validParamList = rcp(
new ParameterList());
85 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
86 validParamList->set<
Scalar > (
"Damping factor", 1.0,
"Damping/Scaling factor in BGS");
87 validParamList->set<
LocalOrdinal > (
"Sweeps", 1,
"Number of BGS sweeps (default = 1)");
89 return validParamList;
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
96 size_t myPos = Teuchos::as<size_t>(pos);
98 if (myPos < FactManager_.size()) {
100 FactManager_.at(myPos) = FactManager;
101 }
else if(myPos == FactManager_.size()) {
103 FactManager_.push_back(FactManager);
105 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
106 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
109 FactManager_.push_back(FactManager);
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").get());
120 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
121 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
125 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
128 currentLevel.
DeclareInput(
"A",(*it)->GetFactory(
"A").get());
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
137 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
141 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
142 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
143 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
146 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.");
147 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.");
150 rangeMapExtractor_ = bA->getRangeMapExtractor();
151 domainMapExtractor_ = bA->getDomainMapExtractor();
154 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
155 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
159 RCP<const SmootherBase> Smoo = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
160 Inverse_.push_back(Smoo);
163 RCP<Matrix> Aii = currentLevel.
Get< RCP<Matrix> >(
"A",(*it)->GetFactory(
"A").get());
164 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
170 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
177 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
178 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
179 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
180 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpDebugB);
182 if(rcpBDebugB.is_null() ==
false) {
189 if(rcpBDebugX.is_null() ==
false) {
197 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
200 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
201 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
204 bool bCopyResultX =
false;
205 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
207 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
208 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
210 if(bX.is_null() ==
true) {
211 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
216 if(bB.is_null() ==
true) {
217 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
222 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
223 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
226 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
227 if(rbA.is_null() ==
false) {
229 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
232 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
234 Teuchos::RCP<MultiVector> test =
235 buildReorderedBlockedMultiVector(brm, bX);
238 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
240 Teuchos::RCP<const MultiVector> test =
241 buildReorderedBlockedMultiVector(brm, bB);
286 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
287 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
296 if( InitialGuessIsZero==
true )
297 rcpX->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
304 for(
size_t i = 0; i<Inverse_.size(); i++) {
308 residual->update(1.0,*rcpB,0.0);
309 if(InitialGuessIsZero ==
false || i > 0 || run > 0)
310 bA->bgs_apply(*rcpX, *residual, i, Teuchos::NO_TRANS, -1.0, 1.0);
313 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
314 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
315 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
316 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
317 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
320 Inverse_.at(i)->Apply(*tXi, *ri,
false);
323 Xi->update(omega,*tXi,1.0);
328 if (bCopyResultX ==
true) {
329 RCP<MultiVector> Xmerged = bX->Merge();
330 X.update(one, *Xmerged, zero);
334 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
339 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
341 std::ostringstream out;
343 out <<
"{type = " << type_ <<
"}";
347 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
357 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
360 if (verbLevel &
Debug) {
364 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
367 return Teuchos::OrdinalTraits<size_t>::invalid();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
block Gauss-Seidel method for blocked matrices
void DeclareInput(Level ¤tLevel) const
Input.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
std::string description() const
Return a simple one-line description of this object.
virtual ~BlockedGaussSeidelSmoother()
Destructor.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.
RCP< SmootherPrototype > Copy() const
RCP< const ParameterList > GetValidParameterList() const
Input.
BlockedGaussSeidelSmoother()
Constructor.
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 verbLevel to an FancyOStream object out.
void Setup(Level ¤tLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual const Teuchos::ParameterList & GetParameterList() const =0
An exception safe way to call the method 'Level::SetFactoryManager()'.
bool IsSetup() const
Get the state of a smoother prototype.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.