MueLu  Version of the Day
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"
60 #ifdef HAVE_MUELU_TPETRA
61 #include "Thyra_TpetraLinearOp.hpp"
62 #include "Thyra_TpetraThyraWrappers.hpp"
63 #endif
64 #ifdef HAVE_MUELU_EPETRA
65 #include "Thyra_EpetraLinearOp.hpp"
66 #include "Thyra_EpetraThyraWrappers.hpp"
67 #endif
68 
69 #include "Teuchos_Ptr.hpp"
70 #include "Teuchos_TestForException.hpp"
71 #include "Teuchos_Assert.hpp"
72 #include "Teuchos_Time.hpp"
73 
74 #include <Xpetra_CrsMatrixWrap.hpp>
75 #include <Xpetra_CrsMatrix.hpp>
76 #include <Xpetra_Matrix.hpp>
77 #include <Xpetra_ThyraUtils.hpp>
78 
79 #include <MueLu_XpetraOperator_decl.hpp> // todo fix me
80 #include <MueLu_RefMaxwell.hpp>
81 #ifdef HAVE_MUELU_TPETRA
82 #include <MueLu_TpetraOperator.hpp>
83 #include <Xpetra_TpetraOperator.hpp>
84 #include <Xpetra_TpetraHalfPrecisionOperator.hpp>
85 #endif
86 #ifdef HAVE_MUELU_EPETRA
87 #include <MueLu_EpetraOperator.hpp>
88 #include <Xpetra_EpetraOperator.hpp>
89 #endif
90 
91 #include "Thyra_PreconditionerFactoryBase.hpp"
92 
93 #include "Kokkos_DefaultNode.hpp"
94 
95 #include <list>
96 
97 namespace Thyra {
98 
107  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node = KokkosClassic::DefaultNode::DefaultNodeType>
108  class MueLuRefMaxwellPreconditionerFactory : public PreconditionerFactoryBase<Scalar> {
109  public:
110 
113 
115  MueLuRefMaxwellPreconditionerFactory();
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<Teuchos::ParameterList> paramList_;
166 
167  };
168 
169 #ifdef HAVE_MUELU_EPETRA
170 
177  template <>
178  class MueLuRefMaxwellPreconditionerFactory<double,int,int,Xpetra::EpetraNode> : public PreconditionerFactoryBase<double> {
179  public:
180  typedef double Scalar;
181  typedef int LocalOrdinal;
182  typedef int GlobalOrdinal;
183  typedef Xpetra::EpetraNode Node;
184 
187 
189  MueLuRefMaxwellPreconditionerFactory() : paramList_(rcp(new ParameterList())) { }
190 
192 
195 
197  bool isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
198  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
199 
200 #ifdef HAVE_MUELU_TPETRA
201  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
202 #endif
203 
204 #ifdef HAVE_MUELU_EPETRA
205  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
206 #endif
207 
208  return false;
209  }
210 
212  Teuchos::RCP<PreconditionerBase<Scalar> > createPrec() const {
213  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
214  }
215 
217  void initializePrec(const Teuchos::RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc,
218  PreconditionerBase<Scalar>* prec,
219  const ESupportSolveUse /* supportSolveUse */
220  ) const {
221  using Teuchos::rcp_dynamic_cast;
222 
223  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
224  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
225  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
226  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
227  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
228  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
230 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
231  typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
232  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
233  typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
234  typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
235  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
236  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
237  typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
238  typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
239  typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
240  typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
241 #endif
242  Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLuRefMaxwell::initializePrec")));
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  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
260 
261  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
262  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
263 
264  // MueLu needs a non-const object as input
265  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
266  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
267 
268  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
269  RCP<XpMat> A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
270  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
271 
272  // Retrieve concrete preconditioner object
273  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
274  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
275 
276  // extract preconditioner operator
277  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
278  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
279 
280  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
281  // rebuild preconditioner if startingOver == true
282  // reuse preconditioner if startingOver == false
283  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get<bool>("refmaxwell: enable reuse"));
284  const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;
285 
286  RCP<XpOp> xpPrecOp;
287  if (startingOver == true) {
288 
289  // Convert to Xpetra
290  std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "M1", "Ms", "D0", "M0inv"};
291  for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
292  replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
293 
294  paramList.set<bool>("refmaxwell: use as preconditioner", true);
295  if (useHalfPrecision) {
296 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
297 
298  // convert to half precision
299  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
300  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
301  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
302  paramList.remove("Coordinates");
303  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
304  paramList.set("Coordinates",halfCoords);
305  }
306  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
307  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
308  paramList.remove("Nullspace");
309  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
310  paramList.set("Nullspace",halfNullspace);
311  }
312  std::list<std::string> convertMat = {"M1", "Ms", "D0", "M0inv"};
313  for (auto it = convertMat.begin(); it != convertMat.end(); ++it) {
314  if (paramList.isType<RCP<XpMat> >(*it)) {
315  RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
316  paramList.remove(*it);
317  RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
318  paramList.set(*it,halfM);
319  }
320  }
321 
322  // build a new half-precision MueLu RefMaxwell preconditioner
323  RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > halfPrec = rcp(new MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>(halfA, paramList, true));
324  xpPrecOp = rcp(new XpHalfPrecOp(halfPrec));
325 #else
326  TEUCHOS_TEST_FOR_EXCEPT(true);
327 #endif
328  } else
329  {
330  // build a new MueLu RefMaxwell preconditioner
331  RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp(new MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, paramList, true));
332  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
333  }
334  } else {
335  // reuse old MueLu preconditioner stored in MueLu Xpetra operator and put in new matrix
336 
337  RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
338  RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
339 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
340  RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
341  if (!xpHalfPrecOp.is_null()) {
342  RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(), true);
343  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
344  preconditioner->resetMatrix(halfA);
345  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
346  } else
347 #endif
348  {
349  RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(xpOp, true);
350  preconditioner->resetMatrix(A);
351  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
352  }
353  }
354 
355  // wrap preconditioner in thyraPrecOp
356  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
357  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
358 
359  RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
360  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
361 
362  defaultPrec->initializeUnspecified(thyraPrecOp);
363  }
364 
366  void uninitializePrec(PreconditionerBase<Scalar>* prec,
367  Teuchos::RCP<const LinearOpSourceBase<Scalar> >* fwdOp,
368  ESupportSolveUse* supportSolveUse
369  ) const {
370  TEUCHOS_ASSERT(prec);
371 
372  // Retrieve concrete preconditioner object
373  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
374  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
375 
376  if (fwdOp) {
377  // TODO: Implement properly instead of returning default value
378  *fwdOp = Teuchos::null;
379  }
380 
381  if (supportSolveUse) {
382  // TODO: Implement properly instead of returning default value
383  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
384  }
385 
386  defaultPrec->uninitialize();
387  }
388 
390 
393 
395  void setParameterList(const Teuchos::RCP<Teuchos::ParameterList>& paramList) {
396  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
397  paramList_ = paramList;
398  }
400  Teuchos::RCP<Teuchos::ParameterList> unsetParameterList() {
401  RCP<ParameterList> savedParamList = paramList_;
402  paramList_ = Teuchos::null;
403  return savedParamList;
404  }
406  Teuchos::RCP<Teuchos::ParameterList> getNonconstParameterList() { return paramList_; }
408  Teuchos::RCP<const Teuchos::ParameterList> getParameterList() const { return paramList_; }
410  Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const {
411  static RCP<const ParameterList> validPL;
412 
413  if (Teuchos::is_null(validPL))
414  validPL = rcp(new ParameterList());
415 
416  return validPL;
417  }
419 
422 
424  std::string description() const { return "Thyra::MueLuRefMaxwellPreconditionerFactory"; }
425 
426  // ToDo: Add an override of describe(...) to give more detail!
427 
429 
430  private:
431  Teuchos::RCP<Teuchos::ParameterList> paramList_;
432  }; // end specialization for Epetra
433 
434 #endif // HAVE_MUELU_EPETRA
435 
436 } // namespace Thyra
437 
438 #endif // #ifdef HAVE_MUELU_STRATIMIKOS
439 
440 #endif // THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DECL_HPP
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultNode Node
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
MueLu::DefaultScalar Scalar
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.
Kokkos::Compat::KokkosSerialWrapperNode EpetraNode
void getValidParameters(Teuchos::ParameterList &params)