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+
")" );
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;
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);
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_ );
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)) {
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;