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();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
direct solver for nxn blocked matrices
void Setup(Level ¤tLevel)
Setup routine Call the underlaying Setup routine of the nested direct solver once the input block mat...
RCP< const ParameterList > GetValidParameterList() const
Input.
void DeclareInput(Level ¤tLevel) const
Input.
RCP< DirectSolver > s_
Direct solver.
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.
std::string description() const
Return a simple one-line description of this object.
BlockedDirectSolver()
Constructor.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
RCP< MergedBlockedMatrixFactory > MergedAFact_
Factory to generate merged block matrix.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< SmootherPrototype > Copy() const
virtual std::string description() const
Return a simple one-line description of this object.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
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()
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.