MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
26 namespace Thyra {
27 
29 using Teuchos::RCP;
30 using Teuchos::rcp;
31 using Teuchos::rcp_const_cast;
32 using Teuchos::rcp_dynamic_cast;
33 
34 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
35 bool Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
36  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
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;
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 
54  typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
57 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
58  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
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 = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
97  paramList.set<RCP<XpMultVec> >(parameterName, 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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
104  paramList.set<RCP<XpMagMultVec> >(parameterName, 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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
115  paramList.set<RCP<XpMagMultVec> >(parameterName, 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());
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
170 template <class GlobalOrdinal>
171 bool 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;
177  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
184 
185  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
186  typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
187  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
188 
192  typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
195 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
196  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
198 #endif
199 #if defined(HAVE_MUELU_EPETRA)
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 = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
238  paramList.set<RCP<XpMultVec> >(parameterName, 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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
245  paramList.set<RCP<XpMagMultVec> >(parameterName, 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 = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
256  paramList.set<RCP<XpMagMultVec> >(parameterName, 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());
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 
346 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
347 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
348  : paramList_(rcp(new ParameterList())) {}
349 
350 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
351 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() = default;
352 
353 // Overridden from PreconditionerFactoryBase
354 
355 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
356 bool 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 
370 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
371 RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
372  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
373 }
374 
375 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
376 void 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,...)
386  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
387  // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
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;
395  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
397 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
398  typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
399  typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
402  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
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();
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) {
431  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
433 
434  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) == false);
435 
436  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
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));
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));
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  }
459 
460  // Retrieve concrete preconditioner object
461  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
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"};
484  for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
485  Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(paramList, *it);
486 
487  for (int lvlNo = 0; lvlNo < 10; ++lvlNo) {
488  if (paramList.isSublist("level " + std::to_string(lvlNo) + " user data")) {
489  ParameterList& lvlList = paramList.sublist("level " + std::to_string(lvlNo) + " user data");
490  std::list<std::string> convertKeys;
491  for (auto it = lvlList.begin(); it != lvlList.end(); ++it)
492  convertKeys.push_back(lvlList.name(it));
493  for (auto it = convertKeys.begin(); it != convertKeys.end(); ++it)
494  Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(lvlList, *it);
495  }
496  }
497 
498  if (useHalfPrecision) {
499 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
500 
501  // CAG: There is nothing special about the combination double-float,
502  // except that I feel somewhat confident that Trilinos builds
503  // with both scalar types.
504 
505  // convert to half precision
506  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
507  const std::string userName = "user data";
508  Teuchos::ParameterList& userParamList = paramList.sublist(userName);
509  if (userParamList.isType<RCP<XpmMV> >("Coordinates")) {
510  RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >("Coordinates");
511  userParamList.remove("Coordinates");
512  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
513  userParamList.set("Coordinates", halfCoords);
514  }
515  if (userParamList.isType<RCP<XpMV> >("Nullspace")) {
516  RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >("Nullspace");
517  userParamList.remove("Nullspace");
518  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
519  userParamList.set("Nullspace", halfNullspace);
520  }
521  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
522  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
523  paramList.remove("Coordinates");
524  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
525  userParamList.set("Coordinates", halfCoords);
526  }
527  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
528  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
529  paramList.remove("Nullspace");
530  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
531  userParamList.set("Nullspace", halfNullspace);
532  }
533 
534  // build a new half-precision MueLu preconditioner
535 
536  RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(halfA, paramList);
537  RCP<MueLuHalfXpOp> xpOp = rcp(new MueLuHalfXpOp(H));
538  xpPrecOp = rcp(new XpHalfPrecOp(xpOp));
539 #else
541 #endif
542  } else {
543  const std::string userName = "user data";
544  Teuchos::ParameterList& userParamList = paramList.sublist(userName);
545  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
546  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
547  paramList.remove("Coordinates");
548  userParamList.set("Coordinates", coords);
549  }
550  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
551  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
552  paramList.remove("Nullspace");
553  userParamList.set("Nullspace", nullspace);
554  }
555 
556  // build a new MueLu RefMaxwell preconditioner
557  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = MueLu::CreateXpetraPreconditioner(A, paramList);
558  xpPrecOp = rcp(new MueLuXpOp(H));
559  }
560  } else {
561  // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix
562  RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
563  xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(), true);
564 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
565  RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
566  if (!xpHalfPrecOp.is_null()) {
567  RCP<MueLu::Hierarchy<HalfScalar, LocalOrdinal, GlobalOrdinal, Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy();
568  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
569 
571  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
572  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
573  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
574  RCP<MueLu::Level> level0 = H->GetLevel(0);
575  RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >("A");
576  RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0, true);
577 
578  if (!A0.is_null()) {
579  // If a user provided a "number of equations" argument in a parameter list
580  // during the initial setup, we must honor that settings and reuse it for
581  // all consequent setups.
582  halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
583  }
584 
585  // set new matrix
586  level0->Set("A", halfA);
587 
588  H->SetupRe();
589  } else
590 #endif
591  {
592  // get old MueLu hierarchy
593  RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(), true);
594  RCP<MueLu::Hierarchy<Scalar, LocalOrdinal, GlobalOrdinal, Node> > H = xpOp->GetHierarchy();
595  ;
596 
598  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
599  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
600  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
601  RCP<MueLu::Level> level0 = H->GetLevel(0);
602  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
603  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
604 
605  if (!A0.is_null()) {
606  // If a user provided a "number of equations" argument in a parameter list
607  // during the initial setup, we must honor that settings and reuse it for
608  // all consequent setups.
609  A->SetFixedBlockSize(A0->GetFixedBlockSize());
610  }
611 
612  // set new matrix
613  level0->Set("A", A);
614 
615  H->SetupRe();
616  }
617  }
618 
619  // wrap preconditioner in thyraPrecOp
620  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getRangeMap());
621  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::toThyra(xpPrecOp->getDomainMap());
622 
623  RCP<ThyLinOpBase> thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
625 
626  defaultPrec->initializeUnspecified(thyraPrecOp);
627 }
628 
629 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
630 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
631  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
632  TEUCHOS_ASSERT(prec);
633 
634  // Retrieve concrete preconditioner object
635  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
637 
638  if (fwdOp) {
639  // TODO: Implement properly instead of returning default value
640  *fwdOp = Teuchos::null;
641  }
642 
643  if (supportSolveUse) {
644  // TODO: Implement properly instead of returning default value
645  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
646  }
647 
648  defaultPrec->uninitialize();
649 }
650 
651 // Overridden from ParameterListAcceptor
652 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
653 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
655  paramList_ = paramList;
656 }
657 
658 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
659 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
660  return paramList_;
661 }
662 
663 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
664 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
665  RCP<ParameterList> savedParamList = paramList_;
666  paramList_ = Teuchos::null;
667  return savedParamList;
668 }
669 
670 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
671 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
672  return paramList_;
673 }
674 
675 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
676 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getValidParameters() const {
677  static RCP<const ParameterList> validPL;
678 
679  if (Teuchos::is_null(validPL)) {
680  validPL = MueLu::MasterList::List();
681  }
682 
683  return validPL;
684 }
685 
686 // Public functions overridden from Teuchos::Describable
687 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
688 std::string MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::description() const {
689  return "Thyra::MueLuPreconditionerFactory";
690 }
691 } // namespace Thyra
692 
693 #endif // HAVE_MUELU_STRATIMIKOS
694 
695 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
bool is_null(const boost::shared_ptr< T > &p)
T & get(const std::string &name, T def_value)
bool nonnull(const std::shared_ptr< T > &p)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::DefaultNode Node
static RCP< Time > getNewTimer(const std::string &name)
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...
void deep_copy(MultiVector< DS, DL, DG, DN > &dst, const MultiVector< SS, SL, SG, SN > &src)
bool remove(std::string const &name, bool throwIfNotExists=true)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
bool isType(const std::string &name) const
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
static RCP< Matrix > Build(const RCP< const Map > &rowMap, size_t maxNumEntriesPerRow, Xpetra::ProfileType pftype=Xpetra::DynamicProfile)
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
static Teuchos::RCP< const Teuchos::ParameterList > List()
Return a &quot;master&quot; list of all valid parameters and their default values.
RCP< Tpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getOperator()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.