44#include "Thyra_AmesosLinearOpWithSolve.hpp"
45#include "Thyra_EpetraThyraWrappers.hpp"
46#include "Thyra_MultiVectorStdOps.hpp"
47#include "Epetra_MultiVector.h"
48#include "Teuchos_TimeMonitor.hpp"
58 amesosSolverTransp_(Thyra::
NOTRANS),
59 amesosSolverScalar_(1.0)
64 const Teuchos::RCP<
const LinearOpBase<double> > &fwdOp,
65 const Teuchos::RCP<
const LinearOpSourceBase<double> > &fwdOpSrc,
66 const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
67 const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
68 const EOpTransp amesosSolverTransp,
69 const double amesosSolverScalar
72 this->
initialize(fwdOp,fwdOpSrc,epetraLP,amesosSolver,
73 amesosSolverTransp,amesosSolverScalar);
78 const Teuchos::RCP<
const LinearOpBase<double> > &fwdOp,
79 const Teuchos::RCP<
const LinearOpSourceBase<double> > &fwdOpSrc,
80 const Teuchos::RCP<Epetra_LinearProblem> &epetraLP,
81 const Teuchos::RCP<Amesos_BaseSolver> &amesosSolver,
82 const EOpTransp amesosSolverTransp,
83 const double amesosSolverScalar
87 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
88 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
89 TEUCHOS_TEST_FOR_EXCEPT(epetraLP.get()==NULL);
90 TEUCHOS_TEST_FOR_EXCEPT(amesosSolver.get()==NULL);
91 TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetLHS()!=NULL);
92 TEUCHOS_TEST_FOR_EXCEPT(epetraLP->GetRHS()!=NULL);
97 amesosSolver_ = amesosSolver;
98 amesosSolverTransp_ = amesosSolverTransp;
99 amesosSolverScalar_ = amesosSolverScalar;
100 const std::string fwdOpLabel = fwdOp_->getObjectLabel();
101 if(fwdOpLabel.length())
102 this->setObjectLabel(
"lows("+fwdOpLabel+
")" );
106Teuchos::RCP<const LinearOpSourceBase<double> >
109 Teuchos::RCP<const LinearOpSourceBase<double> >
110 _fwdOpSrc = fwdOpSrc_;
111 fwdOpSrc_ = Teuchos::null;
117 Teuchos::RCP<
const LinearOpBase<double> > *fwdOp,
118 Teuchos::RCP<
const LinearOpSourceBase<double> > *fwdOpSrc,
119 Teuchos::RCP<Epetra_LinearProblem> *epetraLP,
120 Teuchos::RCP<Amesos_BaseSolver> *amesosSolver,
121 EOpTransp *amesosSolverTransp,
122 double *amesosSolverScalar
126 if(fwdOp) *fwdOp = fwdOp_;
127 if(fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
128 if(epetraLP) *epetraLP = epetraLP_;
129 if(amesosSolver) *amesosSolver = amesosSolver_;
130 if(amesosSolverTransp) *amesosSolverTransp = amesosSolverTransp_;
131 if(amesosSolverScalar) *amesosSolverScalar = amesosSolverScalar_;
133 fwdOp_ = Teuchos::null;
134 fwdOpSrc_ = Teuchos::null;
135 epetraLP_ = Teuchos::null;
136 amesosSolver_ = Teuchos::null;
138 amesosSolverScalar_ = 0.0;
146Teuchos::RCP< const VectorSpaceBase<double> >
149 return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
153Teuchos::RCP< const VectorSpaceBase<double> >
156 return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
160Teuchos::RCP<const LinearOpBase<double> >
163 return Teuchos::null;
172 std::ostringstream oss;
173 oss << Teuchos::Describable::description();
174 if(!is_null(amesosSolver_)) {
176 <<
"{fwdOp="<<fwdOp_->description()
177 <<
",amesosSolver="<<typeName(*amesosSolver_)<<
"}";
184 Teuchos::FancyOStream &out,
185 const Teuchos::EVerbosityLevel verbLevel
188 using Teuchos::OSTab;
189 using Teuchos::typeName;
190 using Teuchos::describe;
192 case Teuchos::VERB_DEFAULT:
193 case Teuchos::VERB_LOW:
196 case Teuchos::VERB_MEDIUM:
197 case Teuchos::VERB_HIGH:
198 case Teuchos::VERB_EXTREME:
201 << Teuchos::Describable::description() <<
"{"
202 <<
"rangeDim=" << this->
range()->dim()
203 <<
",domainDim="<< this->
domain()->dim() <<
"}\n";
205 if(!is_null(fwdOp_)) {
206 out <<
"fwdOp = " <<
describe(*fwdOp_,verbLevel);
208 if(!is_null(amesosSolver_)) {
209 out <<
"amesosSolver=" << typeName(*amesosSolver_) <<
"\n";
214 TEUCHOS_TEST_FOR_EXCEPT(
true);
227 return ::Thyra::opSupported(*fwdOp_,M_trans);
232 const EOpTransp M_trans,
233 const MultiVectorBase<double> &X,
234 const Ptr<MultiVectorBase<double> > &Y,
239 Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
248 if (Thyra::real_trans(M_trans) == Thyra::NOTRANS) {
259 const bool oldUseTranspose = amesosSolver_->UseTranspose();
260 const bool supportsAdjoint = (amesosSolver_->SetUseTranspose(
true) == 0);
261 amesosSolver_->SetUseTranspose(oldUseTranspose);
262 return supportsAdjoint;
267 EOpTransp ,
const SolveMeasureType&
276 const EOpTransp M_trans,
277 const MultiVectorBase<double> &B,
278 const Ptr<MultiVectorBase<double> > &X,
279 const Ptr<
const SolveCriteria<double> >
282 using Teuchos::rcpFromPtr;
283 using Teuchos::rcpFromRef;
284 using Teuchos::OSTab;
286 Teuchos::Time totalTimer(
"");
287 totalTimer.start(
true);
289 THYRA_FUNC_TIME_MONITOR(
"Stratimikos: AmesosLOWS");
291 Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
292 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
293 OSTab tab = this->getOSTab();
294 if(out.get() &&
static_cast<int>(verbLevel) >
static_cast<int>(Teuchos::VERB_NONE))
295 *out <<
"\nSolving block system using Amesos solver "
296 << typeName(*amesosSolver_) <<
" ...\n\n";
301 const EOpTransp amesosOpTransp = real_trans(trans_trans(amesosSolverTransp_,M_trans));
304 &opRangeMap = ( amesosOpTransp ==
NOTRANS
306 &opDomainMap = ( amesosOpTransp ==
NOTRANS
312 Teuchos::RCP<const Epetra_MultiVector>
313 epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
314 Teuchos::RCP<Epetra_MultiVector>
315 epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
320 epetraLP_->SetLHS(&*epetra_X);
327 const bool oldUseTranspose = amesosSolver_->UseTranspose();
328 amesosSolver_->SetUseTranspose(amesosOpTransp==
TRANS);
329 const int err = amesosSolver_->Solve();
330 TEUCHOS_TEST_FOR_EXCEPTION( 0!=err, CatastrophicSolveFailure,
331 "Error, the function Solve() on the amesos solver of type\n"
332 "\'"<<typeName(*amesosSolver_)<<
"\' failed with error code "<<err<<
"!"
334 amesosSolver_->SetUseTranspose(oldUseTranspose);
339 epetraLP_->SetLHS(NULL);
340 epetraLP_->SetRHS(NULL);
341 epetra_X = Teuchos::null;
342 epetra_B = Teuchos::null;
347 if(amesosSolverScalar_!=1.0)
348 Thyra::scale(1.0/amesosSolverScalar_, X);
353 SolveStatus<double> solveStatus;
354 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
355 solveStatus.achievedTol = SolveStatus<double>::unknownTolerance();
356 solveStatus.message =
357 std::string(
"Solver ")+typeName(*amesosSolver_)+std::string(
" converged!");
362 if(out.get() &&
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_LOW))
364 <<
"\nTotal solve time = "<<totalTimer.totalElapsedTime()<<
" sec\n";
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
SolveStatus< double > solveImpl(const EOpTransp M_trans, const MultiVectorBase< double > &B, const Ptr< MultiVectorBase< double > > &X, const Ptr< const SolveCriteria< double > > solveCriteria) const
virtual void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
Teuchos::RCP< const LinearOpBase< double > > clone() const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
void uninitialize(Teuchos::RCP< const LinearOpBase< double > > *fwdOp=NULL, Teuchos::RCP< const LinearOpSourceBase< double > > *fwdOpSrc=NULL, Teuchos::RCP< Epetra_LinearProblem > *epetraLP=NULL, Teuchos::RCP< Amesos_BaseSolver > *amesosSolver=NULL, EOpTransp *amesosSolverTransp=NULL, double *amesosSolverScalar=NULL)
Uninitialize.
Teuchos::RCP< const VectorSpaceBase< double > > domain() const
void initialize(const Teuchos::RCP< const LinearOpBase< double > > &fwdOp, const Teuchos::RCP< const LinearOpSourceBase< double > > &fwdOpSrc, const Teuchos::RCP< Epetra_LinearProblem > &epetraLP, const Teuchos::RCP< Amesos_BaseSolver > &amesosSolver, const EOpTransp amesosSolverTransp, const double amesosSolverScalar)
First initialization.
virtual bool opSupportedImpl(EOpTransp M_trans) const
virtual bool solveSupportsImpl(EOpTransp M_trans) const
std::string description() const
Teuchos::RCP< const VectorSpaceBase< double > > range() const
Teuchos::RCP< const LinearOpSourceBase< double > > extract_fwdOpSrc()
Extract the LinearOpSourceBase<double> object so that it can be modified.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
AmesosLinearOpWithSolve()
Construct to uninitialized.