MueLu Version of the Day
Loading...
Searching...
No Matches
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// This is not as general as possible, but should be good enough for most builds.
55#if (defined(HAVE_MUELU_TPETRA) && \
56 ((defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && !defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && !defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
57 (!defined(HAVE_TPETRA_INST_DOUBLE) && !defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)) || \
58 (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT) && defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT))))
59# define MUELU_CAN_USE_MIXED_PRECISION
60#endif
61
62
63namespace Thyra {
64
65 using Teuchos::RCP;
66 using Teuchos::rcp;
67 using Teuchos::ParameterList;
68 using Teuchos::rcp_dynamic_cast;
69 using Teuchos::rcp_const_cast;
70
71 // Constructors/initializers/accessors
72
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MueLuRefMaxwellPreconditionerFactory() :
75 paramList_(rcp(new ParameterList()))
76 {}
77
78 // Overridden from PreconditionerFactoryBase
79
80 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
81 bool MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isCompatible(const LinearOpSourceBase<Scalar>& fwdOpSrc) const {
82 const RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc.getOp();
83
84#ifdef HAVE_MUELU_TPETRA
85 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isTpetra(fwdOp)) return true;
86#endif
87
88#ifdef HAVE_MUELU_EPETRA
89 if (Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::isEpetra(fwdOp)) return true;
90#endif
91
92 return false;
93 }
94
95
96 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97 RCP<PreconditionerBase<Scalar> > MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::createPrec() const {
98 return Teuchos::rcp(new DefaultPreconditioner<Scalar>);
99 }
100
101 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
102 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
103 initializePrec(const RCP<const LinearOpSourceBase<Scalar> >& fwdOpSrc, PreconditionerBase<Scalar>* prec, const ESupportSolveUse supportSolveUse) const {
104
105 // we are using typedefs here, since we are using objects from different packages (Xpetra, Thyra,...)
106 typedef Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> XpOp;
107 typedef Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpThyUtils;
108 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMat;
109 typedef Thyra::LinearOpBase<Scalar> ThyLinOpBase;
111#if defined(MUELU_CAN_USE_MIXED_PRECISION)
112 typedef Xpetra::TpetraHalfPrecisionOperator<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpHalfPrecOp;
113 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> XpMV;
114 typedef typename XpHalfPrecOp::HalfScalar HalfScalar;
115 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
116 typedef typename Teuchos::ScalarTraits<Magnitude>::halfPrecision HalfMagnitude;
117 typedef Xpetra::MultiVector<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMV;
118 typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node > XpmMV;
119 typedef Xpetra::MultiVector<HalfMagnitude,LocalOrdinal,GlobalOrdinal,Node > XphmMV;
120 typedef Xpetra::Matrix<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> XphMat;
121#endif
122 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string("ThyraMueLuRefMaxwell::initializePrec")));
123
124 // Check precondition
125 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
126 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
127 TEUCHOS_ASSERT(prec);
128
129 // Create a copy, as we may remove some things from the list
130 ParameterList paramList = *paramList_;
131
132 // Retrieve wrapped concrete Xpetra matrix from FwdOp
133 const RCP<const ThyLinOpBase> fwdOp = fwdOpSrc->getOp();
134 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
135
136 // Check whether it is Epetra/Tpetra
137 bool bIsEpetra = XpThyUtils::isEpetra(fwdOp);
138 bool bIsTpetra = XpThyUtils::isTpetra(fwdOp);
139 TEUCHOS_TEST_FOR_EXCEPT((bIsEpetra == true && bIsTpetra == true));
140
141 // wrap the forward operator as an Xpetra::Matrix that MueLu can work with
142 // MueLu needs a non-const object as input
143 RCP<XpMat> A = XpThyUtils::toXpetra(Teuchos::rcp_const_cast<ThyLinOpBase>(fwdOp));
144 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(A));
145
146 // Retrieve concrete preconditioner object
147 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
148 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
149
150 // extract preconditioner operator
151 RCP<ThyLinOpBase> thyra_precOp = Teuchos::null;
152 thyra_precOp = rcp_dynamic_cast<Thyra::LinearOpBase<Scalar> >(defaultPrec->getNonconstUnspecifiedPrecOp(), true);
153
154 // make a decision whether to (re)build the multigrid preconditioner or reuse the old one
155 // rebuild preconditioner if startingOver == true
156 // reuse preconditioner if startingOver == false
157 const bool startingOver = (thyra_precOp.is_null() || !paramList.isParameter("refmaxwell: enable reuse") || !paramList.get<bool>("refmaxwell: enable reuse"));
158 const bool useHalfPrecision = paramList.get<bool>("half precision", false) && bIsTpetra;
159
160 RCP<XpOp> xpPrecOp;
161 if (startingOver == true) {
162
163 // Convert to Xpetra
164 std::list<std::string> convertXpetra = {"Coordinates", "Nullspace", "M1", "Ms", "D0", "M0inv"};
165 for (auto it = convertXpetra.begin(); it != convertXpetra.end(); ++it)
166 Converters<Scalar,LocalOrdinal,GlobalOrdinal,Node>::replaceWithXpetra(paramList,*it);
167
168 paramList.set<bool>("refmaxwell: use as preconditioner", true);
169 if (useHalfPrecision) {
170#if defined(MUELU_CAN_USE_MIXED_PRECISION)
171
172 // convert to half precision
173 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
174 if (paramList.isType<RCP<XpmMV> >("Coordinates")) {
175 RCP<XpmMV> coords = paramList.get<RCP<XpmMV> >("Coordinates");
176 paramList.remove("Coordinates");
177 RCP<XphmMV> halfCoords = Xpetra::convertToHalfPrecision(coords);
178 paramList.set("Coordinates",halfCoords);
179 }
180 if (paramList.isType<RCP<XpMV> >("Nullspace")) {
181 RCP<XpMV> nullspace = paramList.get<RCP<XpMV> >("Nullspace");
182 paramList.remove("Nullspace");
183 RCP<XphMV> halfNullspace = Xpetra::convertToHalfPrecision(nullspace);
184 paramList.set("Nullspace",halfNullspace);
185 }
186 std::list<std::string> convertMat = {"M1", "Ms", "D0", "M0inv"};
187 for (auto it = convertMat.begin(); it != convertMat.end(); ++it) {
188 if (paramList.isType<RCP<XpMat> >(*it)) {
189 RCP<XpMat> M = paramList.get<RCP<XpMat> >(*it);
190 paramList.remove(*it);
191 RCP<XphMat> halfM = Xpetra::convertToHalfPrecision(M);
192 paramList.set(*it,halfM);
193 }
194 }
195
196 // build a new half-precision MueLu RefMaxwell preconditioner
197 RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > halfPrec = rcp(new MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>(halfA, paramList, true));
198 xpPrecOp = rcp(new XpHalfPrecOp(halfPrec));
199#else
200 TEUCHOS_TEST_FOR_EXCEPT(true);
201#endif
202 } else
203 {
204 // build a new MueLu RefMaxwell preconditioner
205 RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp(new MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>(A, paramList, true));
206 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
207 }
208 } else {
209 // reuse old MueLu preconditioner stored in MueLu Xpetra operator and put in new matrix
210
211 RCP<ThyXpOp> thyXpOp = rcp_dynamic_cast<ThyXpOp>(thyra_precOp, true);
212 RCP<XpOp> xpOp = thyXpOp->getXpetraOperator();
213#if defined(MUELU_CAN_USE_MIXED_PRECISION)
214 RCP<XpHalfPrecOp> xpHalfPrecOp = rcp_dynamic_cast<XpHalfPrecOp>(xpOp);
215 if (!xpHalfPrecOp.is_null()) {
216 RCP<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<HalfScalar,LocalOrdinal,GlobalOrdinal,Node>>(xpHalfPrecOp->GetHalfPrecisionOperator(), true);
217 RCP<XphMat> halfA = Xpetra::convertToHalfPrecision(A);
218 preconditioner->resetMatrix(halfA);
219 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
220 } else
221#endif
222 {
223 RCP<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node> > preconditioner = rcp_dynamic_cast<MueLu::RefMaxwell<Scalar,LocalOrdinal,GlobalOrdinal,Node>>(xpOp, true);
224 preconditioner->resetMatrix(A);
225 xpPrecOp = rcp_dynamic_cast<XpOp>(preconditioner);
226 }
227 }
228
229 // wrap preconditioner in thyraPrecOp
230 RCP<const VectorSpaceBase<Scalar> > thyraRangeSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getRangeMap());
231 RCP<const VectorSpaceBase<Scalar> > thyraDomainSpace = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyra(xpPrecOp->getDomainMap());
232
233 RCP<ThyLinOpBase > thyraPrecOp = Thyra::xpetraLinearOp<Scalar, LocalOrdinal, GlobalOrdinal, Node>(thyraRangeSpace, thyraDomainSpace, xpPrecOp);
234 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(thyraPrecOp));
235
236 defaultPrec->initializeUnspecified(thyraPrecOp);
237
238 }
239
240 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
241 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
242 uninitializePrec(PreconditionerBase<Scalar>* prec, RCP<const LinearOpSourceBase<Scalar> >* fwdOp, ESupportSolveUse* supportSolveUse) const {
243 TEUCHOS_ASSERT(prec);
244
245 // Retrieve concrete preconditioner object
246 const Teuchos::Ptr<DefaultPreconditioner<Scalar> > defaultPrec = Teuchos::ptr(dynamic_cast<DefaultPreconditioner<Scalar> *>(prec));
247 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
248
249 if (fwdOp) {
250 // TODO: Implement properly instead of returning default value
251 *fwdOp = Teuchos::null;
252 }
253
254 if (supportSolveUse) {
255 // TODO: Implement properly instead of returning default value
256 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
257 }
258
259 defaultPrec->uninitialize();
260 }
261
262
263 // Overridden from ParameterListAcceptor
264 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
265 void MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::setParameterList(RCP<ParameterList> const& paramList) {
266 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
267 paramList_ = paramList;
268 }
269
270 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
271 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getNonconstParameterList() {
272 return paramList_;
273 }
274
275 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
276 RCP<ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::unsetParameterList() {
277 RCP<ParameterList> savedParamList = paramList_;
278 paramList_ = Teuchos::null;
279 return savedParamList;
280 }
281
282 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
283 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getParameterList() const {
284 return paramList_;
285 }
286
287 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
288 RCP<const ParameterList> MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::getValidParameters() const {
289 static RCP<const ParameterList> validPL;
290
291 if (Teuchos::is_null(validPL))
292 validPL = rcp(new ParameterList());
293
294 return validPL;
295 }
296
297 // Public functions overridden from Teuchos::Describable
298 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
299 std::string MueLuRefMaxwellPreconditionerFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::description() const {
300 return "Thyra::MueLuRefMaxwellPreconditionerFactory";
301 }
302} // namespace Thyra
303
304#endif // HAVE_MUELU_STRATIMIKOS
305
306#endif // ifdef THYRA_MUELU_REFMAXWELL_PRECONDITIONER_FACTORY_DEF_HPP
Preconditioner (wrapped as a Xpetra::Operator) for Maxwell's equations in curl-curl form.
Concrete Thyra::LinearOpBase subclass for Xpetra::Operator.