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