46 #ifndef MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
47 #define MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
51 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
56 #include <Xpetra_CrsMatrixWrap.hpp>
61 #include "MueLu_Utilities.hpp"
67 #include <Stratimikos_DefaultLinearSolverBuilder.hpp>
68 #include "Teuchos_AbstractFactoryStd.hpp"
70 #include <unordered_map>
74 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::StratimikosSmoother(
const std::string type,
const Teuchos::ParameterList& paramList)
77 std::transform(type_.begin(), type_.end(), type_.begin(), ::toupper);
78 ParameterList& pList =
const_cast<ParameterList&
>(paramList);
80 if (pList.isParameter(
"smoother: recurMgOnFilteredA")) {
81 recurMgOnFilteredA_ =
true;
82 pList.remove(
"smoother: recurMgOnFilteredA");
84 bool isSupported = type_ ==
"STRATIMIKOS";
85 this->declareConstructionOutcome(!isSupported,
"Stratimikos does not provide the smoother '" + type_ +
"'.");
87 SetParameterList(paramList);
90 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
91 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetParameterList(
const Teuchos::ParameterList& paramList) {
92 Factory::SetParameterList(paramList);
95 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel)
const {
97 this->Input(currentLevel,
"A");
100 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
101 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
102 FactoryMonitor m(*
this,
"Setup Smoother", currentLevel);
104 A_ = Factory::Get<RCP<Matrix> >(currentLevel,
"A");
105 SetupStratimikos(currentLevel);
106 SmootherPrototype::IsSetup(
true);
107 this->GetOStream(
Statistics1) << description() << std::endl;
110 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetupStratimikos(Level& currentLevel) {
112 RCP<const Thyra::LinearOpBase<Scalar> > thyraA;
113 if (recurMgOnFilteredA_) {
114 RCP<Matrix> filteredA;
115 ExperimentalDropVertConnections(filteredA, currentLevel);
116 thyraA = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix());
118 thyraA = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A_)->getCrsMatrix());
121 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
122 if (recurMgOnFilteredA_) {
124 Stratimikos::enableMueLu<LocalOrdinal, GlobalOrdinal, Node>(linearSolverBuilder);
126 TEUCHOS_TEST_FOR_EXCEPTION(
true, Exceptions::RuntimeError,
"MueLu::StratimikosSmoother:: must compile with MUELU_RECURMG defined. Unfortunately, cmake does not always produce a proper link.txt file (which sometimes requires libmuelu.a before and after libmuelu-interface.a). After configuring, run script muelu/utils/misc/patchLinkForRecurMG to change link.txt files manually. If you want to create test example, add -DMUELU_RECURMG=ON to cmake arguments.");
130 linearSolverBuilder.setParameterList(rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
133 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
134 solver_ = Thyra::linearOpWithSolve(*solverFactory, thyraA);
135 #ifdef dumpOutRecurMGDebug
137 sprintf(mystring,
"for i in A_[0123456789].m P_[0123456789].m; do T=Xecho $i | sed Xs/.m$/%d.m/XX; mv $i $T; done", (
int)currentLevel.GetLevelID());
147 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
148 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::ExperimentalDropVertConnections(RCP<Matrix>& filteredA, Level& currentLevel) {
186 bool sumDropped =
false;
188 LO dofsPerNode = A_->GetFixedBlockSize();
190 RCP<ParameterList> fillCompleteParams(
new ParameterList);
191 fillCompleteParams->set(
"No Nonlocal Changes",
true);
192 filteredA = MatrixFactory::Build(A_->getCrsGraph());
193 filteredA->resumeFill();
195 ArrayView<const LocalOrdinal> inds;
196 ArrayView<const Scalar> valsA;
197 ArrayView<Scalar> vals;
202 TEUCHOS_TEST_FOR_EXCEPTION((LayerId == NULL) || (VertLineId == NULL), Exceptions::RuntimeError,
"MueLu::StratimikosSmoother:: no line information found on this level. Cannot use recurMgOnFilteredA on this level.");
205 for (
size_t i = 0; i < A_->getRowMap()->getLocalNumElements(); i++) {
206 A_->getLocalRowView(i, inds, valsA);
207 size_t nnz = inds.size();
208 ArrayView<const Scalar> vals1;
209 filteredA->getLocalRowView(i, inds, vals1);
210 vals = ArrayView<Scalar>(
const_cast<Scalar*
>(vals1.getRawPtr()), nnz);
211 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz *
sizeof(
Scalar));
212 size_t inode, jdof, jnode, jdof_offset;
213 inode = i / dofsPerNode;
215 std::unordered_map<LocalOrdinal, LocalOrdinal> umap;
219 for (
size_t j = 0; j < nnz; j++) {
221 jnode = jdof / dofsPerNode;
222 jdof_offset = jdof - jnode * dofsPerNode;
223 if (LayerId[jnode] == LayerId[inode]) umap[dofsPerNode * VertLineId[jnode] + jdof_offset] = j;
228 for (
size_t j = 0; j < nnz; j++) {
230 jnode = jdof / dofsPerNode;
231 jdof_offset = jdof - jnode * dofsPerNode;
232 if (LayerId[jnode] != LayerId[inode]) {
234 if (umap.find(dofsPerNode * VertLineId[jnode + jdof_offset]) != umap.end())
235 vals[umap[dofsPerNode * VertLineId[jnode + jdof_offset]]] += vals[j];
241 filteredA->fillComplete(fillCompleteParams);
244 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
245 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X,
const MultiVector& B,
bool InitialGuessIsZero)
const {
246 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() ==
false, Exceptions::RuntimeError,
"MueLu::StratimikosSmoother::Apply(): Setup() has not been called");
249 if (InitialGuessIsZero) {
251 RCP<Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpFromRef(X)));
252 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpFromRef(B));
253 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraB, thyraX.ptr());
254 RCP<MultiVector> thyXpX = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyraX, X.getMap()->getComm());
259 RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
262 RCP<Thyra::MultiVectorBase<Scalar> > thyraCor = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(Cor));
263 RCP<const Thyra::MultiVectorBase<Scalar> > thyraRes = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(Residual);
264 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraRes, thyraCor.ptr());
265 RCP<MultiVector> thyXpCor = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyraCor, X.getMap()->getComm());
266 X.update(TST::one(), *thyXpCor, TST::one());
270 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
271 RCP<MueLu::SmootherPrototype<double, LocalOrdinal, GlobalOrdinal, Node> > StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Copy()
const {
272 RCP<StratimikosSmoother> smoother =
rcp(
new StratimikosSmoother(*
this));
273 smoother->SetParameterList(this->GetParameterList());
277 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
278 std::string StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::description()
const {
279 std::ostringstream out;
280 if (SmootherPrototype::IsSetup()) {
281 out << solver_->description();
283 out <<
"STRATIMIKOS {type = " << type_ <<
"}";
288 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
293 out0 <<
"Parameter list: " << std::endl;
295 out << this->GetParameterList();
299 if (solver_ != Teuchos::null) {
301 out << *solver_ << std::endl;
304 if (verbLevel &
Debug) {
305 out0 <<
"IsSetup: " <<
Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
307 <<
"RCP<solver_>: " << solver_ << std::endl;
311 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 size_t StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::getNodeSmootherComplexity()
const {
318 #endif // HAVE_MUELU_STRATIMIKOS
319 #endif // MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Print external lib objects.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print additional debugging information.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters (more parameters, more verbose)
std::string toString(const T &t)