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
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 }
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#ifdef HAVE_MUELU_EPETRA
170template <class GlobalOrdinal>
171bool Converters<double, int, GlobalOrdinal, Tpetra::KokkosCompat::KokkosSerialWrapperNode>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
172 typedef double Scalar;
173 typedef int LocalOrdinal;
174 typedef Tpetra::KokkosCompat::KokkosSerialWrapperNode Node;
175 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
176 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
177 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
178 typedef Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpCrsMatWrap;
179 typedef Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpCrsMat;
180 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
181 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMultVec;
182 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpMagMultVec;
183 typedef Xpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpVec;
184
185 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
186 typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
187 typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
188
189 typedef Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> TpCrsMat;
190 typedef Tpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> tOp;
191 typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
192 typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
193 typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
194 typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
195#if defined(MUELU_CAN_USE_MIXED_PRECISION)
196 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
197 typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
198#endif
199#if defined(HAVE_MUELU_EPETRA)
200 typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal, Node> XpEpCrsMat;
201#endif
202
203 if (paramList.isParameter(parameterName)) {
204 if (paramList.isType<RCP<XpMat> >(parameterName))
205 return true;
206 else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
207 RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
208 paramList.remove(parameterName);
209 RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
210 paramList.set<RCP<XpMat> >(parameterName, M);
211 return true;
212 } else if (paramList.isType<RCP<XpMultVec> >(parameterName))
213 return true;
214 else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
215 RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
216 paramList.remove(parameterName);
217 RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
218 paramList.set<RCP<XpMultVec> >(parameterName, X);
219 return true;
220 } else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
221 return true;
222 else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
223 RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
224 paramList.remove(parameterName);
225 RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
226 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
227 return true;
228 } else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
229 RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
230 paramList.remove(parameterName);
231 RCP<XpMat> xM = Xpetra::toXpetra(tM);
232 paramList.set<RCP<XpMat> >(parameterName, xM);
233 return true;
234 } else if (paramList.isType<RCP<tMV> >(parameterName)) {
235 RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
236 paramList.remove(parameterName);
237 RCP<XpMultVec> X = Xpetra::toXpetra(tpetra_X);
238 paramList.set<RCP<XpMultVec> >(parameterName, X);
239 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
240 return true;
241 } else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
242 RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
243 paramList.remove(parameterName);
244 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
245 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
246 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
247 return true;
248 }
249#if defined(MUELU_CAN_USE_MIXED_PRECISION)
250 else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
251 RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
252 paramList.remove(parameterName);
253 RCP<tMagMV> tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(), tpetra_hX->getNumVectors()));
254 Tpetra::deep_copy(*tpetra_X, *tpetra_hX);
255 RCP<XpMagMultVec> X = Xpetra::toXpetra(tpetra_X);
256 paramList.set<RCP<XpMagMultVec> >(parameterName, X);
257 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
258 return true;
259 }
260#endif
261#ifdef HAVE_MUELU_EPETRA
262 else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
263 RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
264 paramList.remove(parameterName);
265 RCP<XpEpCrsMat> xeM = rcp(new XpEpCrsMat(eM));
266 RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM, true);
267 RCP<XpCrsMatWrap> xwM = rcp(new XpCrsMatWrap(xCrsM));
268 RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
269 paramList.set<RCP<XpMat> >(parameterName, xM);
270 return true;
271 } else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
272 RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
273 epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
274 paramList.remove(parameterName);
275 RCP<Xpetra::EpetraMultiVectorT<int, Node> > xpEpX = rcp(new Xpetra::EpetraMultiVectorT<int, Node>(epetra_X));
276 RCP<Xpetra::MultiVector<Scalar, int, int, Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar, int, int, Node> >(xpEpX, true);
277 RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult, true);
278 paramList.set<RCP<XpMultVec> >(parameterName, X);
279 return true;
280 }
281#endif
282 else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName) ||
283 (paramList.isType<RCP<const ThyLinOpBase> >(parameterName) && !rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName)).is_null())) {
284 RCP<const ThyDiagLinOpBase> thyM;
285 if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName))
286 thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
287 else
288 thyM = rcp_dynamic_cast<const ThyDiagLinOpBase>(paramList.get<RCP<const ThyLinOpBase> >(parameterName), true);
289 paramList.remove(parameterName);
290 RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
291
292 RCP<const XpVec> xpDiag;
293 if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
294 RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getConstTpetraVector(diag);
295 if (!tDiag.is_null())
296 xpDiag = Xpetra::toXpetra(tDiag);
297 }
298#ifdef HAVE_MUELU_EPETRA
299 if (xpDiag.is_null()) {
300 RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
301 RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
302 if (!map.is_null()) {
303 RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
304 RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
305 RCP<Xpetra::EpetraVectorT<int, Node> > xpEpDiag = rcp(new Xpetra::EpetraVectorT<int, Node>(nceDiag));
306 xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag, true);
307 }
308 }
309#endif
310 TEUCHOS_ASSERT(!xpDiag.is_null());
311 RCP<XpMat> M = Xpetra::MatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(xpDiag);
312 paramList.set<RCP<XpMat> >(parameterName, M);
313 return true;
314 } else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
315 RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
316 paramList.remove(parameterName);
317 try {
318 RCP<XpMat> M = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
319 paramList.set<RCP<XpMat> >(parameterName, M);
320 } catch (std::exception& e) {
321 RCP<XpOp> M = XpThyUtils::toXpetraOperator(Teuchos::rcp_const_cast<ThyLinOpBase>(thyM));
322 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpOp = rcp_dynamic_cast<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(M, true);
323 RCP<tOp> tO = tpOp->getOperator();
324 RCP<tV> diag;
325 if (tO->hasDiagonal()) {
326 diag = rcp(new tV(tO->getRangeMap()));
327 tO->getLocalDiagCopy(*diag);
328 }
330 RCP<Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tpFOp = rcp(new Xpetra::TpetraOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node>(fTpRow));
331 auto op = rcp_dynamic_cast<XpOp>(tpFOp);
332 paramList.set<RCP<XpOp> >(parameterName, op);
333 }
334 return true;
335 } else {
336 TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type.");
337 return false;
338 }
339 } else
340 return false;
341}
342#endif
343
344// Constructors/initializers/accessors
345
346template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
347MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
348 : paramList_(rcp(new ParameterList())) {}
349
350template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
351MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() = default;
352
353// Overridden from PreconditionerFactoryBase
354
355template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
356bool MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
357 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
358
359 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isTpetra(fwdOp)) return true;
360
361#ifdef HAVE_MUELU_EPETRA
362 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isEpetra(fwdOp)) return true;
363#endif
364
365 if (Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::isBlockedOperator(fwdOp)) return true;
366
367 return false;
368}
369
370template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
371RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
372 return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
373}
374
375template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
377 initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
378 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLu::initializePrec")));
379
380 using Teuchos::rcp_dynamic_cast;
381
382 // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
383 typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> XpMap;
384 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
386 typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
387 // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
388 typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpBlockedCrsMat;
389 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMat;
390 // typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
391 // typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
392 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
394 typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpMV;
395 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
396 typedef Xpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> XpmMV;
397#if defined(MUELU_CAN_USE_MIXED_PRECISION)
398 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
399 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
400 typedef Xpetra::Operator<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfOp;
402 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
403 typedef Xpetra::MultiVector<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMV;
404 typedef Xpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> XphmMV;
405 typedef Xpetra::Matrix<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> XphMat;
406#endif
407
408 // Check precondition
409 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
410 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
411 TEUCHOS_ASSERT(prec);
412
413 // Create a copy, as we may remove some things from the list
414 ParameterList paramList = *paramList_;
415
416 // Retrieve wrapped concrete Xpetra matrix from FwdOp
417 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
418 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
419
420 // Check whether it is Epetra/Tpetra
421 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
422 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
423 bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
424 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
425 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
426 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
427
428 RCP<XpMat> A = Teuchos::null;
429 if (bIsBlocked) {
430 Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
431 Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
432 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
433
434 TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) == false);
435
436 Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
437 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
438
439 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
440 // MueLu needs a non-const object as input
441 RCP<XpMat> A00 = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(b00));
442 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
443
444 RCP<const XpMap> rowmap00 = A00->getRowMap();
445 RCP<const Teuchos::Comm<int> > comm = rowmap00->getComm();
446
447 // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
448 RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
449 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
450
451 // save blocked matrix
452 A = bMat;
453 } else {
454 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
455 // MueLu needs a non-const object as input
456 A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
457 }
458 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
459
460 // Retrieve concrete preconditioner object
461 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
462 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
463
464 // extract preconditioner operator
465 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
466 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
467
468 // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
469 // rebuild preconditioner if startingOver == true
470 // reuse preconditioner if startingOver == false
471 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
472 bool useHalfPrecision = false;
473 if (paramList.isParameter("half precision"))
474 useHalfPrecision = paramList.get<bool>("half precision");
475 else if (paramList.isSublist("Hierarchy") && paramList.sublist("Hierarchy").isParameter("half precision"))
476 useHalfPrecision = paramList.sublist("Hierarchy").get<bool>("half precision");
477 if (useHalfPrecision)
478 TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra, MueLu::Exceptions::RuntimeError, "The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed.");
479
480 RCP<XpOp> xpPrecOp;
481 if (startingOver == true) {
482 // Convert to Xpetra
483 std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "Material", "BlockNumber"};
484 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
485 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
486
487 if (paramList.isSublist("user data")) {
488 auto& userData = paramList.sublist("user data");
489 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
490 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(userData, *it);
491 }
492
493 for (int lvlNo = 0; lvlNo < 10; ++lvlNo) {
494 if (paramList.isSublist("level " + std::to_string(lvlNo) + " user data")) {
495 ParameterList& lvlList = paramList.sublist("level " + std::to_string(lvlNo) + " user data");
496 std::list<std::string> convertKeys;
497 for (auto it = lvlList.begin(); it != lvlList.end(); ++it)
498 convertKeys.push_back(lvlList.name(it));
499 for (auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
500 Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
501 }
502 }
503
504 if (useHalfPrecision) {
505#if defined(MUELU_CAN_USE_MIXED_PRECISION)
506
507 // CAG: There is nothing special about the combination double-float,
508 // except that I feel somewhat confident that Trilinos builds
509 // with both scalar types.
510
511 // convert to half precision
512 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
513 const std::string userName = "user data";
514 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
515 if (userParamList.isType<RCP<XpmMV> >("Coordinates")) {
516 RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >("Coordinates");
517 userParamList.remove("Coordinates");
518 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
519 userParamList.set("Coordinates", halfCoords);
520 }
521 if (userParamList.isType<RCP<XpMV> >("Nullspace")) {
522 RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >("Nullspace");
523 userParamList.remove("Nullspace");
524 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
525 userParamList.set("Nullspace", halfNullspace);
526 }
527 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
528 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
529 paramList.remove("Coordinates");
530 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
531 userParamList.set("Coordinates", halfCoords);
532 }
533 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
534 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
535 paramList.remove("Nullspace");
536 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
537 userParamList.set("Nullspace", halfNullspace);
538 }
539
540 // build a new half-precision MueLu preconditioner
541
542 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(halfA, paramList);
543 RCP<MueLuHalfXpOp> xpOp = rcp(new MueLuHalfXpOp(H));
544 xpPrecOp = rcp(new XpHalfPrecOp(xpOp));
545#else
546 TEUCHOS_TEST_FOR_EXCEPT(true);
547#endif
548 } else {
549 const std::string userName = "user data";
550 Teuchos::ParameterList& userParamList = paramList.sublist(userName);
551 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
552 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
553 paramList.remove("Coordinates");
554 userParamList.set("Coordinates", coords);
555 }
556 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
557 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
558 paramList.remove("Nullspace");
559 userParamList.set("Nullspace", nullspace);
560 }
561
562 // build a new MueLu RefMaxwell preconditioner
563 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(A, paramList);
564 xpPrecOp = rcp(new MueLuXpOp(H));
565 }
566 } else {
567 // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix
568 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
569 xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(), true);
570#if defined(MUELU_CAN_USE_MIXED_PRECISION)
571 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
572 if (!xpHalfPrecOp.is_null()) {
573 RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy();
574 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
575
576 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
577 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
578 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
579 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
580 RCP<MueLu::Level> level0 = H->GetLevel(0);
581 RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >("A");
582 RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0, true);
583
584 if (!A0.is_null()) {
585 // If a user provided a "number of equations" argument in a parameter list
586 // during the initial setup, we must honor that settings and reuse it for
587 // all consequent setups.
588 halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
589 }
590
591 // set new matrix
592 level0->Set("A", halfA);
593
594 H->SetupRe();
595 } else
596#endif
597 {
598 // get old MueLu hierarchy
599 RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(), true);
600 RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
601 ;
602
603 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
604 "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
605 TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
606 "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
607 RCP<MueLu::Level> level0 = H->GetLevel(0);
608 RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
609 RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
610
611 if (!A0.is_null()) {
612 // If a user provided a "number of equations" argument in a parameter list
613 // during the initial setup, we must honor that settings and reuse it for
614 // all consequent setups.
615 A->SetFixedBlockSize(A0->GetFixedBlockSize());
616 }
617
618 // set new matrix
619 level0->Set("A", A);
620
621 H->SetupRe();
622 }
623 }
624
625 // wrap preconditioner in thyraPrecOp
626 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
627 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
628
629 RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
630 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
631
632 defaultPrec->initializeUnspecified(thyraPrecOp);
633}
634
635template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
636void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
637 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
638 TEUCHOS_ASSERT(prec);
639
640 // Retrieve concrete preconditioner object
641 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
642 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
643
644 if (fwdOp) {
645 // TODO: Implement properly instead of returning default value
646 *fwdOp = Teuchos::null;
647 }
648
649 if (supportSolveUse) {
650 // TODO: Implement properly instead of returning default value
651 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
652 }
653
654 defaultPrec->uninitialize();
655}
656
657// Overridden from ParameterListAcceptor
658template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
659void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
660 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
661 paramList_ = paramList;
662}
663
664template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
665RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
666 return paramList_;
667}
668
669template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
670RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
671 RCP<ParameterList> savedParamList = paramList_;
672 paramList_ = Teuchos::null;
673 return savedParamList;
674}
675
676template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
677RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
678 return paramList_;
679}
680
681template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
682RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters() const {
683 static RCP<const ParameterList> validPL;
684
685 if (Teuchos::is_null(validPL)) {
686 validPL = MueLu::MasterList::List();
687 }
688
689 return validPL;
690}
691
692// Public functions overridden from Teuchos::Describable
693template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
694std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
695 return "Thyra::MueLuPreconditionerFactory";
696}
697} // namespace Thyra
698
699#endif // HAVE_MUELU_STRATIMIKOS
700
701#endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultNode Node
Exception throws to report errors in the internal logical of the program.
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a "master" list of all valid parameters and their default values.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
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,...