Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_BelosTpetraPreconditionerFactory_def.hpp
Go to the documentation of this file.
1/*@HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Jennifer A. Loe (jloe@sandia.gov)
38//
39// ***********************************************************************
40//@HEADER
41*/
42#ifndef THYRA_BELOS_TPETRA_PRECONDITIONERFACTORY_DEF_HPP
43#define THYRA_BELOS_TPETRA_PRECONDITIONERFACTORY_DEF_HPP
44
46
47#include "Thyra_DefaultPreconditioner.hpp"
48#include "Thyra_TpetraLinearOp.hpp"
49#include "Thyra_TpetraThyraWrappers.hpp"
50
51#include "BelosTpetraOperator.hpp"
52#include "BelosConfigDefs.hpp"
53#include "BelosLinearProblem.hpp"
54#include "BelosTpetraAdapter.hpp"
55
56#include "Tpetra_MixedScalarMultiplyOp.hpp"
57
58#include "Teuchos_TestForException.hpp"
59#include "Teuchos_Assert.hpp"
60#include "Teuchos_Time.hpp"
61#include "Teuchos_FancyOStream.hpp"
62#include "Teuchos_VerbosityLevel.hpp"
63#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
64
65#include <string>
66
67
68// CAG: This is not entirely ideal, since it disables half-precision
69// altogether when we have e.g. double, float and
70// complex<double>, but not complex<float>, although a
71// double-float preconditioner would be possible.
72#if (!defined(HAVE_TPETRA_INST_DOUBLE) || (defined(HAVE_TPETRA_INST_DOUBLE) && defined(HAVE_TPETRA_INST_FLOAT))) && \
73 (!defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) || (defined(HAVE_TPETRA_INST_COMPLEX_DOUBLE) && defined(HAVE_TPETRA_INST_COMPLEX_FLOAT)))
74# define THYRA_BELOS_PREC_ENABLE_HALF_PRECISION
75#endif
76
77namespace Thyra {
78
79// Constructors/initializers/accessors
80
81
82template <typename MatrixType>
85
86
87// Overridden from PreconditionerFactoryBase
88
89
90template <typename MatrixType>
92 const LinearOpSourceBase<scalar_type> &fwdOpSrc
93 ) const
94{
95 typedef typename MatrixType::local_ordinal_type local_ordinal_type;
96 typedef typename MatrixType::global_ordinal_type global_ordinal_type;
97 typedef typename MatrixType::node_type node_type;
98
99 const Teuchos::RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc.getOp();
100 using TpetraExtractHelper = TpetraOperatorVectorExtraction<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
101 const auto tpetraFwdOp = TpetraExtractHelper::getConstTpetraOperator(fwdOp);
102
103 const Teuchos::RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp,true);
104
105 return Teuchos::nonnull(tpetraFwdMatrix);
106}
107
108
109template <typename MatrixType>
110Teuchos::RCP<PreconditionerBase<typename BelosTpetraPreconditionerFactory<MatrixType>::scalar_type> >
112{
113 return Teuchos::rcp(new DefaultPreconditioner<scalar_type>);
114}
115
116
117template <typename MatrixType>
119 const Teuchos::RCP<const LinearOpSourceBase<scalar_type> > &fwdOpSrc,
120 PreconditionerBase<scalar_type> *prec,
121 const ESupportSolveUse /* supportSolveUse */
122 ) const
123{
124 using Teuchos::rcp;
125 using Teuchos::RCP;
126
127 // Check precondition
128
129 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
130 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
131 TEUCHOS_ASSERT(prec);
132
133 Teuchos::Time totalTimer(""), timer("");
134 totalTimer.start(true);
135
136 const RCP<Teuchos::FancyOStream> out = this->getOStream();
137 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
138 Teuchos::OSTab tab(out);
139 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
140 *out << "\nEntering Thyra::BelosTpetraPreconditionerFactory::initializePrec(...) ...\n";
141 }
142
143 // Retrieve wrapped concrete Tpetra matrix from FwdOp
144
145 const RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc->getOp();
146 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
147
148 typedef typename MatrixType::local_ordinal_type local_ordinal_type;
149 typedef typename MatrixType::global_ordinal_type global_ordinal_type;
150 typedef typename MatrixType::node_type node_type;
151
152 typedef Tpetra::Operator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraLinOp;
153 using TpetraExtractHelper = TpetraOperatorVectorExtraction<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
154 const auto tpetraFwdOp = TpetraExtractHelper::getConstTpetraOperator(fwdOp);
155 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdOp));
156
157 // Belos-specific typedefs:
158 typedef Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraMV;
159 typedef Belos::TpetraOperator<scalar_type, local_ordinal_type, global_ordinal_type, node_type> BelosTpOp;
160 typedef Belos::LinearProblem<scalar_type, TpetraMV, TpetraLinOp> BelosTpLinProb;
161
162#ifdef THYRA_BELOS_PREC_ENABLE_HALF_PRECISION
163 // CAG: There is nothing special about the combination double-float,
164 // except that I feel somewhat confident that Trilinos builds
165 // with both scalar types.
166 typedef typename Teuchos::ScalarTraits<scalar_type>::halfPrecision half_scalar_type;
167 typedef Tpetra::Operator<half_scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraLinOpHalf;
168 typedef Tpetra::MultiVector<half_scalar_type, local_ordinal_type, global_ordinal_type, node_type> TpetraMVHalf;
169 typedef Belos::TpetraOperator<half_scalar_type, local_ordinal_type, global_ordinal_type, node_type> BelosTpOpHalf;
170 typedef Belos::LinearProblem<half_scalar_type, TpetraMVHalf, TpetraLinOpHalf> BelosTpLinProbHalf;
171#endif
172
173 // Retrieve concrete preconditioner object
174
175 const Teuchos::Ptr<DefaultPreconditioner<scalar_type> > defaultPrec =
176 Teuchos::ptr(dynamic_cast<DefaultPreconditioner<scalar_type> *>(prec));
177 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
178
179 // This check needed to address Issue #535.
180 RCP<Teuchos::ParameterList> innerParamList;
181 if (paramList_.is_null ()) {
182 innerParamList = rcp(new Teuchos::ParameterList(*getValidParameters()));
183 }
184 else {
185 innerParamList = paramList_;
186 }
187
188 bool useHalfPrecision = false;
189 if (innerParamList->isParameter("half precision"))
190 useHalfPrecision = Teuchos::getParameter<bool>(*innerParamList, "half precision");
191
192 const std::string solverType = Teuchos::getParameter<std::string>(*innerParamList, "BelosPrec Solver Type");
193 const RCP<Teuchos::ParameterList> packageParamList = Teuchos::rcpFromRef(innerParamList->sublist("BelosPrec Solver Params"));
194
195 // solverTypeUpper is the upper-case version of solverType.
196 std::string solverTypeUpper (solverType);
197 std::transform(solverTypeUpper.begin(), solverTypeUpper.end(),solverTypeUpper.begin(), ::toupper);
198
199 // Create the initial preconditioner
200 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
201 *out << "\nCreating a new BelosTpetra::Preconditioner object...\n";
202 }
203 RCP<LinearOpBase<scalar_type> > thyraPrecOp;
204
205 if (useHalfPrecision) {
206#ifdef THYRA_BELOS_PREC_ENABLE_HALF_PRECISION
207 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
208 Teuchos::OSTab(out).o() << "> Creating half precision preconditioner\n";
209 }
210 const RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp);
211 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(tpetraFwdMatrix));
212 auto tpetraFwdMatrixHalf = tpetraFwdMatrix->template convert<half_scalar_type>();
213
214 RCP<BelosTpLinProbHalf> belosLinProbHalf = rcp(new BelosTpLinProbHalf());
215 belosLinProbHalf->setOperator(tpetraFwdMatrixHalf);
216 RCP<TpetraLinOpHalf> belosOpRCPHalf = rcp(new BelosTpOpHalf(belosLinProbHalf, packageParamList, solverType, true));
217 RCP<TpetraLinOp> wrappedOp = rcp(new Tpetra::MixedScalarMultiplyOp<scalar_type,half_scalar_type,local_ordinal_type,global_ordinal_type,node_type>(belosOpRCPHalf));
218
219 thyraPrecOp = Thyra::createLinearOp(wrappedOp);
220#else
221 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Solver does not have correct precisions enabled to use half precision.")
222#endif
223 } else {
224 // Wrap concrete preconditioner
225 RCP<BelosTpLinProb> belosLinProb = rcp(new BelosTpLinProb());
226 belosLinProb->setOperator(tpetraFwdOp);
227 RCP<TpetraLinOp> belosOpRCP = rcp(new BelosTpOp(belosLinProb, packageParamList, solverType, true));
228 thyraPrecOp = Thyra::createLinearOp(belosOpRCP);
229 }
230 defaultPrec->initializeUnspecified(thyraPrecOp);
231
232 totalTimer.stop();
233 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
234 *out << "\nTotal time in Thyra::BelosTpetraPreconditionerFactory::initializePrec(...) = " << totalTimer.totalElapsedTime() << " sec\n";
235 }
236
237 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
238 *out << "\nLeaving Thyra::BelosTpetraPreconditionerFactory::initializePrec(...) ...\n";
239 }
240}
241
242
243template <typename MatrixType>
245 PreconditionerBase<scalar_type> *prec,
246 Teuchos::RCP<const LinearOpSourceBase<scalar_type> > *fwdOp,
247 ESupportSolveUse *supportSolveUse
248 ) const
249{
250 // Check precondition
251
252 TEUCHOS_ASSERT(prec);
253
254 // Retrieve concrete preconditioner object
255
256 const Teuchos::Ptr<DefaultPreconditioner<scalar_type> > defaultPrec =
257 Teuchos::ptr(dynamic_cast<DefaultPreconditioner<scalar_type> *>(prec));
258 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(defaultPrec));
259
260 if (fwdOp) {
261 // TODO: Implement properly instead of returning default value
262 *fwdOp = Teuchos::null;
263 }
264
265 if (supportSolveUse) {
266 // TODO: Implement properly instead of returning default value
267 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
268 }
269
270 defaultPrec->uninitialize();
271}
272
273
274// Overridden from ParameterListAcceptor
275
276
277template <typename MatrixType>
279 Teuchos::RCP<Teuchos::ParameterList> const& paramList
280 )
281{
282 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
283
284 const auto validParamList = this->getValidParameters();
285 paramList->validateParametersAndSetDefaults(*validParamList, 0);
286 paramList->validateParameters(*validParamList, 0);
287
288 paramList_ = paramList;
289 Teuchos::readVerboseObjectSublist(paramList_.getRawPtr(), this);
290}
291
292
293template <typename MatrixType>
294Teuchos::RCP<Teuchos::ParameterList>
299
300
301template <typename MatrixType>
302Teuchos::RCP<Teuchos::ParameterList>
304{
305 const Teuchos::RCP<Teuchos::ParameterList> savedParamList = paramList_;
306 paramList_ = Teuchos::null;
307 return savedParamList;
308}
309
310
311template <typename MatrixType>
312Teuchos::RCP<const Teuchos::ParameterList>
317
318template <typename MatrixType>
319Teuchos::RCP<const Teuchos::ParameterList>
321{
322 static Teuchos::RCP<Teuchos::ParameterList> validParamList;
323
324 if (Teuchos::is_null(validParamList)) {
325 validParamList = Teuchos::rcp(new Teuchos::ParameterList("BelosPrecTpetra"));
326 validParamList->set(
327 "BelosPrec Solver Type", "GMRES",
328 "Name of Belos solver to be used as a preconditioner. (Use valid names for Belos::SolverFactory.)"
329 );
330 validParamList->set(
331 "half precision", false,
332 "Whether a half-of-standard-precision Belos-solver-as-preconditioner should be built."
333 );
334 validParamList->sublist(
335 "BelosPrec Solver Params", false,
336 "Belos solver settings that are passed onto the Belos solver itself."
337 );
338 Teuchos::setupVerboseObjectSublist(validParamList.getRawPtr());
339 }
340
341 return validParamList;
342}
343
344
345// Public functions overridden from Teuchos::Describable
346
347template <typename MatrixType>
349{
350 return "Thyra::BelosTpetraPreconditionerFactory";
351}
352
353
354} // namespace Thyra
355
356#endif // THYRA_BELOS_TPETRA_PRECONDITIONERFACTORY_DEF_HPP
Concrete preconditioner factory subclass based on Belos. (Yes, Belos solvers can also be used as prec...