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>
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) {
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());
293 LocalOrdinal nSweeps = pL.get<LocalOrdinal>(
"Sweeps");
294 Scalar omega = pL.get<Scalar>(
"Damping factor");
297 for (LocalOrdinal run = 0; run < nSweeps; ++run) {
299 for(
size_t i = 0; i<Inverse_.size(); i++) {
303 residual->update(1.0,*rcpB,0.0);
304 if(InitialGuessIsZero ==
false || i > 0 || run > 0)
305 bA->bgs_apply(*rcpX, *residual, i, Teuchos::NO_TRANS, -1.0, 1.0);
308 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
309 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
310 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
311 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
312 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
315 Inverse_.at(i)->Apply(*tXi, *ri,
false);
318 Xi->update(omega,*tXi,1.0);
323 if (bCopyResultX ==
true) {
324 RCP<MultiVector> Xmerged = bX->Merge();
325 X.update(one, *Xmerged, zero);
329 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
334 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
336 std::ostringstream out;
338 out <<
"{type = " << type_ <<
"}";
342 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
348 LocalOrdinal nSweeps = pL.get<LocalOrdinal>(
"Sweeps");
349 Scalar omega = pL.get<Scalar>(
"Damping factor");
352 out0 <<
"Prec. type: " << type_ <<
" Sweeps: " << nSweeps <<
" damping: " << omega << std::endl;
355 if (verbLevel &
Debug) {
360 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
363 return Teuchos::OrdinalTraits<size_t>::invalid();
Important warning messages (one line)
Exception indicating invalid cast attempted.
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.
Print additional debugging information.
Namespace for MueLu classes and methods.
void Setup(Level ¤tLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
Class that holds all level-specific information.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< SmootherPrototype > Copy() const
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual const Teuchos::ParameterList & GetParameterList() const
RCP< const ParameterList > GetValidParameterList() const
Input.
std::string description() const
Return a simple one-line description of this object.
block Jacobi method for blocked matrices
BlockedJacobiSmoother()
Constructor.
void DeclareInput(Level ¤tLevel) const
Input.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver.
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 print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
virtual ~BlockedJacobiSmoother()
Destructor.
virtual std::string description() const
Return a simple one-line description of this object.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.