42 #include "Thyra_AztecOOLinearOpWithSolve.hpp" 43 #include "Thyra_LinearOpWithSolveHelpers.hpp" 44 #include "Thyra_EpetraThyraWrappers.hpp" 45 #include "Thyra_EpetraOperatorWrapper.hpp" 46 #include "Teuchos_BLAS_types.hpp" 47 #include "Teuchos_TimeMonitor.hpp" 48 #include "Teuchos_Time.hpp" 49 #include "Teuchos_implicit_cast.hpp" 57 Teuchos::ETransp convert( Thyra::EOpTransp trans_in )
59 Teuchos::ETransp trans_out;
62 trans_out = Teuchos::NO_TRANS;
65 trans_out = Teuchos::TRANS;
68 TEUCHOS_TEST_FOR_EXCEPT(
true);
79 class SetAztecSolveState {
82 const Teuchos::RCP<AztecOO> &aztecSolver,
83 const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
84 const Teuchos::EVerbosityLevel verbLevel,
85 const Thyra::SolveMeasureType &solveMeasureType
87 ~SetAztecSolveState();
89 Teuchos::RCP<AztecOO> aztecSolver_;
90 Teuchos::RCP<Teuchos::FancyOStream> fancyOStream_;
91 Teuchos::EVerbosityLevel verbLevel_;
98 SetAztecSolveState::SetAztecSolveState(
99 const Teuchos::RCP<AztecOO> &aztecSolver,
100 const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
101 const Teuchos::EVerbosityLevel verbLevel,
102 const Thyra::SolveMeasureType &solveMeasureType
104 :aztecSolver_(aztecSolver.assert_not_null())
109 verbLevel_ = verbLevel;
110 if ( Teuchos::VERB_NONE != verbLevel_ ) {
111 if (!is_null(fancyOStream)) {
115 fancyOStream_= Teuchos::tab(
116 fancyOStream.create_weak(),
118 Teuchos::implicit_cast<std::string>(
"AZTECOO")
120 aztecSolver_->SetOutputStream(*fancyOStream_);
121 aztecSolver_->SetErrorStream(*fancyOStream_);
131 outputFrequency_ = aztecSolver_->GetAllAztecOptions()[AZ_output];
132 aztecSolver_->SetAztecOption(AZ_output,0);
136 convergenceTest_ = aztecSolver_->GetAztecOption(AZ_conv);
137 if (solveMeasureType.useDefault())
143 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
144 Thyra::SOLVE_MEASURE_NORM_RHS
148 aztecSolver_->SetAztecOption(AZ_conv,AZ_rhs);
152 Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
153 Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL
157 aztecSolver_->SetAztecOption(AZ_conv,AZ_r0);
160 TEUCHOS_TEST_FOR_EXCEPT(
"Invalid solve measure type, you should never get here!");
166 SetAztecSolveState::~SetAztecSolveState()
170 if ( Teuchos::VERB_NONE != verbLevel_ ) {
171 if (!is_null(fancyOStream_)) {
172 aztecSolver_->SetOutputStream(std::cout);
173 aztecSolver_->SetErrorStream(std::cerr);
174 *fancyOStream_ <<
"\n";
178 aztecSolver_->SetAztecOption(AZ_output,outputFrequency_);
182 aztecSolver_->SetAztecOption(AZ_conv,convergenceTest_);
197 const int fwdDefaultMaxIterations_in
198 ,
const double fwdDefaultTol_in
199 ,
const int adjDefaultMaxIterations_in
200 ,
const double adjDefaultTol_in
201 ,
const bool outputEveryRhs_in
203 :fwdDefaultMaxIterations_(fwdDefaultMaxIterations_in)
204 ,fwdDefaultTol_(fwdDefaultTol_in)
205 ,adjDefaultMaxIterations_(adjDefaultMaxIterations_in)
206 ,adjDefaultTol_(adjDefaultTol_in)
207 ,outputEveryRhs_(outputEveryRhs_in)
208 ,isExternalPrec_(false)
209 ,allowInexactFwdSolve_(false)
210 ,allowInexactAdjSolve_(false)
211 ,aztecSolverScalar_(0.0)
216 const RCP<
const LinearOpBase<double> > &fwdOp
217 ,
const RCP<
const LinearOpSourceBase<double> > &fwdOpSrc
218 ,
const RCP<
const PreconditionerBase<double> > &prec
219 ,
const bool isExternalPrec_in
220 ,
const RCP<
const LinearOpSourceBase<double> > &approxFwdOpSrc
221 ,
const RCP<AztecOO> &aztecFwdSolver
222 ,
const bool allowInexactFwdSolve
223 ,
const RCP<AztecOO> &aztecAdjSolver
224 ,
const bool allowInexactAdjSolve
225 ,
const double aztecSolverScalar
229 TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
230 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
231 TEUCHOS_TEST_FOR_EXCEPT(aztecFwdSolver.get()==NULL);
234 fwdOpSrc_ = fwdOpSrc;
235 isExternalPrec_ = isExternalPrec_in;
237 approxFwdOpSrc_ = approxFwdOpSrc;
238 aztecFwdSolver_ = aztecFwdSolver;
239 allowInexactFwdSolve_ = allowInexactFwdSolve;
240 aztecAdjSolver_ = aztecAdjSolver;
241 allowInexactAdjSolve_ = allowInexactAdjSolve;
242 aztecSolverScalar_ = aztecSolverScalar;
243 const std::string fwdOpLabel = fwdOp_->getObjectLabel();
244 if (fwdOpLabel.length())
245 this->setObjectLabel(
"lows("+fwdOpLabel+
")" );
249 RCP<const LinearOpSourceBase<double> >
252 RCP<const LinearOpSourceBase<double> >
253 _fwdOpSrc = fwdOpSrc_;
254 fwdOpSrc_ = Teuchos::null;
259 RCP<const PreconditionerBase<double> >
262 RCP<const PreconditionerBase<double> >
264 prec_ = Teuchos::null;
271 return isExternalPrec_;
275 RCP<const LinearOpSourceBase<double> >
278 RCP<const LinearOpSourceBase<double> >
279 _approxFwdOpSrc = approxFwdOpSrc_;
280 approxFwdOpSrc_ = Teuchos::null;
281 return _approxFwdOpSrc;
286 RCP<
const LinearOpBase<double> > *fwdOp,
287 RCP<
const LinearOpSourceBase<double> > *fwdOpSrc,
288 RCP<
const PreconditionerBase<double> > *prec,
289 bool *isExternalPrec_inout,
290 RCP<
const LinearOpSourceBase<double> > *approxFwdOpSrc,
291 RCP<AztecOO> *aztecFwdSolver,
292 bool *allowInexactFwdSolve,
293 RCP<AztecOO> *aztecAdjSolver,
294 bool *allowInexactAdjSolve,
295 double *aztecSolverScalar
298 if (fwdOp) *fwdOp = fwdOp_;
299 if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
300 if (prec) *prec = prec_;
301 if (isExternalPrec_inout) *isExternalPrec_inout = isExternalPrec_;
302 if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
303 if (aztecFwdSolver) *aztecFwdSolver = aztecFwdSolver_;
304 if (allowInexactFwdSolve) *allowInexactFwdSolve = allowInexactFwdSolve_;
305 if (aztecAdjSolver) *aztecAdjSolver = aztecAdjSolver_;
306 if (allowInexactAdjSolve) *allowInexactAdjSolve = allowInexactAdjSolve_;
307 if (aztecSolverScalar) *aztecSolverScalar = aztecSolverScalar_;
309 fwdOp_ = Teuchos::null;
310 fwdOpSrc_ = Teuchos::null;
311 prec_ = Teuchos::null;
312 isExternalPrec_ =
false;
313 approxFwdOpSrc_ = Teuchos::null;
314 aztecFwdSolver_ = Teuchos::null;
315 allowInexactFwdSolve_ =
false;
316 aztecAdjSolver_ = Teuchos::null;
317 allowInexactAdjSolve_ =
false;
318 aztecSolverScalar_ = 0.0;
325 RCP< const VectorSpaceBase<double> >
328 return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
332 RCP< const VectorSpaceBase<double> >
335 return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
339 RCP<const LinearOpBase<double> >
342 return Teuchos::null;
351 std::ostringstream oss;
352 oss << Teuchos::Describable::description();
355 oss <<
"fwdOp="<<fwdOp_->description()<<
"";
363 Teuchos::FancyOStream &out,
364 const Teuchos::EVerbosityLevel verbLevel
367 using Teuchos::OSTab;
368 using Teuchos::typeName;
369 using Teuchos::describe;
371 case Teuchos::VERB_DEFAULT:
372 case Teuchos::VERB_LOW:
375 case Teuchos::VERB_MEDIUM:
376 case Teuchos::VERB_HIGH:
377 case Teuchos::VERB_EXTREME:
380 << Teuchos::Describable::description() <<
"{" 381 <<
"rangeDim=" << this->
range()->dim()
382 <<
",domainDim="<< this->
domain()->dim() <<
"}\n";
384 if (!is_null(fwdOp_)) {
385 out <<
"fwdOp = " <<
describe(*fwdOp_,verbLevel);
387 if (!is_null(prec_)) {
388 out <<
"prec = " <<
describe(*prec_,verbLevel);
390 if (!is_null(aztecFwdSolver_)) {
391 if (aztecFwdSolver_->GetUserOperator())
394 << typeName(*aztecFwdSolver_->GetUserOperator()) <<
"\n";
395 if (aztecFwdSolver_->GetUserMatrix())
397 <<
"Aztec Fwd Mat = " 398 << typeName(*aztecFwdSolver_->GetUserMatrix()) <<
"\n";
399 if (aztecFwdSolver_->GetPrecOperator())
401 <<
"Aztec Fwd Prec Op = " 402 << typeName(*aztecFwdSolver_->GetPrecOperator()) <<
"\n";
403 if (aztecFwdSolver_->GetPrecMatrix())
405 <<
"Aztec Fwd Prec Mat = " 406 << typeName(*aztecFwdSolver_->GetPrecMatrix()) <<
"\n";
408 if (!is_null(aztecAdjSolver_)) {
409 if (aztecAdjSolver_->GetUserOperator())
412 << typeName(*aztecAdjSolver_->GetUserOperator()) <<
"\n";
413 if (aztecAdjSolver_->GetUserMatrix())
415 <<
"Aztec Adj Mat = " 416 << typeName(*aztecAdjSolver_->GetUserMatrix()) <<
"\n";
417 if (aztecAdjSolver_->GetPrecOperator())
419 <<
"Aztec Adj Prec Op = " 420 << typeName(*aztecAdjSolver_->GetPrecOperator()) <<
"\n";
421 if (aztecAdjSolver_->GetPrecMatrix())
423 <<
"Aztec Adj Prec Mat = " 424 << typeName(*aztecAdjSolver_->GetPrecMatrix()) <<
"\n";
429 TEUCHOS_TEST_FOR_EXCEPT(
true);
442 return ::Thyra::opSupported(*fwdOp_,M_trans);
447 const EOpTransp M_trans,
448 const MultiVectorBase<double> &X,
449 const Ptr<MultiVectorBase<double> > &Y,
454 Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
464 if (real_trans(M_trans)==NOTRANS)
return true;
465 return (nonnull(aztecAdjSolver_));
471 EOpTransp M_trans,
const SolveMeasureType& solveMeasureType
474 if (real_trans(M_trans)==NOTRANS) {
475 if (solveMeasureType.useDefault())
481 SOLVE_MEASURE_NORM_RESIDUAL,
482 SOLVE_MEASURE_NORM_RHS
485 allowInexactFwdSolve_
492 SOLVE_MEASURE_NORM_RESIDUAL,
493 SOLVE_MEASURE_NORM_INIT_RESIDUAL
496 allowInexactFwdSolve_
504 if (aztecAdjSolver_.get()==NULL)
508 else if (solveMeasureType.useDefault())
514 SOLVE_MEASURE_NORM_RESIDUAL,
515 SOLVE_MEASURE_NORM_RHS
518 allowInexactFwdSolve_
525 SOLVE_MEASURE_NORM_RESIDUAL,
526 SOLVE_MEASURE_NORM_INIT_RESIDUAL
529 allowInexactFwdSolve_
545 const EOpTransp M_trans,
546 const MultiVectorBase<double> &B,
547 const Ptr<MultiVectorBase<double> > &X,
548 const Ptr<
const SolveCriteria<double> > solveCriteria
553 using Teuchos::rcpFromRef;
554 using Teuchos::rcpFromPtr;
555 using Teuchos::OSTab;
556 typedef SolveCriteria<double> SC;
557 typedef SolveStatus<double> SS;
559 THYRA_FUNC_TIME_MONITOR(
"Stratimikos: AztecOOLOWS");
560 Teuchos::Time totalTimer(
""), timer(
"");
561 totalTimer.start(
true);
563 RCP<Teuchos::FancyOStream> out = this->getOStream();
564 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
565 OSTab tab = this->getOSTab();
566 if (out.get() &&
static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
567 *out <<
"\nSolving block system using AztecOO ...\n\n";
573 SolveMeasureType solveMeasureType;
574 if (nonnull(solveCriteria)) {
575 solveMeasureType = solveCriteria->solveMeasureType;
576 assertSupportsSolveMeasureType(*
this, M_trans, solveMeasureType);
582 const EOpTransp aztecOpTransp = real_trans(M_trans);
588 aztecSolver = ( aztecOpTransp == NOTRANS ? aztecFwdSolver_ : aztecAdjSolver_ );
589 const Epetra_Operator
590 *aztecOp = aztecSolver->GetUserOperator();
596 &opRangeMap = aztecOp->OperatorRangeMap(),
597 &opDomainMap = aztecOp->OperatorDomainMap();
602 double tol = ( aztecOpTransp==NOTRANS ? fwdDefaultTol() : adjDefaultTol() );
603 int maxIterations = ( aztecOpTransp==NOTRANS
604 ? fwdDefaultMaxIterations() : adjDefaultMaxIterations() );
605 bool isDefaultSolveCriteria =
true;
606 if (nonnull(solveCriteria)) {
607 if ( solveCriteria->requestedTol != SC::unspecifiedTolerance() ) {
608 tol = solveCriteria->requestedTol;
609 isDefaultSolveCriteria =
false;
611 if (nonnull(solveCriteria->extraParameters)) {
612 maxIterations = solveCriteria->extraParameters->get(
"Maximum Iterations",maxIterations);
620 RCP<const Epetra_MultiVector> epetra_B;
621 RCP<Epetra_MultiVector> epetra_X;
623 const EpetraOperatorWrapper* opWrapper =
624 dynamic_cast<const EpetraOperatorWrapper*
>(aztecOp);
626 if (opWrapper == 0) {
627 epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
628 epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
635 int totalIterations = 0;
636 SolveStatus<double> solveStatus;
637 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
638 solveStatus.achievedTol = -1.0;
644 const int m = B.domain()->dim();
646 for(
int j = 0; j < m; ++j ) {
648 THYRA_FUNC_TIME_MONITOR_DIFF(
"Stratimikos: AztecOOLOWS:SingleSolve", SingleSolve);
658 RCP<Epetra_Vector> epetra_b_j;
659 RCP<Epetra_Vector> epetra_x_j;
661 if (opWrapper == 0) {
662 epetra_b_j = rcpFromRef(*const_cast<Epetra_Vector*>((*epetra_B)(j)));
663 epetra_x_j = rcpFromRef(*(*epetra_X)(j));
666 if (is_null(epetra_b_j)) {
667 epetra_b_j = rcp(
new Epetra_Vector(opRangeMap));
668 epetra_x_j = rcp(
new Epetra_Vector(opDomainMap));
670 opWrapper->copyThyraIntoEpetra(*B.col(j), *epetra_b_j);
671 opWrapper->copyThyraIntoEpetra(*X->col(j), *epetra_x_j);
678 aztecSolver->SetRHS(&*epetra_b_j);
679 aztecSolver->SetLHS(&*epetra_x_j);
687 setAztecSolveState(aztecSolver,out,verbLevel,solveMeasureType);
688 aztecSolver->Iterate( maxIterations, tol );
699 if (aztecSolverScalar_ != 1.0)
700 epetra_x_j->Scale(1.0/aztecSolverScalar_);
705 if (opWrapper != 0) {
706 opWrapper->copyEpetraIntoThyra(*epetra_x_j, X->col(j).ptr());
713 const int iterations = aztecSolver->NumIters();
714 const double achievedTol = aztecSolver->ScaledResidual();
715 const double *AZ_status = aztecSolver->GetAztecStatus();
716 std::ostringstream oss;
717 bool converged =
false;
718 if (AZ_status[AZ_why]==AZ_normal) { oss <<
"Aztec returned AZ_normal."; converged =
true; }
719 else if (AZ_status[AZ_why]==AZ_param) oss <<
"Aztec returned AZ_param.";
720 else if (AZ_status[AZ_why]==AZ_breakdown) oss <<
"Aztec returned AZ_breakdown.";
721 else if (AZ_status[AZ_why]==AZ_loss) oss <<
"Aztec returned AZ_loss.";
722 else if (AZ_status[AZ_why]==AZ_ill_cond) oss <<
"Aztec returned AZ_ill_cond.";
723 else if (AZ_status[AZ_why]==AZ_maxits) oss <<
"Aztec returned AZ_maxits.";
724 else oss <<
"Aztec returned an unknown status?";
725 oss <<
" Iterations = " << iterations <<
".";
726 oss <<
" Achieved Tolerance = " << achievedTol <<
".";
727 oss <<
" Total time = " << timer.totalElapsedTime() <<
" sec.";
728 if (out.get() &&
static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE) && outputEveryRhs())
729 Teuchos::OSTab(out).o() <<
"j="<<j<<
": " << oss.str() <<
"\n";
731 solveStatus.achievedTol = TEUCHOS_MAX(solveStatus.achievedTol, achievedTol);
734 totalIterations += iterations;
736 solveStatus.message = oss.str();
737 if ( isDefaultSolveCriteria ) {
738 switch(solveStatus.solveStatus) {
739 case SOLVE_STATUS_UNKNOWN:
742 case SOLVE_STATUS_CONVERGED:
743 solveStatus.solveStatus = ( converged ? SOLVE_STATUS_CONVERGED : SOLVE_STATUS_UNCONVERGED );
745 case SOLVE_STATUS_UNCONVERGED:
749 TEUCHOS_TEST_FOR_EXCEPT(
true);
754 aztecSolver->UnsetLHSRHS();
759 epetra_X = Teuchos::null;
760 epetra_B = Teuchos::null;
766 SolveStatus<double> overallSolveStatus;
767 if (isDefaultSolveCriteria) {
768 overallSolveStatus.solveStatus = SOLVE_STATUS_UNKNOWN;
769 overallSolveStatus.achievedTol = SS::unknownTolerance();
772 overallSolveStatus.solveStatus = solveStatus.solveStatus;
773 overallSolveStatus.achievedTol = solveStatus.achievedTol;
775 std::ostringstream oss;
778 << ( overallSolveStatus.solveStatus==SOLVE_STATUS_CONVERGED ?
"converged" :
"unconverged" )
779 <<
" on m = "<<m<<
" RHSs using " << totalIterations <<
" cumulative iterations" 780 <<
" for an average of " << (totalIterations/m) <<
" iterations/RHS and" 781 <<
" total CPU time of "<<totalTimer.totalElapsedTime()<<
" sec.";
782 overallSolveStatus.message = oss.str();
785 if (overallSolveStatus.extraParameters.is_null()) {
786 overallSolveStatus.extraParameters = Teuchos::parameterList ();
788 overallSolveStatus.extraParameters->set (
"AztecOO/Iteration Count",
791 overallSolveStatus.extraParameters->set (
"Iteration Count",
793 overallSolveStatus.extraParameters->set (
"AztecOO/Achieved Tolerance",
794 overallSolveStatus.achievedTol);
799 if (out.get() &&
static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
801 <<
"\nTotal solve time = "<<totalTimer.totalElapsedTime()<<
" sec\n";
803 return overallSolveStatus;
RCP< const LinearOpBase< double > > clone() const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
RCP< const LinearOpSourceBase< double > > extract_fwdOpSrc()
Extract the forward LinearOpBase<double> object so that it can be modified.
RCP< const VectorSpaceBase< double > > domain() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
bool isExternalPrec() const
Determine if the preconditioner was external or not.
virtual bool solveSupportsImpl(EOpTransp M_trans) const
SolveStatus< double > solveImpl(const EOpTransp M_trans, const MultiVectorBase< double > &B, const Ptr< MultiVectorBase< double > > &X, const Ptr< const SolveCriteria< double > > solveCriteria) const
void uninitialize(RCP< const LinearOpBase< double > > *fwdOp=NULL, RCP< const LinearOpSourceBase< double > > *fwdOpSrc=NULL, RCP< const PreconditionerBase< double > > *prec=NULL, bool *isExternalPrec=NULL, RCP< const LinearOpSourceBase< double > > *approxFwdOpSrc=NULL, RCP< AztecOO > *aztecFwdSolver=NULL, bool *allowInexactFwdSolve=NULL, RCP< AztecOO > *aztecAdjSolver=NULL, bool *allowInexactAdjSolve=NULL, double *aztecSolverScalar=NULL)
Uninitialize.
RCP< const PreconditionerBase< double > > extract_prec()
Extract the preconditioner.
std::string description() const
RCP< const VectorSpaceBase< double > > range() const
virtual void applyImpl(const EOpTransp M_trans, const MultiVectorBase< double > &X, const Ptr< MultiVectorBase< double > > &Y, const double alpha, const double beta) const
void initialize(const RCP< const LinearOpBase< double > > &fwdOp, const RCP< const LinearOpSourceBase< double > > &fwdOpSrc, const RCP< const PreconditionerBase< double > > &prec, const bool isExternalPrec, const RCP< const LinearOpSourceBase< double > > &approxFwdOpSrc, const RCP< AztecOO > &aztecFwdSolver, const bool allowInexactFwdSolve=false, const RCP< AztecOO > &aztecAdjSolver=Teuchos::null, const bool allowInexactAdjSolve=false, const double aztecSolverScalar=1.0)
Sets up this object.
virtual bool opSupportedImpl(EOpTransp M_trans) const
RCP< const LinearOpSourceBase< double > > extract_approxFwdOpSrc()
Extract the approximate forward LinearOpBase<double> object used to build the preconditioner.
AztecOOLinearOpWithSolve(const int fwdDefaultMaxIterations=400, const double fwdDefaultTol=1e-6, const int adjDefaultMaxIterations=400, const double adjDefaultTol=1e-6, const bool outputEveryRhs=false)