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 
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 
25 namespace Thyra {
26 
28 using Teuchos::RCP;
29 using Teuchos::rcp;
30 using Teuchos::rcp_const_cast;
31 using Teuchos::rcp_dynamic_cast;
32 
33 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
34 bool Converters<Scalar, LocalOrdinal, GlobalOrdinal, Node>::replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
35  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
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;
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 
53  typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
56 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
57  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
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);
89  RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tM);
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);
95  RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
96  paramList.set<RCP<XpMultVec> >(parameterName, 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);
102  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
103  paramList.set<RCP<XpMagMultVec> >(parameterName, 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);
113  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
114  paramList.set<RCP<XpMagMultVec> >(parameterName, 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());
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  }
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
169 template <class GlobalOrdinal>
170 bool 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;
176  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
183 
184  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
185  typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
186  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
187 
191  typedef Thyra::TpetraVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> thyTpV;
194 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
195  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
197 #endif
198 #if defined(HAVE_MUELU_EPETRA)
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);
230  RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tM);
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);
236  RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
237  paramList.set<RCP<XpMultVec> >(parameterName, 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);
243  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
244  paramList.set<RCP<XpMagMultVec> >(parameterName, 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);
254  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node>(tpetra_X);
255  paramList.set<RCP<XpMagMultVec> >(parameterName, 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());
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  }
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 
345 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
346 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MueLuPreconditionerFactory()
347  : paramList_(rcp(new ParameterList())) {}
348 
349 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
350 MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~MueLuPreconditionerFactory() = default;
351 
352 // Overridden from PreconditionerFactoryBase
353 
354 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
355 bool 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 
369 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
370 RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::createPrec() const {
371  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
372 }
373 
374 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
375 void 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,...)
385  typedef Xpetra::ThyraUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpThyUtils;
386  // typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
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;
394  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
396 #if defined(MUELU_CAN_USE_MIXED_PRECISION)
397  typedef Xpetra::TpetraHalfPrecisionOperator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpHalfPrecOp;
398  typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
401  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
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();
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) {
430  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
432 
433  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0, 0) == false);
434 
435  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0, 0);
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));
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));
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  }
458 
459  // Retrieve concrete preconditioner object
460  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar>*>(prec));
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
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 
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 
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);
624 
625  defaultPrec->initializeUnspecified(thyraPrecOp);
626 }
627 
628 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
629 void 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));
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
651 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
652 void MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::setParameterList(RCP<ParameterList> const& paramList) {
654  paramList_ = paramList;
655 }
656 
657 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
658 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getNonconstParameterList() {
659  return paramList_;
660 }
661 
662 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
663 RCP<ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::unsetParameterList() {
664  RCP<ParameterList> savedParamList = paramList_;
665  paramList_ = Teuchos::null;
666  return savedParamList;
667 }
668 
669 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
670 RCP<const ParameterList> MueLuPreconditionerFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::getParameterList() const {
671  return paramList_;
672 }
673 
674 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
675 RCP<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
685 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
686 std::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
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)
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.