47#ifndef MUELU_BLOCKEDJACOBISMOOTHER_DEF_HPP_
48#define MUELU_BLOCKEDJACOBISMOOTHER_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 Jacobi"), 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 Jacobi");
87 validParamList->set<
LocalOrdinal>(
"Sweeps", 1,
"Number of 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::BlockedJacobiSmoother::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 Jacobi 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::BlockedJacobiSmoother::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::BlockedJacobiSmoother::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) {
199 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
202 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
203 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
206 bool bCopyResultX =
false;
207 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
209 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
210 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
212 if(bX.is_null() ==
true) {
213 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
218 if(bB.is_null() ==
true) {
219 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
224 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
225 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
228 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
229 if(rbA.is_null() ==
false) {
231 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
234 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
236 Teuchos::RCP<MultiVector> test =
237 buildReorderedBlockedMultiVector(brm, bX);
240 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
242 Teuchos::RCP<const MultiVector> test =
243 buildReorderedBlockedMultiVector(brm, bB);
288 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
289 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
298 residual->update(1.0,*rcpB,0.0);
299 if(InitialGuessIsZero ==
false || run > 0){
300 bA->apply(*rcpX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
303 for(
size_t i = 0; i<Inverse_.size(); i++) {
306 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
307 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
308 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
309 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
310 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
313 Inverse_.at(i)->Apply(*tXi, *ri,
false);
316 if( InitialGuessIsZero && run == 0 ){
317 Xi->update(omega,*tXi,0.0);
319 Xi->update(omega,*tXi,1.0);
325 if (bCopyResultX ==
true) {
326 RCP<MultiVector> Xmerged = bX->Merge();
327 X.update(one, *Xmerged, zero);
331 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
336 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
338 std::ostringstream out;
340 out <<
"{type = " << type_ <<
"}";
344 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
354 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
357 if (verbLevel &
Debug) {
362 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
365 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 Jacobi method for blocked matrices
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver.
RCP< const ParameterList > GetValidParameterList() const
Input.
void Setup(Level ¤tLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
std::string description() const
Return a simple one-line description of this object.
void DeclareInput(Level ¤tLevel) const
Input.
RCP< SmootherPrototype > Copy() const
BlockedJacobiSmoother()
Constructor.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
virtual ~BlockedJacobiSmoother()
Destructor.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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.