MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_StratimikosSmoother_def.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// MueLu: A package for multigrid based preconditioning
4//
5// Copyright 2012 NTESS and the MueLu contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
11#define MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
12
13#include "MueLu_ConfigDefs.hpp"
14
15#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
16
17#include <Teuchos_ParameterList.hpp>
18
19#include <Xpetra_CrsMatrix.hpp>
20#include <Xpetra_CrsMatrixWrap.hpp>
21#include <Xpetra_Matrix.hpp>
22
23#include "MueLu_StratimikosSmoother_decl.hpp"
24#include "MueLu_Level.hpp"
25#include "MueLu_Utilities.hpp"
26#include "MueLu_Monitor.hpp"
27#ifdef MUELU_RECURMG
28#include "Stratimikos_MueLuHelpers.hpp"
29#endif
30
31#include <Stratimikos_DefaultLinearSolverBuilder.hpp>
32#include "Teuchos_AbstractFactoryStd.hpp"
33#include <Teuchos_ParameterList.hpp>
34#include <unordered_map>
35
36namespace MueLu {
37
38template <class LocalOrdinal, class GlobalOrdinal, class Node>
39StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::StratimikosSmoother(const std::string type, const Teuchos::ParameterList& paramList)
40 : type_(type) {
41 std::transform(type_.begin(), type_.end(), type_.begin(), ::toupper);
42 ParameterList& pList = const_cast<ParameterList&>(paramList);
43
44 if (pList.isParameter("smoother: recurMgOnFilteredA")) {
45 recurMgOnFilteredA_ = true;
46 pList.remove("smoother: recurMgOnFilteredA");
47 }
48 bool isSupported = type_ == "STRATIMIKOS";
49 this->declareConstructionOutcome(!isSupported, "Stratimikos does not provide the smoother '" + type_ + "'.");
50 if (isSupported)
51 SetParameterList(paramList);
52}
53
54template <class LocalOrdinal, class GlobalOrdinal, class Node>
55void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetParameterList(const Teuchos::ParameterList& paramList) {
56 Factory::SetParameterList(paramList);
57}
58
59template <class LocalOrdinal, class GlobalOrdinal, class Node>
60void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
61 this->Input(currentLevel, "A");
62}
63
64template <class LocalOrdinal, class GlobalOrdinal, class Node>
65void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
66 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
67
68 A_ = Factory::Get<RCP<Matrix> >(currentLevel, "A");
69 SetupStratimikos(currentLevel);
70 SmootherPrototype::IsSetup(true);
71 this->GetOStream(Statistics1) << description() << std::endl;
72}
73
74template <class LocalOrdinal, class GlobalOrdinal, class Node>
75void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetupStratimikos(Level& currentLevel) {
76 RCP<const Thyra::LinearOpBase<Scalar> > thyraA;
77 if (recurMgOnFilteredA_) {
78 RCP<Matrix> filteredA;
79 ExperimentalDropVertConnections(filteredA, currentLevel);
80 thyraA = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(toCrsMatrix(filteredA));
81 } else
82 thyraA = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(toCrsMatrix(A_));
83
84 // Build Stratimikos solver
85 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
86 if (recurMgOnFilteredA_) {
87#ifdef MUELU_RECURMG
88 Stratimikos::enableMueLu<LocalOrdinal, GlobalOrdinal, Node>(linearSolverBuilder);
89#else
90 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.");
91#endif
92 }
93
94 linearSolverBuilder.setParameterList(rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
95
96 // Build a new "solver factory" according to the previously specified parameter list.
97 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
98 solver_ = Thyra::linearOpWithSolve(*solverFactory, thyraA);
99#ifdef dumpOutRecurMGDebug
100 char mystring[120];
101 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());
102 fflush(stdout);
103 mystring[50] = '`';
104 mystring[65] = '"';
105 mystring[76] = '"';
106 mystring[77] = '`';
107 system(mystring);
108#endif
109}
110
111template <class LocalOrdinal, class GlobalOrdinal, class Node>
112void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::ExperimentalDropVertConnections(RCP<Matrix>& filteredA, Level& currentLevel) {
113 // strip out the veritcal connections.
114 //
115 // There is some code, which is currently turned off (via sumDropped). That
116 // makes things more complicated. I want to keep it for now, and so it is
117 // explained here. The basic idea is to try and maintain the character
118 // of the horizontal stencil by summing dropped entries appropriately.
119 // However, this does not correspond to plane relexation, and so I am
120 // not sure it is really justified. Anyway, the picture below
121 // gives a situation with a 15-pt stencil
122 //
123 // a
124 // /
125 // /
126 // b ----c----- d
127 // /
128 // /
129 // e
130 // f dropped a & l summed into f
131 // / dropped b & m summed into g
132 // / dropped c & n summed into i
133 // g ----i----- j dropped d & o summed into j
134 // / dropped e & p summed into k
135 // /
136 // k
137 // l
138 // /
139 // /
140 // m ----n----- o
141 // /
142 // /
143 // p
144 // To do this, we use umap to record locations within the middle layer associated
145 // with each line ID (e.g. g corresponds to the 13th line). Then, when dropping (in a 2nd pass) we
146 // use lineId and umap to find where dropped entries should be summed (e.g., b corresponds to the
147 // 13th line and umap[13] points at location for g).
148 // using TST = typename Teuchos::ScalarTraits<SC>;
149
150 bool sumDropped = false;
151
152 LO dofsPerNode = A_->GetFixedBlockSize();
153
154 RCP<ParameterList> fillCompleteParams(new ParameterList); // This code taken from Build method
155 fillCompleteParams->set("No Nonlocal Changes", true); // within MueLu_FilteredAFactory_def
156 filteredA = MatrixFactory::Build(A_->getCrsGraph());
157 filteredA->resumeFill();
158
159 ArrayView<const LocalOrdinal> inds;
160 ArrayView<const Scalar> valsA;
161 ArrayView<Scalar> vals;
162 Teuchos::ArrayRCP<LocalOrdinal> TVertLineId = Factory::Get<Teuchos::ArrayRCP<LocalOrdinal> >(currentLevel, "LineDetection_VertLineIds");
163 Teuchos::ArrayRCP<LocalOrdinal> TLayerId = Factory::Get<Teuchos::ArrayRCP<LocalOrdinal> >(currentLevel, "LineDetection_Layers");
164 LocalOrdinal* VertLineId = TVertLineId.getRawPtr();
165 LocalOrdinal* LayerId = TLayerId.getRawPtr();
166 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.");
167
168 Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
169 for (size_t i = 0; i < A_->getRowMap()->getLocalNumElements(); i++) {
170 A_->getLocalRowView(i, inds, valsA);
171 size_t nnz = inds.size();
172 ArrayView<const Scalar> vals1;
173 filteredA->getLocalRowView(i, inds, vals1);
174 vals = ArrayView<Scalar>(const_cast<Scalar*>(vals1.getRawPtr()), nnz);
175 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz * sizeof(Scalar));
176 size_t inode, jdof, jnode, jdof_offset;
177 inode = i / dofsPerNode;
178
179 std::unordered_map<LocalOrdinal, LocalOrdinal> umap; // umap[k] indicates where the dropped entry
180 // corresponding to kth line should be added
181 // within the row. See comments above.
182 if (sumDropped) {
183 for (size_t j = 0; j < nnz; j++) {
184 jdof = inds[j];
185 jnode = jdof / dofsPerNode;
186 jdof_offset = jdof - jnode * dofsPerNode;
187 if (LayerId[jnode] == LayerId[inode]) umap[dofsPerNode * VertLineId[jnode] + jdof_offset] = j;
188 }
189 }
190
191 // drop non-middle layer entries. When sumDropped=true, sum dropped entries to corresponding mid-layer entry
192 for (size_t j = 0; j < nnz; j++) {
193 jdof = inds[j];
194 jnode = jdof / dofsPerNode;
195 jdof_offset = jdof - jnode * dofsPerNode;
196 if (LayerId[jnode] != LayerId[inode]) {
197 if (sumDropped) {
198 if (umap.find(dofsPerNode * VertLineId[jnode + jdof_offset]) != umap.end())
199 vals[umap[dofsPerNode * VertLineId[jnode + jdof_offset]]] += vals[j];
200 }
201 vals[j] = ZERO;
202 }
203 }
204 }
205 filteredA->fillComplete(fillCompleteParams);
206}
207
208template <class LocalOrdinal, class GlobalOrdinal, class Node>
209void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
210 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::StratimikosSmoother::Apply(): Setup() has not been called");
211
212 // Apply
213 if (InitialGuessIsZero) {
214 X.putScalar(0.0);
215 RCP<Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpFromRef(X)));
216 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(rcpFromRef(B));
217 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraB, thyraX.ptr());
218 RCP<MultiVector> thyXpX = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyraX, X.getMap()->getComm());
219 X = *thyXpX;
220
221 } else {
222 typedef Teuchos::ScalarTraits<Scalar> TST;
223 RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
224
225 RCP<MultiVector> Cor = Xpetra::MultiVectorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(X.getMap(), X.getNumVectors(), true);
226 RCP<Thyra::MultiVectorBase<Scalar> > thyraCor = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(Cor));
227 RCP<const Thyra::MultiVectorBase<Scalar> > thyraRes = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyraMultiVector(Residual);
228 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraRes, thyraCor.ptr());
229 RCP<MultiVector> thyXpCor = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toXpetra(thyraCor, X.getMap()->getComm());
230 X.update(TST::one(), *thyXpCor, TST::one());
231 }
232}
233
234template <class LocalOrdinal, class GlobalOrdinal, class Node>
235RCP<MueLu::SmootherPrototype<double, LocalOrdinal, GlobalOrdinal, Node> > StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
236 RCP<StratimikosSmoother> smoother = rcp(new StratimikosSmoother(*this));
237 smoother->SetParameterList(this->GetParameterList());
238 return smoother;
239}
240
241template <class LocalOrdinal, class GlobalOrdinal, class Node>
242std::string StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::description() const {
243 std::ostringstream out;
244 if (SmootherPrototype::IsSetup()) {
245 out << solver_->description();
246 } else {
247 out << "STRATIMIKOS {type = " << type_ << "}";
248 }
249 return out.str();
250}
251
252template <class LocalOrdinal, class GlobalOrdinal, class Node>
253void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
254 MUELU_DESCRIBE;
255
256 if (verbLevel & Parameters1) {
257 out0 << "Parameter list: " << std::endl;
258 Teuchos::OSTab tab2(out);
259 out << this->GetParameterList();
260 }
261
262 if (verbLevel & External)
263 if (solver_ != Teuchos::null) {
264 Teuchos::OSTab tab2(out);
265 out << *solver_ << std::endl;
266 }
267
268 if (verbLevel & Debug) {
269 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
270 << "-" << std::endl
271 << "RCP<solver_>: " << solver_ << std::endl;
272 }
273}
274
275template <class LocalOrdinal, class GlobalOrdinal, class Node>
276size_t StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::getNodeSmootherComplexity() const {
277 return Teuchos::OrdinalTraits<size_t>::invalid();
278}
279
280} // namespace MueLu
281
282#endif // HAVE_MUELU_STRATIMIKOS
283#endif // MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
MUELU_DESCRIBE
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Definition MueLu_BaseClass.hpp:44
MueLu_ConfigDefs.hpp
MueLu_Level.hpp
MueLu_Monitor.hpp
MueLu_StratimikosSmoother_decl.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition MueLu_UseDefaultTypes.hpp:13
Scalar
MueLu::DefaultScalar Scalar
Definition MueLu_UseDefaultTypes.hpp:12
Stratimikos_MueLuHelpers.hpp
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42