MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_MueLuRefMaxwellPreconditionerFactory_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_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP
48 #define THYRA_MUELU_REFMAXWELL_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_DiagonalLinearOpBase.hpp"
58 #include "Thyra_XpetraLinearOp.hpp"
59 #ifdef HAVE_MUELU_TPETRA
60 #include "Thyra_TpetraLinearOp.hpp"
61 #include "Thyra_TpetraThyraWrappers.hpp"
62 #endif
63 #ifdef HAVE_MUELU_EPETRA
64 #include "Thyra_EpetraLinearOp.hpp"
65 #include "Thyra_EpetraThyraWrappers.hpp"
66 #endif
67 
68 #include "Teuchos_Ptr.hpp"
69 #include "Teuchos_TestForException.hpp"
70 #include "Teuchos_Assert.hpp"
71 #include "Teuchos_Time.hpp"
72 
73 #include <Xpetra_CrsMatrixWrap.hpp>
74 #include <Xpetra_CrsMatrix.hpp>
75 #include <Xpetra_Matrix.hpp>
76 #include <Xpetra_ThyraUtils.hpp>
77 
78 #include <MueLu_Hierarchy.hpp>
79 #include <MueLu_HierarchyManager.hpp>
80 #include <MueLu_HierarchyUtils.hpp>
81 #include <MueLu_Utilities.hpp>
82 #include <MueLu_ParameterListInterpreter.hpp>
83 #include <MueLu_MLParameterListInterpreter.hpp>
84 #include <MueLu_MasterList.hpp>
85 #include <MueLu_XpetraOperator_decl.hpp> // todo fix me
86 #include <MueLu_RefMaxwell.hpp>
87 #ifdef HAVE_MUELU_TPETRA
88 #include <MueLu_TpetraOperator.hpp>
89 #endif
90 #ifdef HAVE_MUELU_EPETRA
91 #include <MueLu_EpetraOperator.hpp>
92 #include <Xpetra_EpetraOperator.hpp>
93 #endif
94 
95 #include "Thyra_PreconditionerFactoryBase.hpp"
96 
97 #include "Kokkos_DefaultNode.hpp"
98 
99 
100 namespace Thyra {
101 
110  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node = KokkosClassic::DefaultNode::DefaultNodeType>
111  class MueLuRefMaxwellPreconditionerFactory : public PreconditionerFactoryBase<Scalar> {
112  public:
113 
116 
118  MueLuRefMaxwellPreconditionerFactory();
120 
123 
125  bool isCompatible(const LinearOpSourceBase<Scalar>& fwdOp) const;
127  Teuchos::RCP<PreconditionerBase<Scalar> > createPrec() const;
129  void initializePrec(const Teuchos::RCP<const LinearOpSourceBase<Scalar> >& fwdOp,
130  PreconditionerBase<Scalar>* prec,
131  const ESupportSolveUse supportSolveUse
132  ) const;
134  void uninitializePrec(PreconditionerBase<Scalar>* prec,
135  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
136  ESupportSolveUse* supportSolveUse
137  ) const;
138 
140 
143 
145  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList);
147  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList();
149  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList();
151  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const;
153  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
155 
158 
160  std::string description() const;
161 
162  // ToDo: Add an override of describe(...) to give more detail!
163 
165 
166  private:
167 
168  Teuchos::RCP<Teuchos::ParameterList> paramList_;
169 
170  };
171 
172 #ifdef HAVE_MUELU_EPETRA
173 
180  template <>
181  class MueLuRefMaxwellPreconditionerFactory<double,int,int,Xpetra::EpetraNode> : public PreconditionerFactoryBase<double> {
182  public:
183  typedef double Scalar;
184  typedef int LocalOrdinal;
185  typedef int GlobalOrdinal;
186  typedef Xpetra::EpetraNode Node;
187 
190 
192  MueLuRefMaxwellPreconditionerFactory() : paramList_(rcp(new ParameterList())) { }
193 
195 
198 
200  bool isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
201  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
202 
203 #ifdef HAVE_MUELU_TPETRA
204  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
205 #endif
206 
207 #ifdef HAVE_MUELU_EPETRA
208  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
209 #endif
210 
211  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isBlockedOperator(fwdOp)) return true;
212 
213  return false;
214  }
215 
217  Teuchos::RCP<PreconditionerBase<Scalar> > createPrec() const {
218  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
219  }
220 
222  void initializePrec(const Teuchos::RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc,
223  PreconditionerBase<Scalar>* prec,
224  const ESupportSolveUse /* supportSolveUse */
225  ) const {
226  using Teuchos::rcp_dynamic_cast;
227 
228  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
229  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> XpMap;
230  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
231  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
232  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMatWrap;
233  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
234  typedef Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpBlockedCrsMat;
235  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
236  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMultVec;
237  typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> XpMultVecDouble;
238  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
239  typedef Thyra::DiagonalLinearOpBase<Scalar> ThyDiagLinOpBase;
240 #if defined(HAVE_MUELU_EPETRA)
241  typedef Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> XpEpCrsMat;
242 #endif
243 
244  // Check precondition
245  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
246  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
247  TEUCHOS_ASSERT(prec);
248 
249  // Create a copy, as we may remove some things from the list
250  ParameterList paramList = *paramList_;
251 
252  // Retrieve wrapped concrete Xpetra matrix from FwdOp
253  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
254  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
255 
256  // Check whether it is Epetra/Tpetra
257  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
258  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
259  bool bIsBlocked = XpThyUtils::isBlockedOperator(fwdOp);
260  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
261  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == bIsTpetra) && bIsBlocked == false);
262  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra != bIsTpetra) && bIsBlocked == true);
263 
264  RCP<XpMat> A = Teuchos::null;
265  if(bIsBlocked) {
266  Teuchos::RCP<const Thyra::BlockedLinearOpBase<Scalar> > ThyBlockedOp =
267  Teuchos::rcp_dynamic_cast<const Thyra::BlockedLinearOpBase<Scalar> >(fwdOp);
268  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(ThyBlockedOp));
269 
270  TEUCHOS_TEST_FOR_EXCEPT(ThyBlockedOp->blockExists(0,0)==false);
271 
272  Teuchos::RCP<const LinearOpBase<Scalar> > b00 = ThyBlockedOp->getBlock(0,0);
273  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(b00));
274 
275  RCP<const XpCrsMat > xpetraFwdCrsMat00 = XpThyUtils::toXpetra(b00);
276  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat00));
277 
278  // MueLu needs a non-const object as input
279  RCP<XpCrsMat> xpetraFwdCrsMatNonConst00 = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat00);
280  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst00));
281 
282  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
283  RCP<XpMat> A00 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst00));
284  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A00));
285 
286  RCP<const XpMap> rowmap00 = A00->getRowMap();
287  RCP< const Teuchos::Comm< int > > comm = rowmap00->getComm();
288 
289  // create a Xpetra::BlockedCrsMatrix which derives from Xpetra::Matrix that MueLu can work with
290  RCP<XpBlockedCrsMat> bMat = Teuchos::rcp(new XpBlockedCrsMat(ThyBlockedOp, comm));
291  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(bMat));
292 
293  // save blocked matrix
294  A = bMat;
295  } else {
296  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
297  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
298 
299  // MueLu needs a non-const object as input
300  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
301  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
302 
303  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
304  A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
305  }
306  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
307 
308  // Retrieve concrete preconditioner object
309  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
310  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
311 
312  // extract preconditioner operator
313  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
314  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
315 
316  // Variable for RefMaxwell preconditioner: either build a new one or reuse the existing preconditioner
317  RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = Teuchos::null;
318 
319  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
320  // rebuild preconditioner if startingOver == true
321  // reuse preconditioner if startingOver == false
322  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("reuse: type") || paramList.get<std::string>("reuse: type") == "none");
323 
324  if (startingOver == true) {
325  // extract coordinates from parameter list
326  Teuchos::RCP<XpMultVecDouble> coordinates = Teuchos::null;
327  coordinates = MueLu::Utilities<Scalar,LocalOrdinal,GlobalOrdinal,Node>::ExtractCoordinatesFromParameterList(paramList);
328  paramList.set<RCP<XpMultVecDouble> >("Coordinates", coordinates);
329 
330  // TODO check for Xpetra or Thyra vectors?
331 #ifdef HAVE_MUELU_TPETRA
332  if (bIsTpetra) {
333 #if ((defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT))) || \
334  (!defined(EPETRA_HAVE_OMP) && (defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))))
335  typedef Tpetra::Vector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tV;
336  typedef Tpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> tMV;
337  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> TpCrsMat;
338  if (paramList.isType<Teuchos::RCP<tMV> >("Nullspace")) {
339  RCP<tMV> tpetra_nullspace = paramList.get<RCP<tMV> >("Nullspace");
340  paramList.remove("Nullspace");
341  RCP<XpMultVec> nullspace = MueLu::TpetraMultiVector_To_XpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>(tpetra_nullspace);
342  paramList.set<RCP<XpMultVec> >("Nullspace", nullspace);
343  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(nullspace));
344  }
345 
346  if (paramList.isParameter("M1")) {
347  if (paramList.isType<Teuchos::RCP<TpCrsMat> >("M1")) {
348  RCP<TpCrsMat> tM1 = paramList.get<RCP<TpCrsMat> >("M1");
349  paramList.remove("M1");
350  RCP<XpCrsMat> xM1 = rcp_dynamic_cast<XpCrsMat>(tM1, true);
351  paramList.set<RCP<XpCrsMat> >("M1", xM1);
352  } else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("M1")) {
353  RCP<const ThyLinOpBase> thyM1 = paramList.get<RCP<const ThyLinOpBase> >("M1");
354  paramList.remove("M1");
355  RCP<const XpCrsMat> crsM1 = XpThyUtils::toXpetra(thyM1);
356  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1));
357  // MueLu needs a non-const object as input
358  RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
359  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1NonConst));
360  // wrap as an Xpetra::Matrix that MueLu can work with
361  RCP<XpMat> M1 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
362  paramList.set<RCP<XpMat> >("M1", M1);
363  } else if (paramList.isType<Teuchos::RCP<XpMat> >("M1")) {
364  // do nothing
365  } else
366  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter M1 has wrong type.");
367  } else
368  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix M1.");
369 
370  if (paramList.isParameter("Ms")) {
371  if (paramList.isType<Teuchos::RCP<TpCrsMat> >("Ms")) {
372  RCP<TpCrsMat> tMs = paramList.get<RCP<TpCrsMat> >("Ms");
373  paramList.remove("Ms");
374  RCP<XpCrsMat> xMs = rcp_dynamic_cast<XpCrsMat>(tMs, true);
375  paramList.set<RCP<XpCrsMat> >("Ms", xMs);
376  } else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("Ms")) {
377  RCP<const ThyLinOpBase> thyMs = paramList.get<RCP<const ThyLinOpBase> >("Ms");
378  paramList.remove("Ms");
379  RCP<const XpCrsMat> crsMs = XpThyUtils::toXpetra(thyMs);
380  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMs));
381  // MueLu needs a non-const object as input
382  RCP<XpCrsMat> crsMsNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsMs);
383  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMsNonConst));
384  // wrap as an Xpetra::Matrix that MueLu can work with
385  RCP<XpMat> Ms = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsMsNonConst));
386  paramList.set<RCP<XpMat> >("Ms", Ms);
387  } else if (paramList.isType<Teuchos::RCP<XpMat> >("Ms")) {
388  // do nothing
389  } else
390  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter Ms has wrong type.");
391  }
392 
393  if (paramList.isParameter("D0")) {
394  if(paramList.isType<Teuchos::RCP<TpCrsMat> >("D0")) {
395  RCP<TpCrsMat> tD0 = paramList.get<RCP<TpCrsMat> >("D0");
396  paramList.remove("D0");
397  RCP<XpCrsMat> xD0 = rcp_dynamic_cast<XpCrsMat>(tD0, true);
398  paramList.set<RCP<XpCrsMat> >("D0", xD0);
399  } else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("D0")) {
400  RCP<const ThyLinOpBase> thyD0 = paramList.get<RCP<const ThyLinOpBase> >("D0");
401  paramList.remove("D0");
402  RCP<const XpCrsMat> crsD0 = XpThyUtils::toXpetra(thyD0);
403  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0));
404  // MueLu needs a non-const object as input
405  RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
406  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0NonConst));
407  // wrap as an Xpetra::Matrix that MueLu can work with
408  RCP<XpMat> D0 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
409  paramList.set<RCP<XpMat> >("D0", D0);
410  } else if (paramList.isType<Teuchos::RCP<XpMat> >("D0")) {
411  // do nothing
412  } else
413  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter D0 has wrong type.");
414  } else
415  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix D0.");
416 
417  if (paramList.isParameter("M0inv")) {
418  if (paramList.isType<Teuchos::RCP<TpCrsMat> >("M0inv")) {
419  RCP<TpCrsMat> tM0inv = paramList.get<RCP<TpCrsMat> >("M0inv");
420  paramList.remove("M0inv");
421  RCP<XpCrsMat> xM0inv = rcp_dynamic_cast<XpCrsMat>(tM0inv, true);
422  paramList.set<RCP<XpCrsMat> >("M0inv", xM0inv);
423  } else if (paramList.isType<Teuchos::RCP<const ThyDiagLinOpBase> >("M0inv")) {
424  RCP<const ThyDiagLinOpBase> thyM0inv = paramList.get<RCP<const ThyDiagLinOpBase> >("M0inv");
425  paramList.remove("M0inv");
426  RCP<const Thyra::VectorBase<Scalar> > diag = thyM0inv->getDiag();
427  RCP<const Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > ttDiag = rcp_dynamic_cast<const Thyra::TpetraVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(diag);
428  RCP<const tV> tDiag = Thyra::TpetraOperatorVectorExtraction<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getConstTpetraVector(diag);
429  RCP<XpMat> M0inv = Xpetra::MatrixFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Xpetra::toXpetra(tDiag));
430  paramList.set<RCP<XpMat> >("M0inv", M0inv);
431  } else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("M0inv")) {
432  RCP<const ThyLinOpBase> thyM0inv = paramList.get<RCP<const ThyLinOpBase> >("M0inv");
433  paramList.remove("M0inv");
434  RCP<const XpCrsMat> crsM0inv = XpThyUtils::toXpetra(thyM0inv);
435  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0inv));
436  // MueLu needs a non-const object as input
437  RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
438  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0invNonConst));
439  // wrap as an Xpetra::Matrix that MueLu can work with
440  RCP<XpMat> M0inv = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
441  paramList.set<RCP<XpMat> >("M0inv", M0inv);
442  } else if (paramList.isType<Teuchos::RCP<XpMat> >("M0inv")) {
443  // do nothing
444  } else
445  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter M0inv has wrong type.");
446  } else
447  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix M0inv.");
448 #else
449  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError,
450  "Thyra::MueLuRefMaxwellPreconditionerFactory: Tpetra does not support GO=int and or EpetraNode.");
451 #endif
452  }
453 #endif
454 #ifdef HAVE_MUELU_EPETRA
455  if (bIsEpetra) {
456  if (paramList.isType<RCP<Epetra_MultiVector> >("Nullspace")) {
457  RCP<Epetra_MultiVector> epetra_nullspace = Teuchos::null;
458  epetra_nullspace = paramList.get<RCP<Epetra_MultiVector> >("Nullspace");
459  paramList.remove("Nullspace");
460  RCP<Xpetra::EpetraMultiVectorT<int,Node> > xpEpNullspace = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<int,Node>(epetra_nullspace));
461  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, true);
462  RCP<XpMultVec> nullspace = rcp_dynamic_cast<XpMultVec>(xpEpNullspaceMult, true);
463  paramList.set<RCP<XpMultVec> >("Nullspace", nullspace);
464  }
465 
466  if (paramList.isParameter("M1")) {
467  if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >("M1")) {
468  RCP<Epetra_CrsMatrix> eM1 = paramList.get<RCP<Epetra_CrsMatrix> >("M1");
469  paramList.remove("M1");
470  RCP<XpEpCrsMat> xeM1 = Teuchos::rcp(new XpEpCrsMat(eM1));
471  RCP<XpCrsMat> xCrsM1 = rcp_dynamic_cast<XpCrsMat>(xeM1, true);
472  RCP<XpCrsMatWrap> xwM1 = Teuchos::rcp(new XpCrsMatWrap(xCrsM1));
473  RCP<XpMat> xM1 = rcp_dynamic_cast<XpMat>(xwM1);
474  paramList.set<RCP<XpMat> >("M1", xM1);
475  }
476  else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("M1")) {
477  RCP<const ThyLinOpBase> thyM1 = paramList.get<RCP<const ThyLinOpBase> >("M1");
478  paramList.remove("M1");
479  RCP<const XpCrsMat> crsM1 = XpThyUtils::toXpetra(thyM1);
480  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1));
481  // MueLu needs a non-const object as input
482  RCP<XpCrsMat> crsM1NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM1);
483  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM1NonConst));
484  // wrap as an Xpetra::Matrix that MueLu can work with
485  RCP<XpMat> M1 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM1NonConst));
486  paramList.set<RCP<XpMat> >("M1", M1);
487  } else if (paramList.isType<Teuchos::RCP<XpMat> >("M1")) {
488  // do nothing
489  } else
490  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter M1 has wrong type.");
491  } else
492  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix M1.");
493 
494  if (paramList.isParameter("Ms")) {
495  if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >("Ms")) {
496  RCP<Epetra_CrsMatrix> eMs = paramList.get<RCP<Epetra_CrsMatrix> >("Ms");
497  paramList.remove("Ms");
498  RCP<XpEpCrsMat> xeMs = Teuchos::rcp(new XpEpCrsMat(eMs));
499  RCP<XpCrsMat> xCrsMs = rcp_dynamic_cast<XpCrsMat>(xeMs, true);
500  RCP<XpCrsMatWrap> xwMs = Teuchos::rcp(new XpCrsMatWrap(xCrsMs));
501  RCP<XpMat> xMs = rcp_dynamic_cast<XpMat>(xwMs);
502  paramList.set<RCP<XpMat> >("Ms", xMs);
503  }
504  else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("Ms")) {
505  RCP<const ThyLinOpBase> thyMs = paramList.get<RCP<const ThyLinOpBase> >("Ms");
506  paramList.remove("Ms");
507  RCP<const XpCrsMat> crsMs = XpThyUtils::toXpetra(thyMs);
508  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMs));
509  // MueLu needs a non-const object as input
510  RCP<XpCrsMat> crsMsNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsMs);
511  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsMsNonConst));
512  // wrap as an Xpetra::Matrix that MueLu can work with
513  RCP<XpMat> Ms = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsMsNonConst));
514  paramList.set<RCP<XpMat> >("Ms", Ms);
515  } else if (paramList.isType<Teuchos::RCP<XpMat> >("Ms")) {
516  // do nothing
517  } else
518  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter Ms has wrong type.");
519  }
520 
521  if (paramList.isParameter("D0")) {
522  if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >("D0")) {
523  RCP<Epetra_CrsMatrix> eD0 = paramList.get<RCP<Epetra_CrsMatrix> >("D0");
524  paramList.remove("D0");
525  RCP<XpEpCrsMat> xeD0 = Teuchos::rcp(new XpEpCrsMat(eD0));
526  RCP<XpCrsMat> xCrsD0 = rcp_dynamic_cast<XpCrsMat>(xeD0, true);
527  RCP<XpCrsMatWrap> xwD0 = Teuchos::rcp(new XpCrsMatWrap(xCrsD0));
528  RCP<XpMat> xD0 = rcp_dynamic_cast<XpMat>(xwD0);
529  paramList.set<RCP<XpMat> >("D0", xD0);
530  }
531  else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("D0")) {
532  RCP<const ThyLinOpBase> thyD0 = paramList.get<RCP<const ThyLinOpBase> >("D0");
533  paramList.remove("D0");
534  RCP<const XpCrsMat> crsD0 = XpThyUtils::toXpetra(thyD0);
535  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0));
536  // MueLu needs a non-const object as input
537  RCP<XpCrsMat> crsD0NonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsD0);
538  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsD0NonConst));
539  // wrap as an Xpetra::Matrix that MueLu can work with
540  RCP<XpMat> D0 = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsD0NonConst));
541  paramList.set<RCP<XpMat> >("D0", D0);
542  } else if (paramList.isType<Teuchos::RCP<XpMat> >("D0")) {
543  // do nothing
544  } else
545  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter D0 has wrong type.");
546  } else
547  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix D0.");
548 
549  if (paramList.isParameter("M0inv")) {
550  if (paramList.isType<Teuchos::RCP<Epetra_CrsMatrix> >("M0inv")) {
551  RCP<Epetra_CrsMatrix> eM0inv = paramList.get<RCP<Epetra_CrsMatrix> >("M0inv");
552  paramList.remove("M0inv");
553  RCP<XpEpCrsMat> xeM0inv = Teuchos::rcp(new XpEpCrsMat(eM0inv));
554  RCP<XpCrsMat> xCrsM0inv = rcp_dynamic_cast<XpCrsMat>(xeM0inv, true);
555  RCP<XpCrsMatWrap> xwM0inv = Teuchos::rcp(new XpCrsMatWrap(xCrsM0inv));
556  RCP<XpMat> xM0inv = rcp_dynamic_cast<XpMat>(xwM0inv);
557  paramList.set<RCP<XpMat> >("M0inv", xM0inv);
558  }
559  else if (paramList.isType<Teuchos::RCP<const ThyDiagLinOpBase> >("M0inv")) {
560  RCP<const ThyDiagLinOpBase> thyM0inv = paramList.get<RCP<const ThyDiagLinOpBase> >("M0inv");
561  paramList.remove("M0inv");
562 
563  RCP<const Teuchos::Comm<int> > comm = A->getDomainMap()->getComm();
564  RCP<const Epetra_Map> map = Thyra::get_Epetra_Map(*(thyM0inv->range()), Xpetra::toEpetra(comm));
565  // RCP<XpMap> map = XpThyUtils::toXpetra(thyM0inv->range(), comm);
566  RCP<const Thyra::VectorBase<double> > diag = thyM0inv->getDiag();
567  RCP<const Epetra_Vector> eDiag = Thyra::get_Epetra_Vector(*map, diag);
568  RCP<Epetra_Vector> nceDiag = Teuchos::rcp_const_cast<Epetra_Vector>(eDiag);
569  RCP<Xpetra::EpetraVectorT<int,Node> > xpEpDiag = Teuchos::rcp(new Xpetra::EpetraVectorT<int,Node>(nceDiag));
570  RCP<const Xpetra::Vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,Node> > xpDiag = rcp_dynamic_cast<const Xpetra::Vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,int,int,Node> >(xpEpDiag, true);
571  RCP<XpMat> M0inv = Xpetra::MatrixFactory<double,int,int,Node>::Build(xpDiag);
572  paramList.set<RCP<XpMat> >("M0inv", M0inv);
573  } else if (paramList.isType<Teuchos::RCP<const ThyLinOpBase> >("M0inv")) {
574  RCP<const ThyLinOpBase> thyM0inv = paramList.get<RCP<const ThyLinOpBase> >("M0inv");
575  paramList.remove("M0inv");
576  RCP<const XpCrsMat> crsM0inv = XpThyUtils::toXpetra(thyM0inv);
577  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0inv));
578  // MueLu needs a non-const object as input
579  RCP<XpCrsMat> crsM0invNonConst = Teuchos::rcp_const_cast<XpCrsMat>(crsM0inv);
580  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(crsM0invNonConst));
581  // wrap as an Xpetra::Matrix that MueLu can work with
582  RCP<XpMat> M0inv = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(crsM0invNonConst));
583  paramList.set<RCP<XpMat> >("M0inv", M0inv);
584  } else if (paramList.isType<Teuchos::RCP<XpMat> >("M0inv")) {
585  // do nothing
586  } else
587  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Parameter M0inv has wrong type.");
588  } else
589  TEUCHOS_TEST_FOR_EXCEPTION(true, MueLu::Exceptions::RuntimeError, "Need to specify matrix M0inv.");
590  }
591 #endif
592  // build a new MueLu RefMaxwell preconditioner
593  paramList.set<bool>("refmaxwell: use as preconditioner", true);
594  preconditioner = rcp(new MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, paramList, true));
595 
596  } else {
597  // reuse old MueLu preconditioner stored in MueLu Xpetra operator and put in new matrix
598  preconditioner->resetMatrix(A);
599  }
600 
601  // wrap preconditioner in thyraPrecOp
602  RCP<ThyLinOpBase > thyraPrecOp = Teuchos::null;
603  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getRangeMap());
604  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(preconditioner->getDomainMap());
605 
606  RCP<XpOp> xpOp = Teuchos::rcp_dynamic_cast<XpOp>(preconditioner);
607  thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace,xpOp);
608 
609  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
610 
611  defaultPrec->initializeUnspecified(thyraPrecOp);
612  }
613 
615  void uninitializePrec(PreconditionerBase<Scalar>* prec,
616  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
617  ESupportSolveUse* supportSolveUse
618  ) const {
619  TEUCHOS_ASSERT(prec);
620 
621  // Retrieve concrete preconditioner object
622  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
623  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
624 
625  if (fwdOp) {
626  // TODO: Implement properly instead of returning default value
627  *fwdOp = Teuchos::null;
628  }
629 
630  if (supportSolveUse) {
631  // TODO: Implement properly instead of returning default value
632  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
633  }
634 
635  defaultPrec->uninitialize();
636  }
637 
639 
642 
644  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
645  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
646  paramList_ = paramList;
647  }
649  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
650  RCP<ParameterList> savedParamList = paramList_;
651  paramList_ = Teuchos::null;
652  return savedParamList;
653  }
655  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() { return paramList_; }
657  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const { return paramList_; }
659  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
660  static RCP<const ParameterList> validPL;
661 
662  if (Teuchos::is_null(validPL))
663  validPL = rcp(new ParameterList());
664 
665  return validPL;
666  }
668 
671 
673  std::string description() const { return "Thyra::MueLuRefMaxwellPreconditionerFactory"; }
674 
675  // ToDo: Add an override of describe(...) to give more detail!
676 
678 
679  private:
680  Teuchos::RCP<Teuchos::ParameterList> paramList_;
681  }; // end specialization for Epetra
682 
683 #endif // HAVE_MUELU_EPETRA
684 
685 } // namespace Thyra
686 
687 #endif // #ifdef HAVE_MUELU_STRATIMIKOS
688 
689 #endif // THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP
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::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
Epetra_Vector
MueLu_ConfigDefs.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:52
MueLu_MasterList.hpp
MueLu::RefMaxwell
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
Definition: MueLu_RefMaxwell_fwd.hpp:54
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_EpetraOperator.hpp
Teuchos_TestForException.hpp
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