42#ifndef THYRA_BELOS_TPETRA_PRECONDITIONERFACTORY_DEF_HPP
43#define THYRA_BELOS_TPETRA_PRECONDITIONERFACTORY_DEF_HPP
47#include "Thyra_DefaultPreconditioner.hpp"
48#include "Thyra_TpetraLinearOp.hpp"
49#include "Thyra_TpetraThyraWrappers.hpp"
51#include "BelosTpetraOperator.hpp"
52#include "BelosConfigDefs.hpp"
53#include "BelosLinearProblem.hpp"
54#include "BelosTpetraAdapter.hpp"
56#include "Tpetra_MixedScalarMultiplyOp.hpp"
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"
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
82template <
typename MatrixType>
90template <
typename MatrixType>
92 const LinearOpSourceBase<scalar_type> &fwdOpSrc
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;
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);
103 const Teuchos::RCP<const MatrixType> tpetraFwdMatrix = Teuchos::rcp_dynamic_cast<const MatrixType>(tpetraFwdOp,
true);
105 return Teuchos::nonnull(tpetraFwdMatrix);
109template <
typename MatrixType>
110Teuchos::RCP<PreconditionerBase<typename BelosTpetraPreconditionerFactory<MatrixType>::scalar_type> >
113 return Teuchos::rcp(
new DefaultPreconditioner<scalar_type>);
117template <
typename MatrixType>
119 const Teuchos::RCP<
const LinearOpSourceBase<scalar_type> > &fwdOpSrc,
120 PreconditionerBase<scalar_type> *prec,
121 const ESupportSolveUse
129 TEUCHOS_ASSERT(Teuchos::nonnull(fwdOpSrc));
130 TEUCHOS_ASSERT(this->isCompatible(*fwdOpSrc));
131 TEUCHOS_ASSERT(prec);
133 Teuchos::Time totalTimer(
""), timer(
"");
134 totalTimer.start(
true);
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";
145 const RCP<const LinearOpBase<scalar_type> > fwdOp = fwdOpSrc->getOp();
146 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(fwdOp));
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;
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));
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;
162#ifdef THYRA_BELOS_PREC_ENABLE_HALF_PRECISION
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;
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));
180 RCP<Teuchos::ParameterList> innerParamList;
181 if (paramList_.is_null ()) {
182 innerParamList = rcp(
new Teuchos::ParameterList(*getValidParameters()));
185 innerParamList = paramList_;
188 bool useHalfPrecision =
false;
189 if (innerParamList->isParameter(
"half precision"))
190 useHalfPrecision = Teuchos::getParameter<bool>(*innerParamList,
"half precision");
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"));
196 std::string solverTypeUpper (solverType);
197 std::transform(solverTypeUpper.begin(), solverTypeUpper.end(),solverTypeUpper.begin(), ::toupper);
200 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
201 *out <<
"\nCreating a new BelosTpetra::Preconditioner object...\n";
203 RCP<LinearOpBase<scalar_type> > thyraPrecOp;
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";
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>();
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));
219 thyraPrecOp = Thyra::createLinearOp(wrappedOp);
221 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Solver does not have correct precisions enabled to use half precision.")
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);
230 defaultPrec->initializeUnspecified(thyraPrecOp);
233 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_LOW)) {
234 *out <<
"\nTotal time in Thyra::BelosTpetraPreconditionerFactory::initializePrec(...) = " << totalTimer.totalElapsedTime() <<
" sec\n";
237 if (Teuchos::nonnull(out) && Teuchos::includesVerbLevel(verbLevel, Teuchos::VERB_MEDIUM)) {
238 *out <<
"\nLeaving Thyra::BelosTpetraPreconditionerFactory::initializePrec(...) ...\n";
243template <
typename MatrixType>
245 PreconditionerBase<scalar_type> *prec,
246 Teuchos::RCP<
const LinearOpSourceBase<scalar_type> > *fwdOp,
247 ESupportSolveUse *supportSolveUse
252 TEUCHOS_ASSERT(prec);
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));
262 *fwdOp = Teuchos::null;
265 if (supportSolveUse) {
267 *supportSolveUse = Thyra::SUPPORT_SOLVE_UNSPECIFIED;
270 defaultPrec->uninitialize();
277template <
typename MatrixType>
279 Teuchos::RCP<Teuchos::ParameterList>
const& paramList
282 TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(paramList));
284 const auto validParamList = this->getValidParameters();
285 paramList->validateParametersAndSetDefaults(*validParamList, 0);
286 paramList->validateParameters(*validParamList, 0);
288 paramList_ = paramList;
289 Teuchos::readVerboseObjectSublist(paramList_.getRawPtr(),
this);
293template <
typename MatrixType>
294Teuchos::RCP<Teuchos::ParameterList>
301template <
typename MatrixType>
302Teuchos::RCP<Teuchos::ParameterList>
305 const Teuchos::RCP<Teuchos::ParameterList> savedParamList = paramList_;
306 paramList_ = Teuchos::null;
307 return savedParamList;
311template <
typename MatrixType>
312Teuchos::RCP<const Teuchos::ParameterList>
318template <
typename MatrixType>
319Teuchos::RCP<const Teuchos::ParameterList>
322 static Teuchos::RCP<Teuchos::ParameterList> validParamList;
324 if (Teuchos::is_null(validParamList)) {
325 validParamList = Teuchos::rcp(
new Teuchos::ParameterList(
"BelosPrecTpetra"));
327 "BelosPrec Solver Type",
"GMRES",
328 "Name of Belos solver to be used as a preconditioner. (Use valid names for Belos::SolverFactory.)"
331 "half precision",
false,
332 "Whether a half-of-standard-precision Belos-solver-as-preconditioner should be built."
334 validParamList->sublist(
335 "BelosPrec Solver Params",
false,
336 "Belos solver settings that are passed onto the Belos solver itself."
338 Teuchos::setupVerboseObjectSublist(validParamList.getRawPtr());
341 return validParamList;
347template <
typename MatrixType>
350 return "Thyra::BelosTpetraPreconditionerFactory";
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > ¶mList)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
bool isCompatible(const LinearOpSourceBase< scalar_type > &fwdOp) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void uninitializePrec(PreconditionerBase< scalar_type > *prec, Teuchos::RCP< const LinearOpSourceBase< scalar_type > > *fwdOp, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
std::string description() const
void initializePrec(const Teuchos::RCP< const LinearOpSourceBase< scalar_type > > &fwdOp, PreconditionerBase< scalar_type > *prec, const ESupportSolveUse supportSolveUse) const
BelosTpetraPreconditionerFactory()
Teuchos::RCP< PreconditionerBase< scalar_type > > createPrec() const