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_BlockReorderManager.hpp> 62 #include <Xpetra_Matrix.hpp> 63 #include <Xpetra_BlockedCrsMatrix.hpp> 64 #include <Xpetra_ReorderedBlockedCrsMatrix.hpp> 65 #include <Xpetra_ReorderedBlockedMultiVector.hpp> 66 #include <Xpetra_MultiVectorFactory.hpp> 70 #include "MueLu_Utilities.hpp" 72 #include "MueLu_HierarchyUtils.hpp" 77 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 : type_(
"blocked GaussSeidel"), A_(
Teuchos::null)
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
89 RCP<ParameterList> validParamList = rcp(
new ParameterList());
91 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
92 validParamList->set< Scalar > (
"Damping factor", 1.0,
"Damping/Scaling factor in BGS");
93 validParamList->set< LocalOrdinal > (
"Sweeps", 1,
"Number of BGS sweeps (default = 1)");
95 return validParamList;
98 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
102 size_t myPos = Teuchos::as<size_t>(pos);
104 if (myPos < FactManager_.size()) {
106 FactManager_.at(myPos) = FactManager;
107 }
else if( myPos == FactManager_.size()) {
109 FactManager_.push_back(FactManager);
111 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
112 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
115 FactManager_.push_back(FactManager);
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
123 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").
get());
126 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
127 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
131 currentLevel.
DeclareInput(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
134 currentLevel.
DeclareInput(
"A",(*it)->GetFactory(
"A").get());
140 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
146 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
150 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
151 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
152 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
155 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.");
156 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.");
159 rangeMapExtractor_ = bA->getRangeMapExtractor();
160 domainMapExtractor_ = bA->getDomainMapExtractor();
163 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
164 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
168 RCP<const SmootherBase> Smoo = currentLevel.Get< RCP<SmootherBase> >(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
169 Inverse_.push_back(Smoo);
172 RCP<Matrix> Aii = currentLevel.Get< RCP<Matrix> >(
"A",(*it)->GetFactory(
"A").get());
173 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
179 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
184 #if 0 // def HAVE_MUELU_DEBUG 186 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
187 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
188 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
189 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpDebugB);
191 if(rcpBDebugB.is_null() ==
false) {
198 if(rcpBDebugX.is_null() ==
false) {
208 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
211 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
212 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
215 bool bCopyResultX =
false;
216 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
218 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
219 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
221 if(bX.is_null() ==
true) {
222 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
227 if(bB.is_null() ==
true) {
228 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
233 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
234 bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
237 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
238 if(rbA.is_null() ==
false) {
240 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
243 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
245 Teuchos::RCP<MultiVector> test =
246 buildReorderedBlockedMultiVector(brm, bX);
249 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
251 Teuchos::RCP<const MultiVector> test =
252 buildReorderedBlockedMultiVector(brm, bB);
297 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
298 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
302 LocalOrdinal nSweeps = pL.get<LocalOrdinal>(
"Sweeps");
303 Scalar omega = pL.get<Scalar>(
"Damping factor");
306 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
309 for(
size_t i = 0; i<Inverse_.size(); i++) {
313 residual->update(1.0,*rcpB,0.0);
314 if(InitialGuessIsZero ==
false || i > 0 || run > 0)
315 bA->bgs_apply(*rcpX, *residual, i, Teuchos::NO_TRANS, -1.0, 1.0);
319 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
320 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
321 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
322 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
323 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
326 Inverse_.at(i)->Apply(*tXi, *ri,
false);
329 Xi->update(omega,*tXi,1.0);
332 domainMapExtractor_->InsertVector(Xi, i, rcpX, bDomainThyraMode);
336 if (bCopyResultX ==
true) {
337 RCP<MultiVector> Xmerged = bX->Merge();
338 X.update(one, *Xmerged, zero);
342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
347 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
349 std::ostringstream out;
351 out <<
"{type = " << type_ <<
"}";
355 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
361 LocalOrdinal nSweeps = pL.get<LocalOrdinal>(
"Sweeps");
362 Scalar omega = pL.get<Scalar>(
"Damping factor");
365 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
368 if (verbLevel &
Debug) {
372 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
375 return Teuchos::OrdinalTraits<size_t>::invalid();
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.
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.
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.
#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.