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);
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 );
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;