MueLu  Version of the Day
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 #ifdef HAVE_MUELU_STRATIMIKOS
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>
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
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 
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<double, 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;
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<double,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))))
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;
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))))
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<double,int,int,Node> > xpEpNullspaceMult = rcp_dynamic_cast<Xpetra::MultiVector<double,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  H = MueLu::CreateXpetraPreconditioner(A, paramList, coordinates, nullspace);
375 
376  } else {
377  // reuse old MueLu hierarchy stored in MueLu Tpetra/Epetra operator and put in new matrix
378 
379  // get old MueLu hierarchy
380 #if defined(HAVE_MUELU_TPETRA)
381  if (bIsTpetra) {
382 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
383  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
384  RCP<ThyTpLinOp> tpetr_precOp = rcp_dynamic_cast<ThyTpLinOp>(thyra_precOp);
385  RCP<MueTpOp> muelu_precOp = rcp_dynamic_cast<MueTpOp>(tpetr_precOp->getTpetraOperator(),true);
386 
387  H = muelu_precOp->GetHierarchy();
388 #else
389  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
390  "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
391 #endif
392  }
393 #endif
394 #if defined(HAVE_MUELU_EPETRA)// && defined(HAVE_MUELU_SERIAL)
395  if (bIsEpetra) {
396  RCP<ThyEpLinOp> epetr_precOp = rcp_dynamic_cast<ThyEpLinOp>(thyra_precOp);
397  RCP<MueEpOp> muelu_precOp = rcp_dynamic_cast<MueEpOp>(epetr_precOp->epetra_op(),true);
398 
399  H = rcp_dynamic_cast<MueLu::Hierarchy<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_precOp->GetHierarchy());
400  }
401 #endif
402  // TODO add the blocked matrix case here...
403 
404  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetNumLevels(), MueLu::Exceptions::RuntimeError,
405  "Thyra::MueLuPreconditionerFactory: Hierarchy has no levels in it");
406  TEUCHOS_TEST_FOR_EXCEPTION(!H->GetLevel(0)->IsAvailable("A"), MueLu::Exceptions::RuntimeError,
407  "Thyra::MueLuPreconditionerFactory: Hierarchy has no fine level operator");
408  RCP<MueLu::Level> level0 = H->GetLevel(0);
409  RCP<XpOp> O0 = level0->Get<RCP<XpOp> >("A");
410  RCP<XpMat> A0 = rcp_dynamic_cast<XpMat>(O0);
411 
412  if (!A0.is_null()) {
413  // If a user provided a "number of equations" argument in a parameter list
414  // during the initial setup, we must honor that settings and reuse it for
415  // all consequent setups.
416  A->SetFixedBlockSize(A0->GetFixedBlockSize());
417  }
418 
419  // set new matrix
420  level0->Set("A", A);
421 
422  H->SetupRe();
423  }
424 
425  // wrap hierarchy H in thyraPrecOp
426  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
427 #if defined(HAVE_MUELU_TPETRA)
428  if (bIsTpetra) {
429 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
430  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
431  RCP<MueTpOp> muelu_tpetraOp = rcp(new MueTpOp(H));
432  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_tpetraOp));
433  RCP<TpOp> tpOp = Teuchos::rcp_dynamic_cast<TpOp>(muelu_tpetraOp);
434  thyraPrecOp = Thyra::createLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(tpOp);
435 #else
436  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
437  "Thyra::MueLuPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
438 #endif
439  }
440 #endif
441 
442 #if defined(HAVE_MUELU_EPETRA)
443  if (bIsEpetra) {
444  RCP<MueLu::Hierarchy<double,int,int,Xpetra::EpetraNode> > epetraH =
446  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(epetraH), MueLu::Exceptions::RuntimeError,
447  "Thyra::MueLuPreconditionerFactory: Failed to cast Hierarchy to Hierarchy<double,int,int,Xpetra::EpetraNode>. Epetra runs only on the Serial node.");
448  RCP<MueEpOp> muelu_epetraOp = rcp(new MueEpOp(epetraH));
449  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(muelu_epetraOp));
450  // attach fwdOp to muelu_epetraOp to guarantee that it will not go away
451  set_extra_data(fwdOp,"IFPF::fwdOp", Teuchos::inOutArg(muelu_epetraOp), Teuchos::POST_DESTROY,false);
452  RCP<ThyEpLinOp> thyra_epetraOp = Thyra::nonconstEpetraLinearOp(muelu_epetraOp, NOTRANS, EPETRA_OP_APPLY_APPLY_INVERSE, EPETRA_OP_ADJOINT_UNSUPPORTED);
453  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyra_epetraOp));
454  thyraPrecOp = rcp_dynamic_cast<ThyLinOpBase>(thyra_epetraOp);
455  }
456 #endif
457 
458  if(bIsBlocked) {
459  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::nonnull(thyraPrecOp));
460 
462  const RCP<MueXpOp> muelu_xpetraOp = rcp(new MueXpOp(H));
463 
464  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getRangeMap());
465  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(muelu_xpetraOp->getDomainMap());
466 
467  RCP <Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xpOp = Teuchos::rcp_dynamic_cast<Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(muelu_xpetraOp);
468  thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
469  }
470 
471  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
472 
473  defaultPrec->initializeUnspecified(thyraPrecOp);
474  }
475 
477  void uninitializePrec(PreconditionerBase<Scalar>* prec,
478  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
479  ESupportSolveUse* supportSolveUse
480  ) const {
481  TEUCHOS_ASSERT(prec);
482 
483  // Retrieve concrete preconditioner object
484  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
485  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
486 
487  if (fwdOp) {
488  // TODO: Implement properly instead of returning default value
489  *fwdOp = Teuchos::null;
490  }
491 
492  if (supportSolveUse) {
493  // TODO: Implement properly instead of returning default value
494  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
495  }
496 
497  defaultPrec->uninitialize();
498  }
499 
501 
504 
506  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
507  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
508  paramList_ = paramList;
509  }
511  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
512  RCP<ParameterList> savedParamList = paramList_;
513  paramList_ = Teuchos::null;
514  return savedParamList;
515  }
517  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() { return paramList_; }
519  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const { return paramList_; }
521  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
522  static RCP<const ParameterList> validPL;
523 
524  if (Teuchos::is_null(validPL))
525  validPL = rcp(new ParameterList());
526 
527  return validPL;
528  }
530 
533 
535  std::string description() const { return "Thyra::MueLuPreconditionerFactory"; }
536 
537  // ToDo: Add an override of describe(...) to give more detail!
538 
540 
541  private:
542  Teuchos::RCP<Teuchos::ParameterList> paramList_;
543  }; // end specialization for Epetra
544 
545 #endif // HAVE_MUELU_EPETRA
546 
547 } // namespace Thyra
548 
549 #endif // #ifdef HAVE_MUELU_STRATIMIKOS
550 
551 #endif // THYRA_MUELU_PRECONDITIONER_FACTORY_DECL_HPP
RCP< Level > & GetLevel(const int levelID=0)
Retrieve a certain level from hierarchy.
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords=Teuchos::null, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace=Teuchos::null)
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOp, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOp) const
Teuchos::RCP< Teuchos::ParameterList > paramList_
Concrete preconditioner factory subclass for Thyra based on MueLu.Add support for MueLu preconditione...
void uninitializePrec(PreconditionerBase< Scalar > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< PreconditionerBase< Scalar > > createPrec() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &paramList)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Wraps an existing MueLu::Hierarchy as a Tpetra::Operator.
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, PreconditionerBase< Scalar > *prec, const ESupportSolveUse supportSolveUse) const
Exception throws to report errors in the internal logical of the program.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
Wraps an existing MueLu::Hierarchy as a Xpetra::Operator.