MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_StratimikosSmoother_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
47#define MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
48
49#include "MueLu_ConfigDefs.hpp"
50
51#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_THYRA)
52
53#include <Teuchos_ParameterList.hpp>
54
55#include <Xpetra_CrsMatrix.hpp>
56#include <Xpetra_CrsMatrixWrap.hpp>
57#include <Xpetra_Matrix.hpp>
58
60#include "MueLu_Level.hpp"
62#include "MueLu_Utilities.hpp"
63#include "MueLu_Monitor.hpp"
64#ifdef MUELU_RECURMG
66#endif
67
68#include <Stratimikos_DefaultLinearSolverBuilder.hpp>
69#include "Teuchos_AbstractFactoryStd.hpp"
70#if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2)
71#include <Thyra_Ifpack2PreconditionerFactory.hpp>
72#endif
73#include <Teuchos_ParameterList.hpp>
74#include <unordered_map>
75
76
77namespace MueLu {
78
79 template <class LocalOrdinal, class GlobalOrdinal, class Node>
80 StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::StratimikosSmoother(const std::string type, const Teuchos::ParameterList& paramList)
81 : type_(type)
82 {
83 std::transform(type_.begin(), type_.end(), type_.begin(), ::toupper);
84 ParameterList& pList = const_cast<ParameterList&>(paramList);
85
86 if (pList.isParameter("smoother: recurMgOnFilteredA")) {
87 recurMgOnFilteredA_ = true;
88 pList.remove("smoother: recurMgOnFilteredA");
89 }
90 bool isSupported = type_ == "STRATIMIKOS";
91 this->declareConstructionOutcome(!isSupported, "Stratimikos does not provide the smoother '" + type_ + "'.");
92 if (isSupported)
93 SetParameterList(paramList);
94 }
95
96 template <class LocalOrdinal, class GlobalOrdinal, class Node>
97 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetParameterList(const Teuchos::ParameterList& paramList) {
98 Factory::SetParameterList(paramList);
99 }
100
101 template <class LocalOrdinal, class GlobalOrdinal, class Node>
102 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& currentLevel) const {
103 this->Input(currentLevel, "A");
104 }
105
106 template <class LocalOrdinal, class GlobalOrdinal, class Node>
107 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Setup(Level& currentLevel) {
108 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
109
110 A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
111 SetupStratimikos(currentLevel);
112 SmootherPrototype::IsSetup(true);
113 this->GetOStream(Statistics1) << description() << std::endl;
114 }
115
116 template <class LocalOrdinal, class GlobalOrdinal, class Node>
117 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::SetupStratimikos(Level& currentLevel) {
118
119 RCP<const Thyra::LinearOpBase<Scalar> > thyraA;
120 if (recurMgOnFilteredA_) {
121 RCP<Matrix> filteredA;
122 ExperimentalDropVertConnections(filteredA, currentLevel);
123 thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(filteredA)->getCrsMatrix());
124 }
125 else thyraA = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(A_)->getCrsMatrix());
126
127 // Build Stratimikos solver
128 Stratimikos::DefaultLinearSolverBuilder linearSolverBuilder;
129 if (recurMgOnFilteredA_) {
130#ifdef MUELU_RECURMG
131 Stratimikos::enableMueLu<LocalOrdinal,GlobalOrdinal,Node>(linearSolverBuilder);
132#else
133 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.");
134#endif
135 }
136
137 typedef Thyra::PreconditionerFactoryBase<Scalar> Base;
138#if defined(HAVE_MUELU_THYRATPETRAADAPTERS) && defined(HAVE_MUELU_IFPACK2)
139 typedef Thyra::Ifpack2PreconditionerFactory<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > Impl;
140 linearSolverBuilder.setPreconditioningStrategyFactory(Teuchos::abstractFactoryStd<Base, Impl>(), "Ifpack2");
141#endif
142 linearSolverBuilder.setParameterList(rcpFromRef(const_cast<ParameterList&>(this->GetParameterList())));
143
144
145 // Build a new "solver factory" according to the previously specified parameter list.
146 RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > solverFactory = Thyra::createLinearSolveStrategy(linearSolverBuilder);
147 solver_ = Thyra::linearOpWithSolve(*solverFactory, thyraA);
148#ifdef dumpOutRecurMGDebug
149char mystring[120];
150sprintf(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()); fflush(stdout);
151mystring[50]='`'; mystring[65]='"'; mystring[76]='"'; mystring[77]='`';
152system(mystring);
153#endif
154 }
155
156 template <class LocalOrdinal, class GlobalOrdinal, class Node>
157 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::ExperimentalDropVertConnections(RCP<Matrix> & filteredA, Level& currentLevel) {
158
159 // strip out the veritcal connections.
160 //
161 // There is some code, which is currently turned off (via sumDropped). That
162 // makes things more complicated. I want to keep it for now, and so it is
163 // explained here. The basic idea is to try and maintain the character
164 // of the horizontal stencil by summing dropped entries appropriately.
165 // However, this does not correspond to plane relexation, and so I am
166 // not sure it is really justified. Anyway, the picture below
167 // gives a situation with a 15-pt stencil
168 //
169 // a
170 // /
171 // /
172 // b ----c----- d
173 // /
174 // /
175 // e
176 // f dropped a & l summed into f
177 // / dropped b & m summed into g
178 // / dropped c & n summed into i
179 // g ----i----- j dropped d & o summed into j
180 // / dropped e & p summed into k
181 // /
182 // k
183 // l
184 // /
185 // /
186 // m ----n----- o
187 // /
188 // /
189 // p
190 // To do this, we use umap to record locations within the middle layer associated
191 // with each line ID (e.g. g corresponds to the 13th line). Then, when dropping (in a 2nd pass) we
192 // use lineId and umap to find where dropped entries should be summed (e.g., b corresponds to the
193 // 13th line and umap[13] points at location for g).
194 // using TST = typename Teuchos::ScalarTraits<SC>;
195
196 bool sumDropped = false;
197
198 LO dofsPerNode = A_->GetFixedBlockSize();
199
200
201 RCP<ParameterList> fillCompleteParams(new ParameterList); // This code taken from Build method
202 fillCompleteParams->set("No Nonlocal Changes", true); // within MueLu_FilteredAFactory_def
203 filteredA = MatrixFactory::Build(A_->getCrsGraph());
204 filteredA->resumeFill();
205
206 ArrayView<const LocalOrdinal> inds;
207 ArrayView<const Scalar> valsA;
208 ArrayView<Scalar> vals;
209 Teuchos::ArrayRCP<LocalOrdinal> TVertLineId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel, "LineDetection_VertLineIds");
210 Teuchos::ArrayRCP<LocalOrdinal> TLayerId = Factory::Get< Teuchos::ArrayRCP<LocalOrdinal> > (currentLevel, "LineDetection_Layers");
211 LocalOrdinal* VertLineId = TVertLineId.getRawPtr();
212 LocalOrdinal* LayerId = TLayerId.getRawPtr();
213 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.");
214
215 Scalar ZERO = Teuchos::ScalarTraits<Scalar>::zero();
216 for (size_t i = 0; i < A_->getRowMap()->getLocalNumElements(); i++) {
217 A_->getLocalRowView(i, inds, valsA);
218 size_t nnz = inds.size();
219 ArrayView<const Scalar> vals1;
220 filteredA->getLocalRowView(i, inds, vals1);
221 vals = ArrayView<Scalar>(const_cast<Scalar*>(vals1.getRawPtr()), nnz);
222 memcpy(vals.getRawPtr(), valsA.getRawPtr(), nnz*sizeof(Scalar));
223 size_t inode, jdof, jnode, jdof_offset;
224 inode = i/dofsPerNode;
225
226 std::unordered_map<LocalOrdinal, LocalOrdinal> umap; // umap[k] indicates where the dropped entry
227 // corresponding to kth line should be added
228 // within the row. See comments above.
229 if (sumDropped) {
230 for (size_t j = 0; j < nnz; j++) {
231 jdof = inds[j];
232 jnode = jdof/dofsPerNode;
233 jdof_offset = jdof - jnode*dofsPerNode;
234 if ( LayerId[jnode] == LayerId[inode] ) umap[dofsPerNode*VertLineId[jnode]+jdof_offset] = j;
235 }
236 }
237
238 // drop non-middle layer entries. When sumDropped=true, sum dropped entries to corresponding mid-layer entry
239 for (size_t j = 0; j < nnz; j++) {
240 jdof = inds[j];
241 jnode = jdof/dofsPerNode;
242 jdof_offset = jdof - jnode*dofsPerNode;
243 if ( LayerId[jnode] != LayerId[inode] ) {
244 if (sumDropped) {
245 if (umap.find(dofsPerNode*VertLineId[jnode + jdof_offset]) != umap.end())
246 vals[umap[dofsPerNode*VertLineId[jnode + jdof_offset]]] += vals[j];
247 }
248 vals[j] = ZERO;
249 }
250 }
251 }
252 filteredA->fillComplete(fillCompleteParams);
253
254 }
255
256 template <class LocalOrdinal, class GlobalOrdinal, class Node>
257 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
258 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::StratimikosSmoother::Apply(): Setup() has not been called");
259
260 // Apply
261 if (InitialGuessIsZero) {
262 X.putScalar(0.0);
263 RCP< Thyra::MultiVectorBase<Scalar> > thyraX = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(X)));
264 RCP<const Thyra::MultiVectorBase<Scalar> > thyraB = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(rcpFromRef(B));
265 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraB, thyraX.ptr());
266 RCP<MultiVector> thyXpX = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraX, X.getMap()->getComm());
267 X = *thyXpX;
268
269 } else {
270 typedef Teuchos::ScalarTraits<Scalar> TST;
271 RCP<MultiVector> Residual = Utilities::Residual(*A_, X, B);
272
273 RCP<MultiVector> Cor = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(X.getMap(),X.getNumVectors(), true);
274 RCP< Thyra::MultiVectorBase<Scalar> > thyraCor = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Cor));
275 RCP<const Thyra::MultiVectorBase<Scalar> > thyraRes = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(Residual);
276 Thyra::SolveStatus<Scalar> status = Thyra::solve<Scalar>(*solver_, Thyra::NOTRANS, *thyraRes, thyraCor.ptr());
277 RCP<MultiVector> thyXpCor = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraCor, X.getMap()->getComm());
278 X.update(TST::one(), *thyXpCor, TST::one());
279 }
280
281 }
282
283
284 template <class LocalOrdinal, class GlobalOrdinal, class Node>
285 RCP<MueLu::SmootherPrototype<double, LocalOrdinal, GlobalOrdinal, Node> > StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::Copy() const {
286 RCP<StratimikosSmoother> smoother = rcp(new StratimikosSmoother(*this) );
287 smoother->SetParameterList(this->GetParameterList());
288 return smoother;
289 }
290
291 template <class LocalOrdinal, class GlobalOrdinal, class Node>
292 std::string StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::description() const {
293 std::ostringstream out;
294 if (SmootherPrototype::IsSetup()) {
295 out << solver_->description();
296 } else {
297 out << "STRATIMIKOS {type = " << type_ << "}";
298 }
299 return out.str();
300 }
301
302 template <class LocalOrdinal, class GlobalOrdinal, class Node>
303 void StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
305
306 if (verbLevel & Parameters1) {
307 out0 << "Parameter list: " << std::endl;
308 Teuchos::OSTab tab2(out);
309 out << this->GetParameterList();
310 }
311
312 if (verbLevel & External)
313 if (solver_ != Teuchos::null) {
314 Teuchos::OSTab tab2(out);
315 out << *solver_ << std::endl;
316 }
317
318 if (verbLevel & Debug) {
319 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
320 << "-" << std::endl
321 << "RCP<solver_>: " << solver_ << std::endl;
322 }
323 }
324
325 template <class LocalOrdinal, class GlobalOrdinal, class Node>
326 size_t StratimikosSmoother<double, LocalOrdinal, GlobalOrdinal, Node>::getNodeSmootherComplexity() const {
327 return Teuchos::OrdinalTraits<size_t>::invalid();
328 }
329
330
331} // namespace MueLu
332
333#endif // HAVE_MUELU_STRATIMIKOS
334#endif // MUELU_STRATIMIKOSSMOOTHER_DEF_HPP
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
Namespace for MueLu classes and methods.