Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Thyra_Amesos2LinearOpWithSolveFactory_def.hpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stratimikos: Thyra-based strategies for linear solvers
6// Copyright (2006) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
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 Roscoe A. Bartlett (rabartl@sandia.gov)
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44#ifndef THYRA_AMESOS2_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
45#define THYRA_AMESOS2_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
46
48
49#include "Thyra_Amesos2LinearOpWithSolve.hpp"
50#include "Amesos2.hpp"
51#include "Amesos2_Details_LinearSolverFactory.hpp"
52#include "Amesos2_Version.hpp"
53#include "Amesos2_Factory.hpp"
54#include "Thyra_TpetraLinearOp.hpp"
55#include "Thyra_TpetraThyraWrappers.hpp"
56#include "Thyra_DefaultDiagonalLinearOp.hpp"
57#include "Teuchos_dyn_cast.hpp"
58#include "Teuchos_TimeMonitor.hpp"
59#include "Teuchos_TypeNameTraits.hpp"
60#include "Teuchos_TypeTraits.hpp"
61#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
62
63namespace Thyra {
64
65
66// Parameter names for Paramter List
67
68template<typename Scalar>
70 = "Solver Type";
71
72template<typename Scalar>
74 = "Refactorization Policy";
75
76template<typename Scalar>
78 = "Throw on Preconditioner Input";
79
80template<typename Scalar>
82 = "Amesos2 Settings";
83
84// Constructors/initializers/accessors
85
86template<typename Scalar>
88{
89#ifdef TEUCHOS_DEBUG
90 if(paramList_.get())
91 paramList_->validateParameters(
92 *this->getValidParameters(),0 // Only validate this level for now!
93 );
94#endif
95}
96
97template<typename Scalar>
99 const Amesos2::ESolverType solverType,
100 const Amesos2::ERefactorizationPolicy refactorizationPolicy,
101 const bool throwOnPrecInput
102 )
103 :solverType_(solverType)
104 ,refactorizationPolicy_(refactorizationPolicy)
105 ,throwOnPrecInput_(throwOnPrecInput)
106{
107}
108
109// Overridden from LinearOpWithSolveFactoryBase
110
111template<typename Scalar>
113 const LinearOpSourceBase<Scalar> &fwdOpSrc
114 ) const
115{
116 Teuchos::RCP<const LinearOpBase<Scalar> >
117 fwdOp = fwdOpSrc.getOp();
118 auto tpetraFwdOp = ConverterT::getConstTpetraOperator(fwdOp);
119 if ( ! dynamic_cast<const MAT * >(&*tpetraFwdOp) )
120 return false;
121 return true;
122}
123
124template<typename Scalar>
125RCP<LinearOpWithSolveBase<Scalar> >
127{
128 return Teuchos::rcp(new Amesos2LinearOpWithSolve<Scalar>());
129}
130
131template<typename Scalar>
133 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
134 LinearOpWithSolveBase<Scalar> *Op,
135 const ESupportSolveUse /* supportSolveUse */
136 ) const
137{
138 THYRA_FUNC_TIME_MONITOR("Stratimikos: Amesos2LOWSF");
139
140 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
141 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
142 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
143 RCP<const LinearOpBase<Scalar> > fwdOp = fwdOpSrc->getOp();
144
145 Teuchos::RCP<Teuchos::FancyOStream>
146 out = Teuchos::VerboseObjectBase::getDefaultOStream();
147
148 //
149 // Unwrap and get the forward Tpetra::Operator object
150 //
151 auto tpetraFwdOp = ConverterT::getConstTpetraOperator(fwdOp);
152 auto tpetraCrsMat = Teuchos::rcp_dynamic_cast<const MAT>(tpetraFwdOp);
153 // Get the Amesos2LinearOpWithSolve object
154 Amesos2LinearOpWithSolve<Scalar>
155 *amesos2Op = &Teuchos::dyn_cast<Amesos2LinearOpWithSolve<Scalar>>(*Op);
156
157 //
158 // Determine if we must start over or not
159 //
160 bool startOver = ( amesos2Op->get_amesos2Solver()==Teuchos::null );
161 if (!startOver) {
162 auto oldTpetraFwdOp = ConverterT::getConstTpetraOperator(amesos2Op->get_fwdOp());
163 startOver =
164 (
165 tpetraFwdOp.get() != oldTpetraFwdOp.get()
166 // Assuming that, like Amesos, Amesos2 must start over if the matrix changes
167 );
168 }
169 //
170 // Update the amesos2 solver
171 //
172 if (startOver) {
173 //
174 // This LOWS object has not be initialized yet or is not compatible with the existing
175 //
176 // so this is where we setup everything from the ground up.
177
178 // Create the concrete solver
179 Teuchos::RCP<Solver> amesos2Solver;
180 {
181 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: Amesos2LOWSF:InitConstruct",
182 InitConstruct);
183 switch(solverType_) {
185 amesos2Solver = ::Amesos2::create<MAT,MV>("klu2", tpetraCrsMat);
186 break;
187#ifdef HAVE_AMESOS2_LAPACK
188 case Thyra::Amesos2::LAPACK:
189 amesos2Solver = ::Amesos2::create<MAT,MV>("lapack", tpetraCrsMat);
190 break;
191#endif
192#ifdef HAVE_AMESOS2_SUPERLU
193 case Thyra::Amesos2::SUPERLU:
194 amesos2Solver = ::Amesos2::create<MAT,MV>("superlu", tpetraCrsMat);
195 break;
196#endif
197#ifdef HAVE_AMESOS2_SUPERLUMT
198 case Thyra::Amesos2::SUPERLUMT:
199 amesos2Solver = ::Amesos2::create<MAT,MV>("superlumt", tpetraCrsMat);
200 break;
201#endif
202#ifdef HAVE_AMESOS2_SUPERLUDIST
203 case Thyra::Amesos2::SUPERLUDIST:
204 amesos2Solver = ::Amesos2::create<MAT,MV>("superludist", tpetraCrsMat);
205 break;
206# endif
207#ifdef HAVE_AMESOS2_PARDISO_MKL
208 case Thyra::Amesos2::PARDISO_MKL:
209 amesos2Solver = ::Amesos2::create<MAT,MV>("pardiso_mkl", tpetraCrsMat);
210 break;
211#endif
212#ifdef HAVE_AMESOS2_CHOLMOD
213 case Thyra::Amesos2::CHOLMOD:
214 amesos2Solver = ::Amesos2::create<MAT,MV>("cholmod", tpetraCrsMat);
215 break;
216#endif
217#ifdef HAVE_AMESOS2_BASKER
218 case Thyra::Amesos2::BASKER:
219 amesos2Solver = ::Amesos2::create<MAT,MV>("basker", tpetraCrsMat);
220 break;
221#endif
222#ifdef HAVE_AMESOS2_MUMPS
223 case Thyra::Amesos2::MUMPS:
224 amesos2Solver = ::Amesos2::create<MAT,MV>("mumps", tpetraCrsMat);
225 break;
226#endif
227 default:
228 TEUCHOS_TEST_FOR_EXCEPTION(
229 true, std::logic_error
230 ,"Error, the solver type ID = " << solverType_ << " is invalid!"
231 );
232 }
233 }
234
235 // Do the initial factorization
236 {
237 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: Amesos2LOWSF:Symbolic", Symbolic);
238 amesos2Solver->symbolicFactorization();
239 }
240 {
241 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: Amesos2LOWSF:Factor", Factor);
242 amesos2Solver->numericFactorization();
243 }
244
245 // filter out the Stratimikos adapter parameters and hand
246 // parameters down into the Solver
247 const Teuchos::RCP<Teuchos::ParameterList> dup_list
248 = Teuchos::rcp(new Teuchos::ParameterList(*paramList_));
249 dup_list->remove(SolverType_name);
250 dup_list->remove(RefactorizationPolicy_name);
251 dup_list->remove(ThrowOnPreconditionerInput_name);
252 dup_list->remove("VerboseObject");
253 amesos2Solver->setParameters(dup_list);
254
255 // Initialize the LOWS object and we are done!
256 amesos2Op->initialize(fwdOp,fwdOpSrc,amesos2Solver);
257 }
258 else {
259 //
260 // This LOWS object has already be initialized once so we must just reset
261 // the matrix and refactor it.
262 auto amesos2Solver = amesos2Op->get_amesos2Solver();
263
264 // set
265 amesos2Solver->setA(tpetraCrsMat);
266
267 // Do the initial factorization
268 if(refactorizationPolicy_ == Amesos2::REPIVOT_ON_REFACTORIZATION) {
269 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: Amesos2LOWSF:Symbolic", Symbolic);
270 amesos2Solver->symbolicFactorization();
271 }
272 {
273 THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: Amesos2LOWSF::Factor", Factor);
274 amesos2Solver->numericFactorization();
275 }
276
277 // Initialize the LOWS object and we are done!
278 amesos2Op->initialize(fwdOp,fwdOpSrc,amesos2Solver);
279 }
280 amesos2Op->setOStream(this->getOStream());
281 amesos2Op->setVerbLevel(this->getVerbLevel());
282}
283
284template<typename Scalar>
285bool Amesos2LinearOpWithSolveFactory<Scalar>::supportsPreconditionerInputType(const EPreconditionerInputType /* precOpType */) const
286{
287 return false;
288}
289
290template<typename Scalar>
292 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
293 const RCP<const PreconditionerBase<Scalar> > &/* prec */,
294 LinearOpWithSolveBase<Scalar> *Op,
295 const ESupportSolveUse supportSolveUse
296 ) const
297{
298 TEUCHOS_TEST_FOR_EXCEPTION(
299 this->throwOnPrecInput_, std::logic_error,
300 "Error, the concrete implementation described as \'"<<this->description()
301 <<"\' does not support preconditioners"
302 " and has been configured to throw this exception when the"
303 " initializePreconditionedOp(...) function is called!"
304 );
305 this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the preconditioner!
306}
307
308template<typename Scalar>
309void Amesos2LinearOpWithSolveFactory<Scalar>::initializePreconditionedOp(
310 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
311 const RCP<const LinearOpSourceBase<Scalar> > &/* approxFwdOpSrc */,
312 LinearOpWithSolveBase<Scalar> *Op,
313 const ESupportSolveUse supportSolveUse
314 ) const
315{
316 TEUCHOS_TEST_FOR_EXCEPTION(
317 this->throwOnPrecInput_, std::logic_error,
318 "Error, the concrete implementation described as \'"<<this->description()
319 <<"\' does not support preconditioners"
320 " and has been configured to throw this exception when the"
321 " initializePreconditionedOp(...) function is called!"
322 );
323 this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the preconditioner!
324}
325
326template<typename Scalar>
328 LinearOpWithSolveBase<Scalar> *Op,
329 RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
330 RCP<const PreconditionerBase<Scalar> > *prec,
331 RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
332 ESupportSolveUse * /* supportSolveUse */
333 ) const
334{
335#ifdef TEUCHOS_DEBUG
336 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
337#endif
338 Amesos2LinearOpWithSolve<Scalar>
339 *amesos2Op = &Teuchos::dyn_cast<Amesos2LinearOpWithSolve<Scalar>>(*Op);
340 RCP<const LinearOpSourceBase<Scalar> >
341 _fwdOpSrc = amesos2Op->extract_fwdOpSrc(); // Will be null if uninitialized!
342 if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc; // It is fine if the client does not want this object back!
343 if(prec) *prec = Teuchos::null; // We never keep a preconditioner!
344 if(approxFwdOpSrc) *approxFwdOpSrc = Teuchos::null; // never keep approx fwd op!
345}
346
347// Overridden from ParameterListAcceptor
348
349template<typename Scalar>
351 RCP<Teuchos::ParameterList> const& paramList
352 )
353{
354 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
355 // Only validate this level for now here (expect Amesos2 to do its own
356 // validation?)
357 paramList->validateParameters(*this->getValidParameters(),0);
358 paramList_ = paramList;
359 solverType_ =
360 Amesos2::solverTypeNameToEnumMap.get<Amesos2::ESolverType>(
361 paramList_->get(
362 SolverType_name
363 ,Amesos2::toString(solverType_)
364 )
365 ,paramList_->name()+"->"+SolverType_name
366 );
367 refactorizationPolicy_ =
368 Amesos2::refactorizationPolicyNameToEnumMap.get<Amesos2::ERefactorizationPolicy>(
369 paramList_->get(
370 RefactorizationPolicy_name
371 ,Amesos2::toString(refactorizationPolicy_)
372 )
373 ,paramList_->name()+"->"+RefactorizationPolicy_name
374 );
375 throwOnPrecInput_ = paramList_->get(ThrowOnPreconditionerInput_name,throwOnPrecInput_);
376 Teuchos::readVerboseObjectSublist(&*paramList_,this);
377}
378
379template<typename Scalar>
380RCP<Teuchos::ParameterList>
385
386template<typename Scalar>
387RCP<Teuchos::ParameterList>
389{
390 RCP<Teuchos::ParameterList> _paramList = paramList_;
391 paramList_ = Teuchos::null;
392 return _paramList;
393}
394
395template<typename Scalar>
396RCP<const Teuchos::ParameterList>
398{
399 return paramList_;
400}
401
402template<typename Scalar>
403RCP<const Teuchos::ParameterList>
405{
406 return generateAndGetValidParameters();
407}
408
409// Public functions overridden from Teuchos::Describable
410
411template<typename Scalar>
413{
414 std::ostringstream oss;
415 oss << "Thyra::Amesos2LinearOpWithSolveFactory{";
416 oss << "solverType=" << toString(solverType_);
417 oss << "}";
418 return oss.str();
419}
420
421// private
422
423template<typename Scalar>
424RCP<const Teuchos::ParameterList>
426{
427 static RCP<Teuchos::ParameterList> validParamList;
428 if (validParamList.get()==NULL) {
429 validParamList = Teuchos::rcp(new Teuchos::ParameterList("Amesos2"));
430 validParamList->set(SolverType_name, Thyra::Amesos2::solverTypeNames[0]);
431 validParamList->set(RefactorizationPolicy_name,
432 Amesos2::toString(Amesos2::REPIVOT_ON_REFACTORIZATION));
433 validParamList->set(ThrowOnPreconditionerInput_name,bool(true));
434 Teuchos::setupVerboseObjectSublist(&*validParamList);
435 }
436 return validParamList;
437}
438
439} // namespace Thyra
440
441#endif // THYRA_AMESOS2_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
typename Amesos2LinearOpWithSolve< Scalar >::MAT MAT
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
ERefactorizationPolicy
The policy used on refactoring a matrix.
const char * solverTypeNames[numSolverTypes]