MueLu  Version of the Day
Thyra_MueLuPreconditionerFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
47 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
48 
50 
51 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
52 
53 namespace Thyra {
54 
55  using Teuchos::RCP;
56  using Teuchos::rcp;
57  using Teuchos::ParameterList;
58  using Teuchos::rcp_dynamic_cast;
59  using Teuchos::rcp_const_cast;
60 
61 
62  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  static bool replaceWithXpetra(ParameterList& paramList, std::string parameterName) {
64  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
65  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
66  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
67  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
68  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
69  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
70  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
71  typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> XpMagMultVec;
72  typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpVec;
73 
74  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
75  typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
77  typedef Thyra::SpmdVectorSpaceBase<Scalar> ThyVSBase;
78 
79 #ifdef HAVE_MUELU_TPETRA
80  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
81  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
82  typedef Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> thyTpV;
83  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
84  typedef Tpetra::MultiVector<Magnitude, LocalOrdinal, GlobalOrdinal, Node> tMagMV;
85 # if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
86  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
87  typedef Tpetra::MultiVector<HalfMagnitude, LocalOrdinal, GlobalOrdinal, Node> tHalfMagMV;
88 # endif
89 #endif
90 #if defined(HAVE_MUELU_EPETRA)
91  typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
92 #endif
93 
94  if (paramList.isParameter(parameterName)) {
95  if (paramList.isType<RCP<XpMat> >(parameterName))
96  return true;
97  else if (paramList.isType<RCP<const XpMat> >(parameterName)) {
98  RCP<const XpMat> constM = paramList.get<RCP<const XpMat> >(parameterName);
99  paramList.remove(parameterName);
100  RCP<XpMat> M = rcp_const_cast<XpMat>(constM);
101  paramList.set<RCP<XpMat> >(parameterName, M);
102  return true;
103  }
104  else if (paramList.isType<RCP<XpMultVec> >(parameterName))
105  return true;
106  else if (paramList.isType<RCP<const XpMultVec> >(parameterName)) {
107  RCP<const XpMultVec> constX = paramList.get<RCP<const XpMultVec> >(parameterName);
108  paramList.remove(parameterName);
109  RCP<XpMultVec> X = rcp_const_cast<XpMultVec>(constX);
110  paramList.set<RCP<XpMultVec> >(parameterName, X);
111  return true;
112  }
113  else if (paramList.isType<RCP<XpMagMultVec> >(parameterName))
114  return true;
115  else if (paramList.isType<RCP<const XpMagMultVec> >(parameterName)) {
116  RCP<const XpMagMultVec> constX = paramList.get<RCP<const XpMagMultVec> >(parameterName);
117  paramList.remove(parameterName);
118  RCP<XpMagMultVec> X = rcp_const_cast<XpMagMultVec>(constX);
119  paramList.set<RCP<XpMagMultVec> >(parameterName, X);
120  return true;
121  }
122 #ifdef HAVE_MUELU_TPETRA
123  else if (paramList.isType<RCP<TpCrsMat> >(parameterName)) {
124  RCP<TpCrsMat> tM = paramList.get<RCP<TpCrsMat> >(parameterName);
125  paramList.remove(parameterName);
126  RCP<XpMat> xM = MueLu::TpetraCrs_To_XpetraMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tM);
127  paramList.set<RCP<XpMat> >(parameterName, xM);
128  return true;
129  } else if (paramList.isType<RCP<tMV> >(parameterName)) {
130  RCP<tMV> tpetra_X = paramList.get<RCP<tMV> >(parameterName);
131  paramList.remove(parameterName);
132  RCP<XpMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
133  paramList.set<RCP<XpMultVec> >(parameterName, X);
134  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
135  return true;
136  } else if (paramList.isType<RCP<tMagMV> >(parameterName)) {
137  RCP<tMagMV> tpetra_X = paramList.get<RCP<tMagMV> >(parameterName);
138  paramList.remove(parameterName);
139  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
140  paramList.set<RCP<XpMagMultVec> >(parameterName, X);
141  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
142  return true;
143  }
144 # if defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
145  else if (paramList.isType<RCP<tHalfMagMV> >(parameterName)) {
146  RCP<tHalfMagMV> tpetra_hX = paramList.get<RCP<tHalfMagMV> >(parameterName);
147  paramList.remove(parameterName);
148  RCP<tMagMV> tpetra_X = rcp(new tMagMV(tpetra_hX->getMap(),tpetra_hX->getNumVectors()));
149  Tpetra::deep_copy(*tpetra_X,*tpetra_hX);
150  RCP<XpMagMultVec> X = MueLu::TpetraMultiVector_To_XpetraMultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node>(tpetra_X);
151  paramList.set<RCP<XpMagMultVec> >(parameterName, X);
152  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(X));
153  return true;
154  }
155 # endif
156 #endif
157 #ifdef HAVE_MUELU_EPETRA
158  else if (paramList.isType<RCP<Epetra_CrsMatrix> >(parameterName)) {
159  RCP<Epetra_CrsMatrix> eM = paramList.get<RCP<Epetra_CrsMatrix> >(parameterName);
160  paramList.remove(parameterName);
161  RCP<XpEpCrsMat> xeM = rcp(new XpEpCrsMat(eM));
162  RCP<XpCrsMat> xCrsM = rcp_dynamic_cast<XpCrsMat>(xeM, true);
163  RCP<XpCrsMatWrap> xwM = rcp(new XpCrsMatWrap(xCrsM));
164  RCP<XpMat> xM = rcp_dynamic_cast<XpMat>(xwM);
165  paramList.set<RCP<XpMat> >(parameterName, xM);
166  return true;
167  } else if (paramList.isType<RCP<Epetra_MultiVector> >(parameterName)) {
168  RCP<Epetra_MultiVector> epetra_X = Teuchos::null;
169  epetra_X = paramList.get<RCP<Epetra_MultiVector> >(parameterName);
170  paramList.remove(parameterName);
171  RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpX = rcp(new Xpetra::EpetraMultiVectorT<int,Node>(epetra_X));
172  RCP<Xpetra::MultiVector<Scalar,int,int,Node> > xpEpXMult = rcp_dynamic_cast<Xpetra::MultiVector<Scalar,int,int,Node> >(xpEpX, true);
173  RCP<XpMultVec> X = rcp_dynamic_cast<XpMultVec>(xpEpXMult, true);
174  paramList.set<RCP<XpMultVec> >(parameterName, X);
175  return true;
176  }
177 #endif
178  else if (paramList.isType<RCP<const ThyLinOpBase> >(parameterName)) {
179  RCP<const ThyLinOpBase> thyM = paramList.get<RCP<const ThyLinOpBase> >(parameterName);
180  paramList.remove(parameterName);
181  RCP<const XpCrsMat> crsM = XpThyUtils::toXpetra(thyM);
182  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM));
183  // MueLu needs a non-const object as input
184  RCP<XpCrsMat> crsMNonConst = rcp_const_cast<XpCrsMat>(crsM);
185  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMNonConst));
186  // wrap as an Xpetra::Matrix that MueLu can work with
187  RCP<XpMat> M = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsMNonConst));
188  paramList.set<RCP<XpMat> >(parameterName, M);
189  return true;
190  } else if (paramList.isType<RCP<const ThyDiagLinOpBase> >(parameterName)) {
191  RCP<const ThyDiagLinOpBase> thyM = paramList.get<RCP<const ThyDiagLinOpBase> >(parameterName);
192  paramList.remove(parameterName);
193  RCP<const Thyra::VectorBase<Scalar> > diag = thyM->getDiag();
194 
195  RCP<const XpVec> xpDiag;
196 #ifdef HAVE_MUELU_TPETRA
197  if (!rcp_dynamic_cast<const thyTpV>(diag).is_null()) {
198  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
199  if (!tDiag.is_null())
200  xpDiag = Xpetra::toXpetra(tDiag);
201  }
202 #endif
203 #ifdef HAVE_MUELU_EPETRA
204  if (xpDiag.is_null()) {
205  RCP<const Epetra_Comm> comm = Thyra::get_Epetra_Comm(*rcp_dynamic_cast<const ThyVSBase>(thyM->range())->getComm());
206  RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM->range()), comm);
207  if (!map.is_null()) {
208  RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
209  RCP<Epetra_Vector> nceDiag = rcp_const_cast<Epetra_Vector>(eDiag);
210  RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = rcp(new Xpetra::EpetraVectorT<int,Node>(nceDiag));
211  xpDiag = rcp_dynamic_cast<XpVec>(xpEpDiag, true);
212  }
213  }
214 #endif
215  TEUCHOS_ASSERT(!xpDiag.is_null());
216  RCP<XpMat> M = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(xpDiag);
217  paramList.set<RCP<XpMat> >(parameterName, M);
218  return true;
219  }
220  else {
221  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter " << parameterName << " has wrong type.");
222  return false;
223  }
224  } else
225  return false;
226  }
227 
228  // Constructors/initializers/accessors
229 
230  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
231  MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuPreconditionerFactory() :
232  paramList_(rcp(new ParameterList()))
233  {}
234 
235  // Overridden from PreconditionerFactoryBase
236 
237  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
238  bool MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
239  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
240 
241 #ifdef HAVE_MUELU_TPETRA
242  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
243 #endif
244 
245 #ifdef HAVE_MUELU_EPETRA
246  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
247 #endif
248 
249 
250  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp)) return true;
251 
252  return false;
253  }
254 
255 
256  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
257  RCP<PreconditionerBase<Scalar> > MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec() const {
258  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
259  }
260 
261  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
262  void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
263  initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
264  Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLu::initializePrec")));
265 
266  using Teuchos::rcp_dynamic_cast;
267 
268  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
269  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
270  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
272  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
273  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
274  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
275  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
276  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
277  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
278  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
280  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
281  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
282  typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
283 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
284  typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
285  typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
286  typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
288  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
289  typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
290  typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
291  typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
292 #endif
293 
294 
295  // Check precondition
296  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
297  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
298  TEUCHOS_ASSERT(prec);
299 
300  // Create a copy, as we may remove some things from the list
301  ParameterList paramList = *paramList_;
302 
303  // Retrieve wrapped concrete Xpetra matrix from FwdOp
304  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
305  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
306 
307  // Check whether it is Epetra/Tpetra
308  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
309  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
310  bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
311  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
312  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
313  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
314 
315  RCP<XpMat> A = Teuchos::null;
316  if(bIsBlocked) {
317  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
318  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
319  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
320 
321  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
322 
323  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
324  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
325 
326  RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
327  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
328 
329  // MueLu needs a non-const object as input
330  RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
331  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
332 
333  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
334  RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
335  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
336 
337  RCP<const XpMap> rowmap00 = A00->getRowMap();
338  RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
339 
340  // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
341  RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
342  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
343 
344  // save blocked matrix
345  A = bMat;
346  } else {
347  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
348  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
349 
350  // MueLu needs a non-const object as input
351  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
352  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
353 
354  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
355  A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
356  }
357  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
358 
359  // Retrieve concrete preconditioner object
360  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
361  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
362 
363  // extract preconditioner operator
364  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
365  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
366 
367  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
368  // rebuild preconditioner if startingOver == true
369  // reuse preconditioner if startingOver == false
370  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
371  bool useHalfPrecision = false;
372  if (paramList.isParameter("half precision"))
373  useHalfPrecision = paramList.get<bool>("half precision");
374  else if (paramList.isSublist("Hierarchy") && paramList.sublist("Hierarchy").isParameter("half precision"))
375  useHalfPrecision = paramList.sublist("Hierarchy").get<bool>("half precision");
376  if (useHalfPrecision)
377  TEUCHOS_TEST_FOR_EXCEPTION(!bIsTpetra, MueLu::Exceptions::RuntimeError, "The only scalar type Epetra knows is double, so a half precision preconditioner cannot be constructed.");
378 
379  RCP<XpOp> xpPrecOp;
380  if (startingOver == true) {
381  // Convert to Xpetra
382  std::list<std::string> convertXpetra = {"Coordinates", "Nullspace"};
383  for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
384  replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
385 
386  if (useHalfPrecision) {
387 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
388 
389  // CAG: There is nothing special about the combination double-float,
390  // except that I feel somewhat confident that Trilinos builds
391  // with both scalar types.
392 
393  // convert to half precision
394  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
395  const std::string userName = "user data";
396  Teuchos::ParameterList& userParamList = paramList.sublist(userName);
397  if (userParamList.isType<RCP<XpmMV> >("Coordinates")) {
398  RCP<XpmMV> coords = userParamList.get<RCP<XpmMV> >("Coordinates");
399  userParamList.remove("Coordinates");
400  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
401  userParamList.set("Coordinates",halfCoords);
402  }
403  if (userParamList.isType<RCP<XpMV> >("Nullspace")) {
404  RCP<XpMV> nullspace = userParamList.get<RCP<XpMV> >("Nullspace");
405  userParamList.remove("Nullspace");
406  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
407  userParamList.set("Nullspace",halfNullspace);
408  }
409  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
410  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
411  paramList.remove("Coordinates");
412  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
413  userParamList.set("Coordinates",halfCoords);
414  }
415  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
416  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
417  paramList.remove("Nullspace");
418  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
419  userParamList.set("Nullspace",halfNullspace);
420  }
421 
422 
423  // build a new half-precision MueLu preconditioner
424 
425  RCP<MueLu::Hierarchy<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > H = MueLu::CreateXpetraPreconditioner(halfA, paramList);
426  RCP<MueLuHalfXpOp> xpOp = rcp(new MueLuHalfXpOp(H));
427  xpPrecOp = rcp(new XpHalfPrecOp(xpOp));
428 #else
429  TEUCHOS_TEST_FOR_EXCEPT(true);
430 #endif
431  } else
432  {
433  const std::string userName = "user data";
434  Teuchos::ParameterList& userParamList = paramList.sublist(userName);
435  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
436  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
437  paramList.remove("Coordinates");
438  userParamList.set("Coordinates",coords);
439  }
440  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
441  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
442  paramList.remove("Nullspace");
443  userParamList.set("Nullspace",nullspace);
444  }
445 
446  // build a new MueLu RefMaxwell preconditioner
447  RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = MueLu::CreateXpetraPreconditioner(A, paramList);
448  xpPrecOp = rcp(new MueLuXpOp(H));
449  }
450  } else {
451  // reuse old MueLu hierarchy stored in MueLu Xpetra operator and put in new matrix
452  RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
453  xpPrecOp = rcp_dynamic_cast<XpOp>(thyXpOp->getXpetraOperator(), true);
454 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
455  RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpPrecOp);
456  if (!xpHalfPrecOp.is_null()) {
457  RCP<MueLu::Hierarchy<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > H = rcp_dynamic_cast<MueLuHalfXpOp>(xpHalfPrecOp->GetHalfPrecisionOperator(), true)->GetHierarchy();
458  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
459 
460  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
461  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
462  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
463  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
464  RCP<MueLu::Level> level0 = H->GetLevel(0);
465  RCP<XpHalfOp> O0 = level0->Get<RCP<XpHalfOp> >("A");
466  RCP<XphMat> A0 = rcp_dynamic_cast<XphMat>(O0, true);
467 
468  if (!A0.is_null()) {
469  // If a user provided a "number of equations" argument in a parameter list
470  // during the initial setup, we must honor that settings and reuse it for
471  // all consequent setups.
472  halfA->SetFixedBlockSize(A0->GetFixedBlockSize());
473  }
474 
475  // set new matrix
476  level0->Set("A", halfA);
477 
478  H->SetupRe();
479  } else
480 #endif
481  {
482  // get old MueLu hierarchy
483  RCP<MueLuXpOp> xpOp = rcp_dynamic_cast<MueLuXpOp>(thyXpOp->getXpetraOperator(), true);
484  RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = xpOp->GetHierarchy();;
485 
486  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
487  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
488  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
489  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
490  RCP<MueLu::Level> level0 = H->GetLevel(0);
491  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
492  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
493 
494  if (!A0.is_null()) {
495  // If a user provided a "number of equations" argument in a parameter list
496  // during the initial setup, we must honor that settings and reuse it for
497  // all consequent setups.
498  A->SetFixedBlockSize(A0->GetFixedBlockSize());
499  }
500 
501  // set new matrix
502  level0->Set("A", A);
503 
504  H->SetupRe();
505  }
506  }
507 
508  // wrap preconditioner in thyraPrecOp
509  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
510  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
511 
512  RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
513  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
514 
515  defaultPrec->initializeUnspecified(thyraPrecOp);
516 
517  }
518 
519  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
520  void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
521  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
522  TEUCHOS_ASSERT(prec);
523 
524  // Retrieve concrete preconditioner object
525  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
526  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
527 
528  if (fwdOp) {
529  // TODO: Implement properly instead of returning default value
530  *fwdOp = Teuchos::null;
531  }
532 
533  if (supportSolveUse) {
534  // TODO: Implement properly instead of returning default value
535  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
536  }
537 
538  defaultPrec->uninitialize();
539  }
540 
541 
542  // Overridden from ParameterListAcceptor
543  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
544  void MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList> const& paramList) {
545  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
546  paramList_ = paramList;
547  }
548 
549  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
550  RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
551  return paramList_;
552  }
553 
554  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
555  RCP<ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
556  RCP<ParameterList> savedParamList = paramList_;
557  paramList_ = Teuchos::null;
558  return savedParamList;
559  }
560 
561  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
562  RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList() const {
563  return paramList_;
564  }
565 
566  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
567  RCP<const ParameterList> MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters() const {
568  static RCP<const ParameterList> validPL;
569 
570  if (Teuchos::is_null(validPL))
571  validPL = rcp(new ParameterList());
572 
573  return validPL;
574  }
575 
576  // Public functions overridden from Teuchos::Describable
577  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
578  std::string MueLuPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
579  return "Thyra::MueLuPreconditionerFactory";
580  }
581 } // namespace Thyra
582 
583 #endif // HAVE_MUELU_STRATIMIKOS
584 
585 #endif // ifdef THYRA_MUELU_PRECONDITIONER_FACTORY_DEF_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultNode Node
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
Exception throws to report errors in the internal logical of the program.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.