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>
175 #if 0 // def HAVE_MUELU_DEBUG 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();
Important warning messages (one line)
Exception indicating invalid cast attempted.
virtual ~BlockedGaussSeidelSmoother()
Destructor.
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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.
MueLu::DefaultScalar Scalar
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...
virtual const Teuchos::ParameterList & GetParameterList() const
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.
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.
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()'.
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.