53 #ifndef MUELU_BLOCKEDDIRECTSOLVER_DEF_HPP_ 54 #define MUELU_BLOCKEDDIRECTSOLVER_DEF_HPP_ 56 #include "Teuchos_ArrayViewDecl.hpp" 57 #include "Teuchos_ScalarTraits.hpp" 61 #include <Xpetra_Matrix.hpp> 62 #include <Xpetra_BlockedCrsMatrix.hpp> 63 #include <Xpetra_MultiVectorFactory.hpp> 65 #include "MueLu_BlockedDirectSolver.hpp" 66 #include "MueLu_MergedBlockedMatrixFactory.hpp" 68 #include "MueLu_Utilities.hpp" 70 #include "MueLu_HierarchyUtils.hpp" 72 #include "MueLu_DirectSolver.hpp" 76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 : type_(
"blocked direct solver")
82 Teuchos::ParameterList params;
86 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 RCP<ParameterList> validParamList = rcp(
new ParameterList());
90 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 currentLevel.
DeclareInput(
"A",this->GetFactory(
"A").
get());
103 MergedAFact_->SetFactory(
"A", this->GetFactory(
"A"));
106 s_->SetFactory(
"A",MergedAFact_);
107 s_->DeclareInput(currentLevel);
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
114 FactoryMonitor m(*
this,
"Setup BlockedDirectSolver", currentLevel);
115 if (this->IsSetup() ==
true)
116 this->GetOStream(
Warnings0) <<
"MueLu::BlockedDirectSolver::Setup(): Setup() has already been called";
119 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
120 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
122 "MueLu::BlockedDirectSolver::Build: input matrix A is not of type BlockedCrsMatrix.");
124 s_->Setup(currentLevel);
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
132 "MueLu::BlockedDirectSolver::Apply(): Setup() has not been called");
134 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
135 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
136 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
137 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<
const BlockedMultiVector>(rcpB);
139 #ifdef HAVE_MUELU_DEBUG 140 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
141 if(bB.is_null() ==
false) {
144 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullRangeMap()->isSameAs(*(B.getMap())) ==
false,
Exceptions::RuntimeError,
"MueLu::BlockedDirectSolver::Apply(): The map of RHS vector B is not the same as range map of the blocked operator A. Please check the map of B and A.");
146 if(bX.is_null() ==
false) {
149 TEUCHOS_TEST_FOR_EXCEPTION(bA->getFullDomainMap()->isSameAs(*(X.getMap())) ==
false,
Exceptions::RuntimeError,
"MueLu::BlockedDirectSolver::Apply(): The map of the solution vector X is not the same as domain map of the blocked operator A. Please check the map of X and A.");
153 if (bB.is_null() ==
true && bX.is_null() ==
true) {
155 s_->Apply(X, B, InitialGuessIsZero);
156 }
else if(bB.is_null() ==
false && bX.is_null() ==
false) {
158 RCP<MultiVector> mergedX = bX->Merge();
159 RCP<const MultiVector> mergedB = bB->Merge();
160 s_->Apply(*mergedX, *mergedB, InitialGuessIsZero);
161 RCP<MultiVector> xx = Teuchos::rcp(
new BlockedMultiVector(bX->getBlockedMap(),mergedX));
162 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
163 X.update(one,*xx,zero);
164 }
else if (bB.is_null() ==
true && bX.is_null() ==
false) {
166 RCP<MultiVector> mergedX = bX->Merge();
167 s_->Apply(*mergedX, B, InitialGuessIsZero);
168 RCP<MultiVector> xx = Teuchos::rcp(
new BlockedMultiVector(bX->getBlockedMap(),mergedX));
169 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
170 X.update(one,*xx,zero);
171 }
else if (bB.is_null() ==
false && bX.is_null() ==
true) {
173 RCP<const MultiVector> mergedB = bB->Merge();
174 s_->Apply(X, *mergedB, InitialGuessIsZero);
178 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
179 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
184 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
186 std::ostringstream out;
188 out <<
"{type = " << type_ <<
"}";
192 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
197 out0 <<
"Prec. type: " << type_ << std::endl;
199 if (verbLevel &
Debug)
203 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
206 return Teuchos::OrdinalTraits<size_t>::invalid();
Important warning messages (one line)
Exception indicating invalid cast attempted.
void Setup(Level ¤tLevel)
Setup routine Call the underlaying Setup routine of the nested direct solver once the input block mat...
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.
BlockedDirectSolver()
Constructor.
Timer to be used in factories. Similar to Monitor but with additional timers.
void DeclareInput(Level ¤tLevel) const
Input.
Print additional debugging information.
Namespace for MueLu classes and methods.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
RCP< const ParameterList > GetValidParameterList() const
Input.
RCP< MergedBlockedMatrixFactory > MergedAFact_
Factory to generate merged block matrix.
Class that holds all level-specific information.
direct solver for nxn blocked matrices
#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 to an FancyOStream object.
Exception throws to report errors in the internal logical of the program.
RCP< SmootherPrototype > Copy() const
RCP< DirectSolver > s_
Direct solver.
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.
std::string description() const
Return a simple one-line description of this object.