MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_MueLuPreconditionerFactory_decl.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 // Tobias Wiesner (tawiesn@sandia.gov)
43 //
44 // ***********************************************************************
45 //
46 // @HEADER
47 #ifndef THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP
48 #define THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP
49 
50 #include <MueLu_ConfigDefs.hpp>
51 
52 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
53 
54 // Stratimikos needs Thyra, so we don't need special guards for Thyra here
55 #include "Thyra_DefaultPreconditioner.hpp"
56 #include "Thyra_BlockedLinearOpBase.hpp"
57 #include "Thyra_XpetraLinearOp.hpp"
58 #ifdef HAVE_MUELU_TPETRA
59 #include "Thyra_TpetraLinearOp.hpp"
60 #include "Thyra_TpetraThyraWrappers.hpp"
61 #endif
62 #ifdef HAVE_MUELU_EPETRA
63 #include "Thyra_EpetraLinearOp.hpp"
64 #endif
65 
66 #include "Teuchos_Ptr.hpp"
67 #include "Teuchos_TestForException.hpp"
68 #include "Teuchos_Assert.hpp"
69 #include "Teuchos_Time.hpp"
70 
71 #include <Xpetra_CrsMatrixWrap.hpp>
72 #include <Xpetra_CrsMatrix.hpp>
73 #include <Xpetra_Matrix.hpp>
74 #include <Xpetra_ThyraUtils.hpp>
75 
76 #include <MueLu_Hierarchy.hpp>
77 #include <MueLu_HierarchyManager.hpp>
78 #include <MueLu_HierarchyUtils.hpp>
79 #include <MueLu_Utilities.hpp>
80 #include <MueLu_ParameterListInterpreter.hpp>
81 #include <MueLu_MLParameterListInterpreter.hpp>
82 #include <MueLu_MasterList.hpp>
83 #include <MueLu_XpetraOperator_decl.hpp> // todo fix me
84 #include <MueLu_CreateXpetraPreconditioner.hpp>
85 #ifdef HAVE_MUELU_TPETRA
86 #include <MueLu_TpetraOperator.hpp>
87 #endif
88 #ifdef HAVE_MUELU_EPETRA
89 #include <MueLu_EpetraOperator.hpp>
90 #endif
91 
92 #include "Thyra_PreconditionerFactoryBase.hpp"
93 
94 #include "Kokkos_DefaultNode.hpp"
95 
96 
97 namespace Thyra {
98 
107  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node = KokkosClassic::DefaultNode::DefaultNodeType>
108  class MueLuPreconditionerFactory : public PreconditionerFactoryBase<Scalar> {
109  public:
110 
113 
115  MueLuPreconditionerFactory();
117 
120 
122  bool isCompatible(const LinearOpSourceBase<Scalar>& fwdOp) const;
124  Teuchos::RCP<PreconditionerBase<Scalar> > createPrec() const;
126  void initializePrec(const Teuchos::RCP<const LinearOpSourceBase<Scalar> >& fwdOp,
127  PreconditionerBase<Scalar>* prec,
128  const ESupportSolveUse supportSolveUse
129  ) const;
131  void uninitializePrec(PreconditionerBase<Scalar>* prec,
132  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
133  ESupportSolveUse* supportSolveUse
134  ) const;
135 
137 
140 
142  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList);
144  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
146  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
148  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const;
150  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
152 
155 
157  std::string description() const;
158 
159  // ToDo: Add an override of describe(...) to give more detail!
160 
162 
163  private:
164 
165  //Teuchos::RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > CreateXpetraPreconditioner(Teuchos::RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > op, const Teuchos::ParameterList& paramList, Teuchos::RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType, LocalOrdinal, GlobalOrdinal, Node> > coords, Teuchos::RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > nullspace) const;
166 
167  Teuchos::RCP<Teuchos::ParameterList> paramList_;
168 
169  };
170 
171 #ifdef HAVE_MUELU_EPETRA
172 
179  template <>
180  class MueLuPreconditionerFactory<double,int,int,Xpetra::EpetraNode> : public PreconditionerFactoryBase<double> {
181  public:
182  typedef double Scalar;
183  typedef int LocalOrdinal;
184  typedef int GlobalOrdinal;
185  typedef Xpetra::EpetraNode Node;
186 
189 
191  MueLuPreconditionerFactory() : paramList_(rcp(new ParameterList())) { }
192 
194 
197 
199  bool isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
200  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
201 
202 #ifdef HAVE_MUELU_TPETRA
203  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
204 #endif
205 
206 #ifdef HAVE_MUELU_EPETRA
207  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
208 #endif
209 
210  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp)) return true;
211 
212  return false;
213  }
214 
216  Teuchos::RCP<PreconditionerBase<Scalar> > createPrec() const {
217  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
218  }
219 
221  void initializePrec(const Teuchos::RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc,
222  PreconditionerBase<Scalar>* prec,
223  const ESupportSolveUse /* supportSolveUse */
224  ) const {
225  using Teuchos::rcp_dynamic_cast;
226 
227  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
228  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
229  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
230  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
231  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
232  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
233  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
234  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
235  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
236  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
237 #ifdef HAVE_MUELU_TPETRA
238  // TAW 1/26/2016: We deal with Tpetra objects
239 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
240  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
241  typedef MueLu::TpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> MueTpOp;
242  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpOp;
243  typedef Thyra::TpetraLinearOp<Scalar,LocalOrdinal,GlobalOrdinal,Node> ThyTpLinOp;
244 #endif
245 #endif
246 #if defined(HAVE_MUELU_EPETRA)
247  typedef MueLu::EpetraOperator MueEpOp;
248  typedef Thyra::EpetraLinearOp ThyEpLinOp;
249 #endif
250 
251  //std::cout << "-======---------------------------------" << std::endl;
252  //std::cout << *paramList_ << std::endl;
253  //std::cout << "-======---------------------------------" << std::endl;
254 
255  // Check precondition
256  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
257  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
258  TEUCHOS_ASSERT(prec);
259 
260  // Create a copy, as we may remove some things from the list
261  ParameterList paramList = *paramList_;
262 
263  // Retrieve wrapped concrete Xpetra matrix from FwdOp
264  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
265  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
266 
267  // Check whether it is Epetra/Tpetra
268  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
269  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
270  bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
271  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
272  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
273  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
274 
275  RCP<XpMat> A = Teuchos::null;
276  if(bIsBlocked) {
277  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
278  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
279  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
280 
281  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
282 
283  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
284  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
285 
286  RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
287  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
288 
289  // MueLu needs a non-const object as input
290  RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
291  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
292 
293  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
294  RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
295  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
296 
297  RCP<const XpMap> rowmap00 = A00->getRowMap();
298  RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
299 
300  // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
301  RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
302  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
303 
304  // save blocked matrix
305  A = bMat;
306  } else {
307  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
308  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
309 
310  // MueLu needs a non-const object as input
311  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
312  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
313 
314  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
315  A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
316  }
317  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
318 
319  // Retrieve concrete preconditioner object
320  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
321  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
322 
323  // extract preconditioner operator
324  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
325  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
326 
327  // Variable for multigrid hierarchy: either build a new one or reuse the existing hierarchy
328  RCP<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> > H = Teuchos::null;
329 
330  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
331  // rebuild preconditioner if startingOver == true
332  // reuse preconditioner if startingOver == false
333  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
334 
335  if (startingOver == true) {
336  // extract coordinates from parameter list
337  Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
338  coordinates = MueLu::Utilities<Scalar,LocalOrdinal,GlobalOrdinal,Node>::ExtractCoordinatesFromParameterList(paramList);
339 
340  // TODO check for Xpetra or Thyra vectors?
341  RCP<XpMultVec> nullspace = Teuchos::null;
342 #ifdef HAVE_MUELU_TPETRA
343  if (bIsTpetra) {
344 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
345  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
346  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
347  RCP<tMV> tpetra_nullspace = Teuchos::null;
348  if (paramList.isType<Teuchos::RCP<tMV> >("Nullspace")) {
349  tpetra_nullspace = paramList.get<RCP<tMV> >("Nullspace");
350  paramList.remove("Nullspace");
351  nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
352  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
353  }
354 #else
355  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
356  "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
357 #endif
358  }
359 #endif
360 #ifdef HAVE_MUELU_EPETRA
361  if (bIsEpetra) {
362  RCP<Epetra_MultiVector> epetra_nullspace = Teuchos::null;
363  if (paramList.isType<RCP<Epetra_MultiVector> >("Nullspace")) {
364  epetra_nullspace = paramList.get<RCP<Epetra_MultiVector> >("Nullspace");
365  paramList.remove("Nullspace");
366  RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpNullspace = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<int,Node>(epetra_nullspace));
367  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,Node> > xpEpNullspaceMult = rcp_dynamic_cast<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,Node> >(xpEpNullspace);
368  nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult);
369  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
370  }
371  }
372 #endif
373  // build a new MueLu hierarchy
374  const std::string userName = "user data";
375  Teuchos::ParameterList& userParamList = paramList.sublist(userName);
376  if(Teuchos::nonnull(coordinates)) {
377  userParamList.set<RCP<XpMultVecDouble> >("Coordinates", coordinates);
378  }
379  if(Teuchos::nonnull(nullspace)) {
380  userParamList.set<RCP<XpMultVec> >("Nullspace", nullspace);
381  }
382  H = MueLu::CreateXpetraPreconditioner(A, paramList);
383 
384  } else {
385  // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix
386 
387  // get old MueLu hierarchy
388 #if defined(HAVE_MUELU_TPETRA)
389  if (bIsTpetra) {
390 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
391  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
392  RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
393  RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),true);
394 
395  H = muelu_precOp->GetHierarchy();
396 #else
397  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
398  "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
399 #endif
400  }
401 #endif
402 #if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL)
403  if (bIsEpetra) {
404  RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
405  RCP<MueEpOp> muelu_precOp = rcp_dynamic_cast<MueEpOp>(epetr_precOp->epetra_op(),true);
406 
407  H = rcp_dynamic_cast<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_precOp->GetHierarchy());
408  }
409 #endif
410  // TODO add the blocked matrix case here...
411 
412  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
413  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
414  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
415  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
416  RCP<MueLu::Level> level0 = H->GetLevel(0);
417  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
418  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
419 
420  if (!A0.is_null()) {
421  // If a user provided a "number of equations" argument in a parameter list
422  // during the initial setup, we must honor that settings and reuse it for
423  // all consequent setups.
424  A->SetFixedBlockSize(A0->GetFixedBlockSize());
425  }
426 
427  // set new matrix
428  level0->Set("A", A);
429 
430  H->SetupRe();
431  }
432 
433  // wrap hierarchy H in thyraPrecOp
434  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
435 #if defined(HAVE_MUELU_TPETRA)
436  if (bIsTpetra) {
437 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
438  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
439  RCP<MueTpOp> muelu_tpetraOp = rcp(new MueTpOp(H));
440  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
441  RCP<TpOp> tpOp = Teuchos::rcp_dynamic_cast<TpOp>(muelu_tpetraOp);
442  thyraPrecOp = Thyra::createLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpOp);
443 #else
444  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
445  "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
446 #endif
447  }
448 #endif
449 
450 #if defined(HAVE_MUELU_EPETRA)
451  if (bIsEpetra) {
452  RCP<MueLu::Hierarchy<double,int,int,Xpetra::EpetraNode> > epetraH =
453  rcp_dynamic_cast<MueLu::Hierarchy<double,int,int,Xpetra::EpetraNode> >(H);
454  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(epetraH), MueLu::Exceptions::RuntimeError,
455  "Thyra::MueLuPreconditionerFactory: Failed to cast Hierarchy to Hierarchy<double,int,int,Xpetra::EpetraNode>. Epetra runs only on the Serial node.");
456  RCP<MueEpOp> muelu_epetraOp = rcp(new MueEpOp(epetraH));
457  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_epetraOp));
458  // attach fwdOp to muelu_epetraOp to guarantee that it will not go away
459  set_extra_data(fwdOp,"IFPF::fwdOp", Teuchos::inOutArg(muelu_epetraOp), Teuchos::POST_DESTROY,false);
460  RCP<ThyEpLinOp> thyra_epetraOp = Thyra::nonconstEpetraLinearOp(muelu_epetraOp, NOTRANS, EPETRA_OP_APPLY_APPLY_INVERSE, EPETRA_OP_ADJOINT_UNSUPPORTED);
461  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyra_epetraOp));
462  thyraPrecOp = rcp_dynamic_cast<ThyLinOpBase>(thyra_epetraOp);
463  }
464 #endif
465 
466  if(bIsBlocked) {
467  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp));
468 
469  typedef MueLu::XpetraOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> MueXpOp;
470  const RCP<MueXpOp> muelu_xpetraOp = rcp(new MueXpOp(H));
471 
472  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
473  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
474 
475  RCP <Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp = Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_xpetraOp);
476  thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
477  }
478 
479  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
480 
481  defaultPrec->initializeUnspecified(thyraPrecOp);
482  }
483 
485  void uninitializePrec(PreconditionerBase<Scalar>* prec,
486  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
487  ESupportSolveUse* supportSolveUse
488  ) const {
489  TEUCHOS_ASSERT(prec);
490 
491  // Retrieve concrete preconditioner object
492  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
493  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
494 
495  if (fwdOp) {
496  // TODO: Implement properly instead of returning default value
497  *fwdOp = Teuchos::null;
498  }
499 
500  if (supportSolveUse) {
501  // TODO: Implement properly instead of returning default value
502  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
503  }
504 
505  defaultPrec->uninitialize();
506  }
507 
509 
512 
514  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
515  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
516  paramList_ = paramList;
517  }
519  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
520  RCP<ParameterList> savedParamList = paramList_;
521  paramList_ = Teuchos::null;
522  return savedParamList;
523  }
525  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() { return paramList_; }
527  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const { return paramList_; }
529  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
530  static RCP<const ParameterList> validPL;
531 
532  if (Teuchos::is_null(validPL))
533  validPL = rcp(new ParameterList());
534 
535  return validPL;
536  }
538 
541 
543  std::string description() const { return "Thyra::MueLuPreconditionerFactory"; }
544 
545  // ToDo: Add an override of describe(...) to give more detail!
546 
548 
549  private:
550  Teuchos::RCP<Teuchos::ParameterList> paramList_;
551  }; // end specialization for Epetra
552 
553 #endif // HAVE_MUELU_EPETRA
554 
555 } // namespace Thyra
556 
557 #endif // #ifdef HAVE_MUELU_STRATIMIKOS
558 
559 #endif // THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP
MueLu_CreateXpetraPreconditioner.hpp
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
MueLu::Utilities::ExtractCoordinatesFromParameterList
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
Definition: MueLu_Utilities_def.hpp:557
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::nonnull
bool nonnull(const std::shared_ptr< T > &p)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Teuchos_Time.hpp
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
MueLu_XpetraOperator_decl.hpp
MueLu_ConfigDefs.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:52
MueLu_MasterList.hpp
MueLu_HierarchyManager.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Scalar
MueLu::DefaultScalar Scalar
Definition: MueLu_UseDefaultTypes.hpp:49
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:51
MueLu::CreateXpetraPreconditioner
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...
Definition: MueLu_CreateXpetraPreconditioner.hpp:42
MueLu_EpetraOperator.hpp
Teuchos_TestForException.hpp
MueLu::TpetraOperator
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
Definition: MueLu_TpetraOperator_fwd.hpp:54
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
EpetraNode
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Definition: MueLu_ETI_3arg.hpp:22
TEUCHOS_ASSERT
#define TEUCHOS_ASSERT(assertion_test)
Ifpack2::getValidParameters
void getValidParameters(Teuchos::ParameterList &params)
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Thyra_XpetraLinearOp.hpp
MueLu::Hierarchy
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Definition: MueLu_Hierarchy_decl.hpp:103
MueLu::XpetraOperator
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.
Definition: MueLu_XpetraOperator_decl.hpp:64