MueLu Version of the Day
Loading...
Searching...
No Matches
Thyra_MueLuPreconditionerFactory_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 THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
11#define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
12
13#include "Thyra_MueLuPreconditionerFactory_decl.hpp"
14#include "MueLu_TpetraOperatorAsRowMatrix.hpp"
15#include "MueLu_MasterList.hpp"
16
17#if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
18
19// This is not as general as possible, but should be good enough for most builds.
20#if ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
21 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
22 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)))
23#define MUELU_CAN_USE_MIXED_PRECISION
24#endif
25
26namespace Thyra {
27
28using Teuchos::ParameterList;
29using Teuchos::RCP;
30using Teuchos::rcp;
31using Teuchos::rcp_const_cast;
32using Teuchos::rcp_dynamic_cast;
33
34template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
35bool Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
36 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
37 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
38 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
39 // typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
40 // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
41 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
42 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMultVec;
43 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpMagMultVec;
44 typedef Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpVec;
45
46 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
47 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
48 // typedef Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyXpOp;
49 // typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
50
51 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpCrsMat;
52 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
53 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
54 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
55 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
56 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
57#if defined(MUELU_CAN_USE_MIXED_PRECISION)
58 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
59 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
60#endif
61
62 if (paramList.isParameter(parameterName)) {
63 if (paramList.isType<RCP<XpMat> >(parameterName))
64 return true;
65 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
66 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
67 paramList.remove(parameterName);
68 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
69 paramList.set<RCP<XpMat> >(parameterName, M);
70 return true;
71 } else if (paramList.isType<RCP<XpMultVec> >(parameterName))
72 return true;
73 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
74 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
75 paramList.remove(parameterName);
76 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
77 paramList.set<RCP<XpMultVec> >(parameterName, X);
78 return true;
79 } else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
80 return true;
81 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
82 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
83 paramList.remove(parameterName);
84 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
85 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
86 return true;
87 } else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
88 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
89 paramList.remove(parameterName);
90 RCP<XpMat> xM = Xpetra::toXpetra(tM);
91 paramList.set<RCP<XpMat> >(parameterName, xM);
92 return true;
93 } else if (paramList.isType<RCP<tMV> >(parameterName)) {
94 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
95 paramList.remove(parameterName);
96 RCP<XpMultVec> X = Xpetra::toXpetra(tpetra_X);
97 paramList.set<RCP<XpMultVec> >(parameterName, X);
98 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
99 return true;
100 } else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
101 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
102 paramList.remove(parameterName);
103 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
104 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
105 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
106 return true;
107 }
108#if defined(MUELU_CAN_USE_MIXED_PRECISION)
109 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
110 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
111 paramList.remove(parameterName);
112 RCP<tMagMV> tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(), tpetra_hX->getNumVectors()));
113 Tpetra::deep_copy(*tpetra_X, *tpetra_hX);
114 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
115 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
116 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
117 return true;
118 }
119#endif
120 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
121 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).is_null())) {
122 RCP<const ThyDiagLinOpBase> thyM;
123 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
124 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
125 else
126 thyM = rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName), true);
127 paramList.remove(parameterName);
128 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
129
130 RCP<const XpVec> xpDiag;
131 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
132 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
133 if (!tDiag.is_null())
134 xpDiag = Xpetra::toXpetra(tDiag);
135 }
136 TEUCHOS_ASSERT(!xpDiag.is_null());
137 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
138 paramList.set<RCP<XpMat> >(parameterName, M);
139 return true;
140 } else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
141 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
142 paramList.remove(parameterName);
143 try {
144 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
145 paramList.set<RCP<XpMat> >(parameterName, M);
146 } catch (std::exception& e) {
147 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
148 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(M, true);
149 RCP<tOp> tO = tpOp->getOperator();
150 RCP<tV> diag;
151 if (tO->hasDiagonal()) {
152 diag = rcp(new tV(tO->getRangeMap()));
153 tO->getLocalDiagCopy(*diag);
154 }
155 auto fTpRow = rcp(new MueLu::TpetraOperatorAsRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tO, diag));
156 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpFOp = rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fTpRow));
157 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
158 paramList.set<RCP<XpOp> >(parameterName, op);
159 }
160 return true;
161 } else {
162 TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type.");
163 return false;
164 }
165 } else
166 return false;
167}
168
169// Constructors/initializers/accessors
170
171template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
172MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
173 : paramList_(rcp(new ParameterList())) {}
174
175template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
176MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() = default;
177
178// Overridden from PreconditionerFactoryBase
179
180template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
181bool MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
182 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
183
184 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp)) return true;
185
186 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isBlockedOperator(fwdOp)) return true;
187
188 return false;
189}
190
191template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
192RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
193 return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
194}
195
196template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
197void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
198 initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
199 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLu::initializePrec")));
200
201 using Teuchos::rcp_dynamic_cast;
202
203 // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
204 typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> XpMap;
205 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
206 typedef MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> MueLuXpOp;
207 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
208 // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
209 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpBlockedCrsMat;
210 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
211 // typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
212 // typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
213 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
214 typedef Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyXpOp;
215 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMV;
216 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
217 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpmMV;
218#if defined(MUELU_CAN_USE_MIXED_PRECISION)
219 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
220 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
221 typedef Xpetra::Operator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfOp;
222 typedef MueLu::XpetraOperator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> MueLuHalfXpOp;
223 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
224 typedef Xpetra::MultiVector<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMV;
225 typedef Xpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> XphmMV;
226 typedef Xpetra::Matrix<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMat;
227#endif
228
229 // Check precondition
230 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
231 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
232 TEUCHOS_ASSERT(prec);
233
234 // Create a copy, as we may remove some things from the list
235 ParameterList paramList = *paramList_;
236
237 // Retrieve wrapped concrete Xpetra matrix from FwdOp
238 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
239 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
240
241 // Check whether it is Tpetra
242 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
243 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
244 TEUCHOS_TEST_FOR_EXCEPT((bIsTpetra == false) && (bIsBlocked == false));
245 TEUCHOS_TEST_FOR_EXCEPT((bIsTpetra == true) && bIsBlocked == true);
246
247 RCP<XpMat> A = Teuchos::null;
248 if (bIsBlocked) {
249 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
250 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
251 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
252
253 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) == false);
254
255 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
256 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
257
258 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
259 // MueLu needs a non-const object as input
260 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
261 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
262
263 RCP<const XpMap> rowmap00 = A00->getRowMap();
264 RCP<const Teuchos::Comm<int> > comm = rowmap00->getComm();
265
266 // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
267 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
268 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
269
270 // save blocked matrix
271 A = bMat;
272 } else {
273 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
274 // MueLu needs a non-const object as input
275 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
276 }
277 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
278
279 // Retrieve concrete preconditioner object
280 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
281 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
282
283 // extract preconditioner operator
284 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
285 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
286
287 // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
288 // rebuild preconditioner if startingOver == true
289 // reuse preconditioner if startingOver == false
290 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
291 bool useHalfPrecision = false;
292 if (paramList.isParameter("half precision"))
293 useHalfPrecision = paramList.get<bool>("half precision");
294 else if (paramList.isSublist("Hierarchy") && paramList.sublist("Hierarchy").isParameter("half precision"))
295 useHalfPrecision = paramList.sublist("Hierarchy").get<bool>("half precision");
296
297 RCP<XpOp> xpPrecOp;
298 if (startingOver == true) {
299 // Convert to Xpetra
300 std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "Material", "BlockNumber"};
301 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
302 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
303
304 if (paramList.isSublist("user data")) {
305 auto& userData = paramList.sublist("user data");
306 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
307 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(userData, *it);
308 }
309
310 for (int lvlNo = 0; lvlNo < 10; ++lvlNo) {
311 if (paramList.isSublist("level " + std::to_string(lvlNo) + " user data")) {
312 ParameterList& lvlList = paramList.sublist("level " + std::to_string(lvlNo) + " user data");
313 std::list<std::string> convertKeys;
314 for (auto it = lvlList.begin(); it != lvlList.end(); ++it)
315 convertKeys.push_back(lvlList.name(it));
316 for (auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
317 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
318 }
319 }
320
321 if (useHalfPrecision) {
322#if defined(MUELU_CAN_USE_MIXED_PRECISION)
323
324 // CAG: There is nothing special about the combination double-float,
325 // except that I feel somewhat confident that Trilinos builds
326 // with both scalar types.
327
328 // convert to half precision
329 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
330 const std::string userName = "user data";
331 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
332 if (userParamList.isType<RCP<XpmMV> >("Coordinates")) {
333 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >("Coordinates");
334 userParamList.remove("Coordinates");
335 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
336 userParamList.set("Coordinates", halfCoords);
337 }
338 if (userParamList.isType<RCP<XpMV> >("Nullspace")) {
339 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >("Nullspace");
340 userParamList.remove("Nullspace");
341 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
342 userParamList.set("Nullspace", halfNullspace);
343 }
344 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
345 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
346 paramList.remove("Coordinates");
347 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
348 userParamList.set("Coordinates", halfCoords);
349 }
350 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
351 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
352 paramList.remove("Nullspace");
353 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
354 userParamList.set("Nullspace", halfNullspace);
355 }
356
357 // build a new half-precision MueLu preconditioner
358
359 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(halfA, paramList);
360 RCP<MueLuHalfXpOp> xpOp = rcp(new MueLuHalfXpOp(H));
361 xpPrecOp = rcp(new XpHalfPrecOp(xpOp));
362#else
363 TEUCHOS_TEST_FOR_EXCEPT(true);
364#endif
365 } else {
366 const std::string userName = "user data";
367 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
368 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
369 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
370 paramList.remove("Coordinates");
371 userParamList.set("Coordinates", coords);
372 }
373 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
374 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
375 paramList.remove("Nullspace");
376 userParamList.set("Nullspace", nullspace);
377 }
378
379 // build a new MueLu RefMaxwell preconditioner
380 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(A, paramList);
381 xpPrecOp = rcp(new MueLuXpOp(H));
382 }
383 } else {
384 // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix
385 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
386 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(), true);
387#if defined(MUELU_CAN_USE_MIXED_PRECISION)
388 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
389 if (!xpHalfPrecOp.is_null()) {
390 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy();
391 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
392
393 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
394 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
395 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
396 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
397 RCP<MueLu::Level> level0 = H->GetLevel(0);
398 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >("A");
399 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0, true);
400
401 if (!A0.is_null()) {
402 // If a user provided a "number of equations" argument in a parameter list
403 // during the initial setup, we must honor that settings and reuse it for
404 // all consequent setups.
405 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
406 }
407
408 // set new matrix
409 level0->Set("A", halfA);
410
411 H->SetupRe();
412 } else
413#endif
414 {
415 // get old MueLu hierarchy
416 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(), true);
417 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
418 ;
419
420 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
421 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
422 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
423 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
424 RCP<MueLu::Level> level0 = H->GetLevel(0);
425 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
426 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
427
428 if (!A0.is_null()) {
429 // If a user provided a "number of equations" argument in a parameter list
430 // during the initial setup, we must honor that settings and reuse it for
431 // all consequent setups.
432 A->SetFixedBlockSize(A0->GetFixedBlockSize());
433 }
434
435 // set new matrix
436 level0->Set("A", A);
437
438 H->SetupRe();
439 }
440 }
441
442 // wrap preconditioner in thyraPrecOp
443 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
444 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
445
446 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
447 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
448
449 defaultPrec->initializeUnspecified(thyraPrecOp);
450}
451
452template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
453void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
454 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
455 TEUCHOS_ASSERT(prec);
456
457 // Retrieve concrete preconditioner object
458 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
459 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
460
461 if (fwdOp) {
462 // TODO: Implement properly instead of returning default value
463 *fwdOp = Teuchos::null;
464 }
465
466 if (supportSolveUse) {
467 // TODO: Implement properly instead of returning default value
468 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
469 }
470
471 defaultPrec->uninitialize();
472}
473
474// Overridden from ParameterListAcceptor
475template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
476void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
477 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
478 paramList_ = paramList;
479}
480
481template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
482RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
483 return paramList_;
484}
485
486template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
487RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
488 RCP<ParameterList> savedParamList = paramList_;
489 paramList_ = Teuchos::null;
490 return savedParamList;
491}
492
493template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
494RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
495 return paramList_;
496}
497
498template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
499RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters() const {
500 static RCP<const ParameterList> validPL;
501
502 if (Teuchos::is_null(validPL)) {
503 validPL = MueLu::MasterList::List();
504 }
505
506 return validPL;
507}
508
509// Public functions overridden from Teuchos::Describable
510template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
511std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
512 return "Thyra::MueLuPreconditionerFactory";
513}
514} // namespace Thyra
515
516#endif // HAVE_MUELU_STRATIMIKOS
517
518#endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
MueLu_MasterList.hpp
MueLu_TpetraOperatorAsRowMatrix.hpp
Thyra_MueLuPreconditionerFactory_decl.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition MueLu_Exceptions.hpp:36
MueLu::MasterList::List
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a "master" list of all valid parameters and their default values.
Definition MueLu_MasterList.cpp:18
MueLu::XpetraOperator
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.
Definition MueLu_XpetraOperator_decl.hpp:28
Thyra::XpetraLinearOp
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
Definition Thyra_XpetraLinearOp_decl.hpp:28
MueLu::CreateXpetraPreconditioner
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix,...
Definition MueLu_CreateXpetraPreconditioner.hpp:45