MueLu  Version of the Day
Thyra_MueLuRefMaxwellPreconditionerFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
47 #define THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
48 
50 #include <list>
51 
52 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA)
53 
54 namespace Thyra {
55 
56  using Teuchos::RCP;
57  using Teuchos::rcp;
58  using Teuchos::ParameterList;
59  using Teuchos::rcp_dynamic_cast;
60  using Teuchos::rcp_const_cast;
61 
62  // Constructors/initializers/accessors
63 
64  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuRefMaxwellPreconditionerFactory() :
66  paramList_(rcp(new ParameterList()))
67  {}
68 
69  // Overridden from PreconditionerFactoryBase
70 
71  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
72  bool MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
73  const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
74 
75 #ifdef HAVE_MUELU_TPETRA
76  if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
77 #endif
78 
79  return false;
80  }
81 
82 
83  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
84  RCP<PreconditionerBase<Scalar> > MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec() const {
85  return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
86  }
87 
88  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
89  void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
90  initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
91 
92  // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
93  typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
94  typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
95  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpCrsMat;
96  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
97  typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
99 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
100  typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
101  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
102  typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
103  typedef Xpetra::Operator<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfOp;
104  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
105  typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
106  typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
107  typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
108  typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
109  typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
110 #endif
111  Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLuRefMaxwell::initializePrec")));
112 
113  // Check precondition
114  TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
115  TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
116  TEUCHOS_ASSERT(prec);
117 
118  // Create a copy, as we may remove some things from the list
119  ParameterList paramList = *paramList_;
120 
121  // Retrieve wrapped concrete Xpetra matrix from FwdOp
122  const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
123  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
124 
125  // Check whether it is Epetra/Tpetra
126  bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
127  bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
128  TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
129 
130  RCP<const XpCrsMat > xpetraFwdCrsMat = XpThyUtils::toXpetra(fwdOp);
131  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMat));
132 
133  // MueLu needs a non-const object as input
134  RCP<XpCrsMat> xpetraFwdCrsMatNonConst = Teuchos::rcp_const_cast<XpCrsMat>(xpetraFwdCrsMat);
135  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(xpetraFwdCrsMatNonConst));
136 
137  // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
138  RCP<XpMat> A = rcp(new Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>(xpetraFwdCrsMatNonConst));
139  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
140 
141  // Retrieve concrete preconditioner object
142  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
143  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
144 
145  // extract preconditioner operator
146  RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
147  thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
148 
149  // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
150  // rebuild preconditioner if startingOver == true
151  // reuse preconditioner if startingOver == false
152  const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get<bool>("refmaxwell: enable reuse"));
153  const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;
154 
155  RCP<XpOp> xpPrecOp;
156  if (startingOver == true) {
157 
158  // Convert to Xpetra
159  std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "M1", "Ms", "D0", "M0inv"};
160  for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
161  replaceWithXpetra<Scalar,LocalOrdinal,GlobalOrdinal,Node>(paramList,*it);
162 
163  paramList.set<bool>("refmaxwell: use as preconditioner", true);
164  if (useHalfPrecision) {
165 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
166 
167  // convert to half precision
168  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
169  if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
170  RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
171  paramList.remove("Coordinates");
172  RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
173  paramList.set("Coordinates",halfCoords);
174  }
175  if (paramList.isType<RCP<XpMV> >("Nullspace")) {
176  RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
177  paramList.remove("Nullspace");
178  RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
179  paramList.set("Nullspace",halfNullspace);
180  }
181  std::list<std::string> convertMat = {"M1", "Ms", "D0", "M0inv"};
182  for (auto it = convertMat.begin(); it != convertMat.end(); ++it) {
183  if (paramList.isType<RCP<XpMat> >(*it)) {
184  RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
185  paramList.remove(*it);
186  RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
187  paramList.set(*it,halfM);
188  }
189  }
190 
191  // build a new half-precision MueLu RefMaxwell preconditioner
192  RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > halfPrec = rcp(new MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>(halfA, paramList, true));
193  xpPrecOp = rcp(new XpHalfPrecOp(halfPrec));
194 #else
195  TEUCHOS_TEST_FOR_EXCEPT(true);
196 #endif
197  } else
198  {
199  // build a new MueLu RefMaxwell preconditioner
200  RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp(new MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, paramList, true));
201  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
202  }
203  } else {
204  // reuse old MueLu preconditioner stored in MueLu Xpetra operator and put in new matrix
205 
206  RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
207  RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
208 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT)
209  RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
210  if (!xpHalfPrecOp.is_null()) {
211  RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(), true);
212  RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
213  preconditioner->resetMatrix(halfA);
214  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
215  } else
216 #endif
217  {
218  RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(xpOp, true);
219  preconditioner->resetMatrix(A);
220  xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
221  }
222  }
223 
224  // wrap preconditioner in thyraPrecOp
225  RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
226  RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
227 
228  RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
229  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
230 
231  defaultPrec->initializeUnspecified(thyraPrecOp);
232 
233  }
234 
235  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
236  void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
237  uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
238  TEUCHOS_ASSERT(prec);
239 
240  // Retrieve concrete preconditioner object
241  const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
242  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
243 
244  if (fwdOp) {
245  // TODO: Implement properly instead of returning default value
246  *fwdOp = Teuchos::null;
247  }
248 
249  if (supportSolveUse) {
250  // TODO: Implement properly instead of returning default value
251  *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
252  }
253 
254  defaultPrec->uninitialize();
255  }
256 
257 
258  // Overridden from ParameterListAcceptor
259  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
260  void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList> const& paramList) {
261  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
262  paramList_ = paramList;
263  }
264 
265  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
266  RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
267  return paramList_;
268  }
269 
270  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
271  RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
272  RCP<ParameterList> savedParamList = paramList_;
273  paramList_ = Teuchos::null;
274  return savedParamList;
275  }
276 
277  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
278  RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList() const {
279  return paramList_;
280  }
281 
282  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
283  RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters() const {
284  static RCP<const ParameterList> validPL;
285 
286  if (Teuchos::is_null(validPL))
287  validPL = rcp(new ParameterList());
288 
289  return validPL;
290  }
291 
292  // Public functions overridden from Teuchos::Describable
293  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
294  std::string MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
295  return "Thyra::MueLuRefMaxwellPreconditionerFactory";
296  }
297 } // namespace Thyra
298 
299 #endif // HAVE_MUELU_STRATIMIKOS
300 
301 #endif // ifdef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell&#39;s equations in curl-curl form...
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.