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<ThyLinOpBase> >(parameterName)) {
141 RCP<ThyLinOpBase> thyM = paramList.get<RCP<ThyLinOpBase> >(parameterName);
142 paramList.remove(parameterName);
143 try {
144 RCP<XpMat> M = XpThyUtils::toXpetra(thyM);
145 paramList.set<RCP<XpMat> >(parameterName, M);
146 } catch (std::exception& e) {
147 RCP<XpOp> M = XpThyUtils::toXpetraOperator(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 if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
162 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
163 paramList.remove(parameterName);
164 try {
165 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
166 paramList.set<RCP<XpMat> >(parameterName, M);
167 } catch (std::exception& e) {
168 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
169 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(M, true);
170 RCP<tOp> tO = tpOp->getOperator();
171 RCP<tV> diag;
172 if (tO->hasDiagonal()) {
173 diag = rcp(new tV(tO->getRangeMap()));
174 tO->getLocalDiagCopy(*diag);
175 }
176 auto fTpRow = rcp(new MueLu::TpetraOperatorAsRowMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tO, diag));
177 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpFOp = rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fTpRow));
178 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
179 paramList.set<RCP<XpOp> >(parameterName, op);
180 }
181 return true;
182 } else {
183 TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type " << paramList.getEntry(parameterName));
184 return false;
185 }
186 } else
187 return false;
188}
189
190// Constructors/initializers/accessors
191
192template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
193MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
194 : paramList_(rcp(new ParameterList())) {}
195
196template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
197MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() = default;
198
199// Overridden from PreconditionerFactoryBase
200
201template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
202bool MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
203 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
204
205 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp)) return true;
206
207 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isBlockedOperator(fwdOp)) return true;
208
209 return false;
210}
211
212template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
213RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
214 return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
215}
216
217template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
218void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
219 initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
220 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLu::initializePrec")));
221
222 using Teuchos::rcp_dynamic_cast;
223
224 // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
225 typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> XpMap;
226 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
227 typedef MueLu::XpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> MueLuXpOp;
228 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
229 // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
230 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpBlockedCrsMat;
231 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
232 // typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
233 // typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
234 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
235 typedef Thyra::XpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> ThyXpOp;
236 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMV;
237 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
238 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpmMV;
239#if defined(MUELU_CAN_USE_MIXED_PRECISION)
240 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
241 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
242 typedef Xpetra::Operator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfOp;
243 typedef MueLu::XpetraOperator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> MueLuHalfXpOp;
244 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
245 typedef Xpetra::MultiVector<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMV;
246 typedef Xpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> XphmMV;
247 typedef Xpetra::Matrix<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMat;
248#endif
249
250 // Check precondition
251 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
252 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
253 TEUCHOS_ASSERT(prec);
254
255 // Create a copy, as we may remove some things from the list
256 ParameterList paramList = *paramList_;
257
258 // Retrieve wrapped concrete Xpetra matrix from FwdOp
259 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
260 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
261
262 // Check whether it is Tpetra
263 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
264 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
265 TEUCHOS_TEST_FOR_EXCEPT((bIsTpetra == false) && (bIsBlocked == false));
266 TEUCHOS_TEST_FOR_EXCEPT((bIsTpetra == true) && bIsBlocked == true);
267
268 RCP<XpMat> A = Teuchos::null;
269 if (bIsBlocked) {
270 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
271 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
272 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
273
274 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) == false);
275
276 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
277 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
278
279 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
280 // MueLu needs a non-const object as input
281 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
282 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
283
284 RCP<const XpMap> rowmap00 = A00->getRowMap();
285 RCP<const Teuchos::Comm<int> > comm = rowmap00->getComm();
286
287 // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
288 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
289 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
290
291 // save blocked matrix
292 A = bMat;
293 } else {
294 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
295 // MueLu needs a non-const object as input
296 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
297 }
298 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
299
300 // Retrieve concrete preconditioner object
301 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
302 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
303
304 // extract preconditioner operator
305 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
306 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
307
308 // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
309 // rebuild preconditioner if startingOver == true
310 // reuse preconditioner if startingOver == false
311 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
312 bool useHalfPrecision = false;
313 if (paramList.isParameter("half precision"))
314 useHalfPrecision = paramList.get<bool>("half precision");
315 else if (paramList.isSublist("Hierarchy") && paramList.sublist("Hierarchy").isParameter("half precision"))
316 useHalfPrecision = paramList.sublist("Hierarchy").get<bool>("half precision");
317
318 RCP<XpOp> xpPrecOp;
319 if (startingOver == true) {
320 // Convert to Xpetra
321 std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "Material", "BlockNumber"};
322 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
323 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
324
325 if (paramList.isSublist("user data")) {
326 auto& userData = paramList.sublist("user data");
327 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
328 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(userData, *it);
329 }
330
331 for (int lvlNo = 0; lvlNo < 10; ++lvlNo) {
332 if (paramList.isSublist("level " + std::to_string(lvlNo) + " user data")) {
333 ParameterList& lvlList = paramList.sublist("level " + std::to_string(lvlNo) + " user data");
334 std::list<std::string> convertKeys;
335 for (auto it = lvlList.begin(); it != lvlList.end(); ++it)
336 convertKeys.push_back(lvlList.name(it));
337 for (auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
338 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
339 }
340 }
341
342 if (useHalfPrecision) {
343#if defined(MUELU_CAN_USE_MIXED_PRECISION)
344
345 // CAG: There is nothing special about the combination double-float,
346 // except that I feel somewhat confident that Trilinos builds
347 // with both scalar types.
348
349 // convert to half precision
350 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
351 const std::string userName = "user data";
352 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
353 if (userParamList.isType<RCP<XpmMV> >("Coordinates")) {
354 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >("Coordinates");
355 userParamList.remove("Coordinates");
356 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
357 userParamList.set("Coordinates", halfCoords);
358 }
359 if (userParamList.isType<RCP<XpMV> >("Nullspace")) {
360 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >("Nullspace");
361 userParamList.remove("Nullspace");
362 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
363 userParamList.set("Nullspace", halfNullspace);
364 }
365 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
366 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
367 paramList.remove("Coordinates");
368 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
369 userParamList.set("Coordinates", halfCoords);
370 }
371 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
372 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
373 paramList.remove("Nullspace");
374 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
375 userParamList.set("Nullspace", halfNullspace);
376 }
377
378 // build a new half-precision MueLu preconditioner
379
380 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(halfA, paramList);
381 RCP<MueLuHalfXpOp> xpOp = rcp(new MueLuHalfXpOp(H));
382 xpPrecOp = rcp(new XpHalfPrecOp(xpOp));
383#else
384 TEUCHOS_TEST_FOR_EXCEPT(true);
385#endif
386 } else {
387 const std::string userName = "user data";
388 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
389 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
390 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
391 paramList.remove("Coordinates");
392 userParamList.set("Coordinates", coords);
393 }
394 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
395 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
396 paramList.remove("Nullspace");
397 userParamList.set("Nullspace", nullspace);
398 }
399
400 // build a new MueLu RefMaxwell preconditioner
401 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(A, paramList);
402 xpPrecOp = rcp(new MueLuXpOp(H));
403 }
404 } else {
405 // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix
406 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
407 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(), true);
408#if defined(MUELU_CAN_USE_MIXED_PRECISION)
409 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
410 if (!xpHalfPrecOp.is_null()) {
411 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy();
412 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
413
414 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
415 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
416 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
417 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
418 RCP<MueLu::Level> level0 = H->GetLevel(0);
419 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >("A");
420 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0, true);
421
422 if (!A0.is_null()) {
423 // If a user provided a "number of equations" argument in a parameter list
424 // during the initial setup, we must honor that settings and reuse it for
425 // all consequent setups.
426 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
427 }
428
429 // set new matrix
430 level0->Set("A", halfA);
431
432 H->SetupRe();
433 } else
434#endif
435 {
436 // get old MueLu hierarchy
437 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(), true);
438 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
439 ;
440
441 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
442 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
443 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
444 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
445 RCP<MueLu::Level> level0 = H->GetLevel(0);
446 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
447 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
448
449 if (!A0.is_null()) {
450 // If a user provided a "number of equations" argument in a parameter list
451 // during the initial setup, we must honor that settings and reuse it for
452 // all consequent setups.
453 A->SetFixedBlockSize(A0->GetFixedBlockSize());
454 }
455
456 // set new matrix
457 level0->Set("A", A);
458
459 H->SetupRe();
460 }
461 }
462
463 // wrap preconditioner in thyraPrecOp
464 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
465 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
466
467 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
468 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
469
470 defaultPrec->initializeUnspecified(thyraPrecOp);
471}
472
473template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
474void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
475 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
476 TEUCHOS_ASSERT(prec);
477
478 // Retrieve concrete preconditioner object
479 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
480 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
481
482 if (fwdOp) {
483 // TODO: Implement properly instead of returning default value
484 *fwdOp = Teuchos::null;
485 }
486
487 if (supportSolveUse) {
488 // TODO: Implement properly instead of returning default value
489 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
490 }
491
492 defaultPrec->uninitialize();
493}
494
495// Overridden from ParameterListAcceptor
496template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
497void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
498 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
499 paramList_ = paramList;
500}
501
502template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
503RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
504 return paramList_;
505}
506
507template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
508RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
509 RCP<ParameterList> savedParamList = paramList_;
510 paramList_ = Teuchos::null;
511 return savedParamList;
512}
513
514template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
515RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
516 return paramList_;
517}
518
519template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
520RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters() const {
521 static RCP<const ParameterList> validPL;
522
523 if (Teuchos::is_null(validPL)) {
524 validPL = MueLu::MasterList::List();
525 }
526
527 return validPL;
528}
529
530// Public functions overridden from Teuchos::Describable
531template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
532std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
533 return "Thyra::MueLuPreconditionerFactory";
534}
535} // namespace Thyra
536
537#endif // HAVE_MUELU_STRATIMIKOS
538
539#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