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