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 = 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);
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);
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);
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 = 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);
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);
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);
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", "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
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 
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 
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);
631 
632  defaultPrec->initializeUnspecified(thyraPrecOp);
633 }
634 
635 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
636 void 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));
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
658 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
659 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
661  paramList_ = paramList;
662 }
663 
664 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
665 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
666  return paramList_;
667 }
668 
669 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
670 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
671  RCP<ParameterList> savedParamList = paramList_;
672  paramList_ = Teuchos::null;
673  return savedParamList;
674 }
675 
676 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
677 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
678  return paramList_;
679 }
680 
681 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
682 RCP<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
693 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
694 std::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
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.