Amesos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Thyra_AmesosLinearOpWithSolveFactory.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Amesos: Direct Sparse Solver Package
6// Copyright (2004) 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// This library is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Lesser General Public License as
13// published by the Free Software Foundation; either version 2.1 of the
14// License, or (at your option) any later version.
15//
16// This library is distributed in the hope that it will be useful, but
17// WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19// Lesser General Public License for more details.
20//
21// You should have received a copy of the GNU Lesser General Public
22// License along with this library; if not, write to the Free Software
23// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
24// USA
25// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
26//
27// ***********************************************************************
28// @HEADER
29*/
30
31#ifndef __sun
32
33#include "Thyra_AmesosLinearOpWithSolveFactory.hpp"
34
35#include "Thyra_AmesosLinearOpWithSolve.hpp"
36#include "Thyra_EpetraOperatorViewExtractorStd.hpp"
37#include "Amesos.h"
38#include "Teuchos_dyn_cast.hpp"
39#include "Teuchos_TimeMonitor.hpp"
40
41#ifdef HAVE_AMESOS_KLU
42#include "Amesos_Klu.h"
43#endif
44#ifdef HAVE_AMESOS_LAPACK
45#include "Amesos_Lapack.h"
46#endif
47#ifdef HAVE_AMESOS_MUMPS
48#include "Amesos_Mumps.h"
49#endif
50#ifdef HAVE_AMESOS_SCALAPACK
51#include "Amesos_Scalapack.h"
52#endif
53#ifdef HAVE_AMESOS_UMFPACK
54#include "Amesos_Umfpack.h"
55#endif
56#ifdef HAVE_AMESOS_SUPERLUDIST
57#include "Amesos_Superludist.h"
58#endif
59#ifdef HAVE_AMESOS_SUPERLU
60#include "Amesos_Superlu.h"
61#endif
62#ifdef HAVE_AMESOS_DSCPACK
63#include "Amesos_Dscpack.h"
64#endif
65#if defined(HAVE_AMESOS_PARDISO) || defined(HAVE_AMESOS_PARDISO_MKL)
66#include "Amesos_Pardiso.h"
67#endif
68#ifdef HAVE_AMESOS_TAUCS
69#include "Amesos_Taucs.h"
70#endif
71#ifdef HAVE_AMESOS_PARAKLETE
72#include "Amesos_Paraklete.h"
73#endif
74
75namespace {
76
77Teuchos::RCP<Teuchos::Time> overallTimer, constructTimer, symbolicTimer, factorTimer;
78
79const std::string epetraFwdOp_str = "epetraFwdOp";
80
81} // namespace
82
83namespace Thyra {
84
85
86// Parameter names for Paramter List
87
88const std::string AmesosLinearOpWithSolveFactory::SolverType_name = "Solver Type";
89
90const std::string AmesosLinearOpWithSolveFactory::RefactorizationPolicy_name = "Refactorization Policy";
91
92const std::string AmesosLinearOpWithSolveFactory::ThrowOnPreconditionerInput_name = "Throw on Preconditioner Input";
93
94const std::string AmesosLinearOpWithSolveFactory::Amesos_Settings_name = "Amesos Settings";
95
96// Constructors/initializers/accessors
97
98AmesosLinearOpWithSolveFactory::~AmesosLinearOpWithSolveFactory()
99{
100#ifdef TEUCHOS_DEBUG
101 if(paramList_.get())
102 paramList_->validateParameters(
103 *this->getValidParameters(),0 // Only validate this level for now!
104 );
105#endif
106}
107
108AmesosLinearOpWithSolveFactory::AmesosLinearOpWithSolveFactory(
109 const Amesos::ESolverType solverType
110 ,const Amesos::ERefactorizationPolicy refactorizationPolicy
111 ,const bool throwOnPrecInput
112 )
113 :epetraFwdOpViewExtractor_(Teuchos::rcp(new EpetraOperatorViewExtractorStd()))
114 ,solverType_(solverType)
115 ,refactorizationPolicy_(refactorizationPolicy)
116 ,throwOnPrecInput_(throwOnPrecInput)
117{
118 initializeTimers();
119}
120
121// Overridden from LinearOpWithSolveFactoryBase
122
123bool AmesosLinearOpWithSolveFactory::isCompatible(
124 const LinearOpSourceBase<double> &fwdOpSrc
125 ) const
126{
127 Teuchos::RCP<const LinearOpBase<double> >
128 fwdOp = fwdOpSrc.getOp();
129 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
130 ETransp epetraFwdOpTransp;
131 EApplyEpetraOpAs epetraFwdOpApplyAs;
132 EAdjointEpetraOp epetraFwdOpAdjointSupport;
133 double epetraFwdOpScalar;
134 epetraFwdOpViewExtractor_->getEpetraOpView(
135 fwdOp
136 ,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
137 );
138 if( !dynamic_cast<const Epetra_RowMatrix*>(&*epetraFwdOp) )
139 return false;
140 return true;
141}
142
143Teuchos::RCP<LinearOpWithSolveBase<double> >
144AmesosLinearOpWithSolveFactory::createOp() const
145{
146 return Teuchos::rcp(new AmesosLinearOpWithSolve());
147}
148
149void AmesosLinearOpWithSolveFactory::initializeOp(
150 const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc
151 ,LinearOpWithSolveBase<double> *Op
152 ,const ESupportSolveUse supportSolveUse
153 ) const
154{
155 Teuchos::TimeMonitor overallTimeMonitor(*overallTimer);
156#ifdef TEUCHOS_DEBUG
157 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
158#endif
159 const Teuchos::RCP<const LinearOpBase<double> >
160 fwdOp = fwdOpSrc->getOp();
161 //
162 // Unwrap and get the forward Epetra_Operator object
163 //
164 Teuchos::RCP<const Epetra_Operator> epetraFwdOp;
165 ETransp epetraFwdOpTransp;
166 EApplyEpetraOpAs epetraFwdOpApplyAs;
167 EAdjointEpetraOp epetraFwdOpAdjointSupport;
168 double epetraFwdOpScalar;
169 epetraFwdOpViewExtractor_->getEpetraOpView(
170 fwdOp,&epetraFwdOp,&epetraFwdOpTransp,&epetraFwdOpApplyAs,&epetraFwdOpAdjointSupport,&epetraFwdOpScalar
171 );
172 // Get the AmesosLinearOpWithSolve object
173 AmesosLinearOpWithSolve
174 *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
175 //
176 // Determine if we must start over or not
177 //
178 bool startOver = ( amesosOp->get_amesosSolver()==Teuchos::null );
179 if(!startOver) {
180 startOver =
181 (
182 epetraFwdOpTransp != amesosOp->get_amesosSolverTransp() ||
183 epetraFwdOp.get() != amesosOp->get_epetraLP()->GetOperator()
184 // We must start over if the matrix object changes. This is a
185 // weakness of Amesos but there is nothing I can do about this right
186 // now!
187 );
188 }
189 //
190 // Update the amesos solver
191 //
192 if(startOver) {
193 //
194 // This LOWS object has not be initialized yet or is not compatible with the existing
195 //
196 // so this is where we setup everything from the ground up.
197 //
198 // Create the linear problem and set the operator with memory of RCP to Epetra_Operator view!
199 Teuchos::RCP<Epetra_LinearProblem>
200 epetraLP = Teuchos::rcp(new Epetra_LinearProblem());
201 epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
202 Teuchos::set_extra_data< Teuchos::RCP<const Epetra_Operator> >( epetraFwdOp,
203 epetraFwdOp_str, Teuchos::inOutArg(epetraLP) );
204 // Create the concrete solver
205 Teuchos::RCP<Amesos_BaseSolver>
206 amesosSolver;
207 {
208 Teuchos::TimeMonitor constructTimeMonitor(*constructTimer);
209 switch(solverType_) {
210 case Thyra::Amesos::LAPACK :
211 amesosSolver = Teuchos::rcp(new Amesos_Lapack(*epetraLP));
212 break;
213#ifdef HAVE_AMESOS_KLU
214 case Thyra::Amesos::KLU :
215 amesosSolver = Teuchos::rcp(new Amesos_Klu(*epetraLP));
216 break;
217#endif
218#ifdef HAVE_AMESOS_MUMPS
219 case Thyra::Amesos::MUMPS :
220 amesosSolver = Teuchos::rcp(new Amesos_Mumps(*epetraLP));
221 break;
222#endif
223#ifdef HAVE_AMESOS_SCALAPACK
224 case Thyra::Amesos::SCALAPACK :
225 amesosSolver = Teuchos::rcp(new Amesos_Scalapack(*epetraLP));
226 break;
227#endif
228#ifdef HAVE_AMESOS_UMFPACK
229 case Thyra::Amesos::UMFPACK :
230 amesosSolver = Teuchos::rcp(new Amesos_Umfpack(*epetraLP));
231 break;
232#endif
233#ifdef HAVE_AMESOS_SUPERLUDIST
234 case Thyra::Amesos::SUPERLUDIST :
235 amesosSolver = Teuchos::rcp(new Amesos_Superludist(*epetraLP));
236 break;
237#endif
238#ifdef HAVE_AMESOS_SUPERLU
239 case Thyra::Amesos::SUPERLU :
240 amesosSolver = Teuchos::rcp(new Amesos_Superlu(*epetraLP));
241 break;
242#endif
243#ifdef HAVE_AMESOS_DSCPACK
244 case Thyra::Amesos::DSCPACK :
245 amesosSolver = Teuchos::rcp(new Amesos_Dscpack(*epetraLP));
246 break;
247#endif
248#ifdef HAVE_AMESOS_PARDISO
249 case Thyra::Amesos::PARDISO :
250 amesosSolver = Teuchos::rcp(new Amesos_Pardiso(*epetraLP));
251 break;
252#endif
253#ifdef HAVE_AMESOS_TAUCS
254 case Thyra::Amesos::TAUCS :
255 amesosSolver = Teuchos::rcp(new Amesos_Taucs(*epetraLP));
256 break;
257#endif
258#ifdef HAVE_AMESOS_PARAKLETE
259 case Thyra::Amesos::PARAKLETE :
260 amesosSolver = Teuchos::rcp(new Amesos_Paraklete(*epetraLP));
261 break;
262#endif
263 default:
264 TEUCHOS_TEST_FOR_EXCEPTION(
265 true, std::logic_error
266 ,"Error, the solver type ID = " << solverType_ << " is invalid!"
267 );
268 }
269 }
270 // Set the parameters
271 if(paramList_.get()) amesosSolver->setParameterList(sublist(paramList_,"Amesos Settings"));
272 // Do the initial factorization
273 {
274 Teuchos::TimeMonitor symbolicTimeMonitor(*symbolicTimer);
275 amesosSolver->SymbolicFactorization();
276 }
277 {
278 Teuchos::TimeMonitor factorTimeMonitor(*factorTimer);
279 amesosSolver->NumericFactorization();
280 }
281 // Initialize the LOWS object and we are done!
282 amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
283 }
284 else {
285 //
286 // This LOWS object has already be initialized once so we must just reset
287 // the matrix and refactor it.
288 //
289 // Get non-const pointers to the linear problem and the amesos solver.
290 // These const-casts are just fine since the amesosOp in non-const.
291 Teuchos::RCP<Epetra_LinearProblem>
292 epetraLP = Teuchos::rcp_const_cast<Epetra_LinearProblem>(amesosOp->get_epetraLP());
293 Teuchos::RCP<Amesos_BaseSolver>
294 amesosSolver = amesosOp->get_amesosSolver();
295 // Reset the forward operator with memory of RCP to Epetra_Operator view!
296 epetraLP->SetOperator(const_cast<Epetra_Operator*>(&*epetraFwdOp));
297 Teuchos::get_extra_data< Teuchos::RCP<const Epetra_Operator> >(epetraLP,epetraFwdOp_str) = epetraFwdOp;
298 // Reset the parameters
299 if(paramList_.get()) amesosSolver->setParameterList(sublist(paramList_,Amesos_Settings_name));
300 // Repivot if asked
301 if(refactorizationPolicy_==Amesos::REPIVOT_ON_REFACTORIZATION) {
302 Teuchos::TimeMonitor symbolicTimeMonitor(*symbolicTimer);
303 amesosSolver->SymbolicFactorization();
304 }
305 {
306 Teuchos::TimeMonitor factorTimeMonitor(*factorTimer);
307 amesosSolver->NumericFactorization();
308 }
309 // Reinitialize the LOWS object and we are done! (we must do this to get the
310 // possibly new transpose and scaling factors back in)
311 amesosOp->initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,epetraFwdOpTransp,epetraFwdOpScalar);
312 }
313 amesosOp->setOStream(this->getOStream());
314 amesosOp->setVerbLevel(this->getVerbLevel());
315}
316
317bool AmesosLinearOpWithSolveFactory::supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
318{
319 return false;
320}
321
322void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
323 const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc
324 ,const Teuchos::RCP<const PreconditionerBase<double> > &prec
325 ,LinearOpWithSolveBase<double> *Op
326 ,const ESupportSolveUse supportSolveUse
327 ) const
328{
329 TEUCHOS_TEST_FOR_EXCEPTION(
330 this->throwOnPrecInput_, std::logic_error
331 ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support preconditioners "
332 "and has been configured to throw this exception when the initializePreconditionedOp(...) function is called!"
333 );
334 this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the preconditioner!
335}
336
337void AmesosLinearOpWithSolveFactory::initializePreconditionedOp(
338 const Teuchos::RCP<const LinearOpSourceBase<double> > &fwdOpSrc
339 ,const Teuchos::RCP<const LinearOpSourceBase<double> > &approxFwdOpSrc
340 ,LinearOpWithSolveBase<double> *Op
341 ,const ESupportSolveUse supportSolveUse
342 ) const
343{
344 TEUCHOS_TEST_FOR_EXCEPTION(
345 this->throwOnPrecInput_, std::logic_error
346 ,"Error, the concrete implementation described as \'"<<this->description()<<"\' does not support preconditioners "
347 "and has been configured to throw this exception when the initializePreconditionedOp(...) function is called!"
348 );
349 this->initializeOp(fwdOpSrc,Op,supportSolveUse); // Ignore the preconditioner!
350}
351
352void AmesosLinearOpWithSolveFactory::uninitializeOp(
353 LinearOpWithSolveBase<double> *Op
354 ,Teuchos::RCP<const LinearOpSourceBase<double> > *fwdOpSrc
355 ,Teuchos::RCP<const PreconditionerBase<double> > *prec
356 ,Teuchos::RCP<const LinearOpSourceBase<double> > *approxFwdOpSrc
357 ,ESupportSolveUse *supportSolveUse
358 ) const
359{
360#ifdef TEUCHOS_DEBUG
361 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
362#endif
363 AmesosLinearOpWithSolve
364 *amesosOp = &Teuchos::dyn_cast<AmesosLinearOpWithSolve>(*Op);
365 Teuchos::RCP<const LinearOpSourceBase<double> >
366 _fwdOpSrc = amesosOp->extract_fwdOpSrc(); // Will be null if uninitialized!
367 if(_fwdOpSrc.get()) {
368 // Erase the Epetra_Operator view of the forward operator!
369 Teuchos::RCP<Epetra_LinearProblem> epetraLP = amesosOp->get_epetraLP();
370 Teuchos::get_extra_data< Teuchos::RCP<const Epetra_Operator> >(
371 epetraLP,epetraFwdOp_str
372 )
373 = Teuchos::null;
374 // Note, we did not erase the address of the operator in
375 // epetraLP->GetOperator() since it seems that the amesos solvers do not
376 // recheck the value of GetProblem()->GetOperator() so you had better not
377 // rest this!
378 }
379 if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc; // It is fine if the client does not want this object back!
380 if(prec) *prec = Teuchos::null; // We never keep a preconditioner!
381 if(approxFwdOpSrc) *approxFwdOpSrc = Teuchos::null; // We never keep an approximate fwd operator!
382}
383
384// Overridden from ParameterListAcceptor
385
386void AmesosLinearOpWithSolveFactory::setParameterList(
387 Teuchos::RCP<Teuchos::ParameterList> const& paramList
388 )
389{
390 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
391 paramList->validateParameters(*this->getValidParameters(),0); // Only validate this level for now!
392 paramList_ = paramList;
393 solverType_ =
394 Amesos::solverTypeNameToEnumMap.get<Amesos::ESolverType>(
395 paramList_->get(
396 SolverType_name
397 ,Amesos::toString(solverType_)
398 )
399 ,paramList_->name()+"->"+SolverType_name
400 );
401 refactorizationPolicy_ =
402 Amesos::refactorizationPolicyNameToEnumMap.get<Amesos::ERefactorizationPolicy>(
403 paramList_->get(
404 RefactorizationPolicy_name
405 ,Amesos::toString(refactorizationPolicy_)
406 )
407 ,paramList_->name()+"->"+RefactorizationPolicy_name
408 );
409 throwOnPrecInput_ = paramList_->get(ThrowOnPreconditionerInput_name,throwOnPrecInput_);
410}
411
412Teuchos::RCP<Teuchos::ParameterList>
413AmesosLinearOpWithSolveFactory::getNonconstParameterList()
414{
415 return paramList_;
416}
417
418Teuchos::RCP<Teuchos::ParameterList>
419AmesosLinearOpWithSolveFactory::unsetParameterList()
420{
421 Teuchos::RCP<Teuchos::ParameterList> _paramList = paramList_;
422 paramList_ = Teuchos::null;
423 return _paramList;
424}
425
426Teuchos::RCP<const Teuchos::ParameterList>
427AmesosLinearOpWithSolveFactory::getParameterList() const
428{
429 return paramList_;
430}
431
432Teuchos::RCP<const Teuchos::ParameterList>
433AmesosLinearOpWithSolveFactory::getValidParameters() const
434{
435 return generateAndGetValidParameters();
436}
437
438// Public functions overridden from Teuchos::Describable
439
440std::string AmesosLinearOpWithSolveFactory::description() const
441{
442 std::ostringstream oss;
443 oss << "Thyra::AmesosLinearOpWithSolveFactory{";
444 oss << "solverType=" << toString(solverType_);
445 oss << "}";
446 return oss.str();
447}
448
449// private
450
451
452void AmesosLinearOpWithSolveFactory::initializeTimers()
453{
454 if(!overallTimer.get()) {
455 overallTimer = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF");
456 constructTimer = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:InitConstruct");
457 symbolicTimer = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:Symbolic");
458 factorTimer = Teuchos::TimeMonitor::getNewTimer("AmesosLOWSF:Factor");
459 }
460}
461
462Teuchos::RCP<const Teuchos::ParameterList>
463AmesosLinearOpWithSolveFactory::generateAndGetValidParameters()
464{
465 static Teuchos::RCP<Teuchos::ParameterList> validParamList;
466 if(validParamList.get()==NULL) {
467 validParamList = Teuchos::rcp(new Teuchos::ParameterList("Amesos"));
468 validParamList->set(
469 SolverType_name
470#ifdef HAVE_AMESOS_KLU
471 ,Amesos::toString(Amesos::KLU)
472#else
473 ,Amesos::toString(Amesos::LAPACK)
474#endif
475 );
476 validParamList->set(RefactorizationPolicy_name,Amesos::toString(Amesos::REPIVOT_ON_REFACTORIZATION));
477 validParamList->set(ThrowOnPreconditionerInput_name,bool(true));
478 validParamList->sublist(Amesos_Settings_name).setParameters(::Amesos::GetValidParameters());
479 }
480 return validParamList;
481}
482
483} // namespace Thyra
484
485#endif // __sun
Amesos_Dscpack: An object-oriented wrapper for Dscpack.
Amesos_Klu: A serial, unblocked code ideal for getting started and for very sparse matrices,...
Definition Amesos_Klu.h:116
Amesos_Lapack: an interface to LAPACK.
Amesos_Mumps: An object-oriented wrapper for the double precision version of MUMPS.
Amesos_Paraklete: A serial, unblocked code ideal for getting started and for very sparse matrices,...
Amesos_Pardiso: Interface to the PARDISO package.
Amesos_Scalapack: A serial and parallel dense solver. For now, we implement only the unsymmetric ScaL...
Amesos_Superlu: Amesos interface to Xioye Li's SuperLU 3.0 serial code.
Amesos_Superludist: An object-oriented wrapper for Superludist.
Amesos_Taucs: An interface to the TAUCS package.
Class Amesos_Umfpack: An object-oriented wrapper for UMFPACK.
static Teuchos::ParameterList GetValidParameters()
Get the list of valid parameters.
Definition Amesos.cpp:303