Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_AztecOOLinearOpWithSolve.cpp
Go to the documentation of this file.
00001 /*@HEADER
00002 // ***********************************************************************
00003 // 
00004 //        AztecOO: An Object-Oriented Aztec Linear Solver Package 
00005 //                 Copyright (2002) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 */
00029 
00030 #include "Thyra_AztecOOLinearOpWithSolve.hpp"
00031 #include "Thyra_LinearOpWithSolveHelpers.hpp"
00032 #include "Thyra_EpetraThyraWrappers.hpp"
00033 #include "Thyra_EpetraOperatorWrapper.hpp"
00034 #include "Teuchos_BLAS_types.hpp"
00035 #include "Teuchos_TimeMonitor.hpp"
00036 #include "Teuchos_Time.hpp"
00037 #include "Teuchos_implicit_cast.hpp"
00038 
00039 
00040 namespace {
00041 
00042 
00043 inline
00044 Teuchos::ETransp convert( Thyra::EOpTransp trans_in )
00045 {
00046   Teuchos::ETransp trans_out;
00047   switch(trans_in) {
00048     case Thyra::NOTRANS:
00049       trans_out = Teuchos::NO_TRANS;
00050       break;
00051     case Thyra::TRANS:
00052       trans_out = Teuchos::TRANS;
00053       break;
00054     default:
00055       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00056   }
00057   return trans_out;
00058 }
00059 
00060 
00061 // This class sets some solve instance specific state and then sets it back to
00062 // the default state on destruction.  But using the destructor to unset the
00063 // state we can be sure that the state is rest correctly even if an exception
00064 // is thrown.
00065 class SetAztecSolveState {
00066 public:
00067   SetAztecSolveState( 
00068     const Teuchos::RCP<AztecOO> &aztecSolver,
00069     const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
00070     const Teuchos::EVerbosityLevel verbLevel,
00071     const Thyra::SolveMeasureType &solveMeasureType
00072     );
00073   ~SetAztecSolveState();
00074 private:
00075   Teuchos::RCP<AztecOO> aztecSolver_;
00076   Teuchos::RCP<Teuchos::FancyOStream> fancyOStream_;
00077   Teuchos::EVerbosityLevel verbLevel_;
00078   int outputFrequency_;
00079   int convergenceTest_;
00080   SetAztecSolveState(); // Not defined and not to be called!
00081 };
00082 
00083 
00084 SetAztecSolveState::SetAztecSolveState( 
00085   const Teuchos::RCP<AztecOO> &aztecSolver,
00086   const Teuchos::RCP<Teuchos::FancyOStream> &fancyOStream,
00087   const Teuchos::EVerbosityLevel verbLevel,
00088   const Thyra::SolveMeasureType &solveMeasureType
00089   )
00090   :aztecSolver_(aztecSolver.assert_not_null())
00091 {
00092   
00093   // Output state
00094   verbLevel_ = verbLevel;
00095   if ( Teuchos::VERB_NONE != verbLevel_ ) {
00096     if (!is_null(fancyOStream)) {
00097       // AztecOO puts in two tabs before it prints anything.  Therefore,
00098       // there is not much that we can do to improve the layout of the
00099       // indentation so just leave it!
00100       fancyOStream_= Teuchos::tab(
00101         fancyOStream.create_weak(),
00102         0, // Don't indent since AztecOO already puts in two tabs (not spaces!)
00103         Teuchos::implicit_cast<std::string>("AZTECOO")
00104         );
00105       aztecSolver_->SetOutputStream(*fancyOStream_);
00106       aztecSolver_->SetErrorStream(*fancyOStream_);
00107       // Note, above we can not save the current output and error streams
00108       // since AztecOO does not define functions to get them.  In the
00109       // future, AztecOO should define these functions if we are to avoid
00110       // treading on each others print statements.  However, since the
00111       // AztecOO object is most likely owned by these Thyra wrappers, this
00112       // should not be a problem.
00113     }
00114   }
00115   else {
00116     outputFrequency_ = aztecSolver_->GetAllAztecOptions()[AZ_output];
00117     aztecSolver_->SetAztecOption(AZ_output,0);
00118   }
00119 
00120   // Convergence test
00121   convergenceTest_ = aztecSolver_->GetAztecOption(AZ_conv);
00122   if (solveMeasureType.useDefault())
00123   {
00124     // Just use the default solve measure type already set!
00125   }
00126   else if (
00127     solveMeasureType(
00128       Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
00129       Thyra::SOLVE_MEASURE_NORM_RHS
00130       )
00131     )
00132   {
00133     aztecSolver_->SetAztecOption(AZ_conv,AZ_rhs);
00134   }
00135   else if (
00136     solveMeasureType(
00137       Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
00138       Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL
00139       )
00140     )
00141   {
00142     aztecSolver_->SetAztecOption(AZ_conv,AZ_r0);
00143   }
00144   else {
00145     TEUCHOS_TEST_FOR_EXCEPT("Invalid solve measure type, you should never get here!");
00146   }
00147 
00148 }
00149 
00150 
00151 SetAztecSolveState::~SetAztecSolveState()
00152 {
00153       
00154   // Output state
00155   if ( Teuchos::VERB_NONE != verbLevel_ ) {
00156     if (!is_null(fancyOStream_)) {
00157       aztecSolver_->SetOutputStream(std::cout);
00158       aztecSolver_->SetErrorStream(std::cerr);
00159       *fancyOStream_ << "\n";
00160     }
00161   }
00162   else {
00163     aztecSolver_->SetAztecOption(AZ_output,outputFrequency_);
00164   }
00165 
00166   // Convergence test
00167   aztecSolver_->SetAztecOption(AZ_conv,convergenceTest_);
00168       
00169 }
00170 
00171 
00172 } // namespace
00173 
00174 
00175 namespace Thyra {
00176 
00177 
00178 // Constructors/initializers/accessors
00179 
00180 
00181 AztecOOLinearOpWithSolve::AztecOOLinearOpWithSolve(
00182   const int fwdDefaultMaxIterations_in
00183   ,const double fwdDefaultTol_in
00184   ,const int adjDefaultMaxIterations_in
00185   ,const double adjDefaultTol_in
00186   ,const bool outputEveryRhs_in
00187   )
00188   :fwdDefaultMaxIterations_(fwdDefaultMaxIterations_in)
00189   ,fwdDefaultTol_(fwdDefaultTol_in)
00190   ,adjDefaultMaxIterations_(adjDefaultMaxIterations_in)
00191   ,adjDefaultTol_(adjDefaultTol_in)
00192   ,outputEveryRhs_(outputEveryRhs_in)
00193   ,isExternalPrec_(false)
00194   ,allowInexactFwdSolve_(false)
00195   ,allowInexactAdjSolve_(false)
00196   ,aztecSolverScalar_(0.0)
00197 {}
00198 
00199 
00200 void AztecOOLinearOpWithSolve::initialize(
00201   const RCP<const LinearOpBase<double> > &fwdOp
00202   ,const RCP<const LinearOpSourceBase<double> > &fwdOpSrc
00203   ,const RCP<const PreconditionerBase<double> > &prec
00204   ,const bool isExternalPrec_in
00205   ,const RCP<const LinearOpSourceBase<double> > &approxFwdOpSrc
00206   ,const RCP<AztecOO> &aztecFwdSolver
00207   ,const bool allowInexactFwdSolve
00208   ,const RCP<AztecOO> &aztecAdjSolver
00209   ,const bool allowInexactAdjSolve
00210   ,const double aztecSolverScalar
00211   )
00212 {
00213 #ifdef TEUCHOS_DEBUG
00214   TEUCHOS_TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00215   TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00216   TEUCHOS_TEST_FOR_EXCEPT(aztecFwdSolver.get()==NULL);
00217 #endif
00218   fwdOp_ = fwdOp;
00219   fwdOpSrc_ = fwdOpSrc;
00220   isExternalPrec_ = isExternalPrec_in;
00221   prec_ = prec;
00222   approxFwdOpSrc_ = approxFwdOpSrc;
00223   aztecFwdSolver_ = aztecFwdSolver;
00224   allowInexactFwdSolve_ = allowInexactFwdSolve;
00225   aztecAdjSolver_ = aztecAdjSolver;
00226   allowInexactAdjSolve_ = allowInexactAdjSolve;
00227   aztecSolverScalar_ = aztecSolverScalar;
00228   const std::string fwdOpLabel = fwdOp_->getObjectLabel();
00229   if (fwdOpLabel.length())
00230     this->setObjectLabel( "lows("+fwdOpLabel+")" );
00231 }
00232 
00233 
00234 RCP<const LinearOpSourceBase<double> >
00235 AztecOOLinearOpWithSolve::extract_fwdOpSrc()
00236 {
00237   RCP<const LinearOpSourceBase<double> >
00238     _fwdOpSrc = fwdOpSrc_;
00239   fwdOpSrc_ = Teuchos::null;
00240   return _fwdOpSrc;
00241 }
00242 
00243 
00244 RCP<const PreconditionerBase<double> >
00245 AztecOOLinearOpWithSolve::extract_prec()
00246 {
00247   RCP<const PreconditionerBase<double> >
00248     _prec = prec_;
00249   prec_ = Teuchos::null;
00250   return _prec;
00251 }
00252 
00253 
00254 bool AztecOOLinearOpWithSolve::isExternalPrec() const
00255 {
00256   return isExternalPrec_;
00257 }
00258 
00259 
00260 RCP<const LinearOpSourceBase<double> >
00261 AztecOOLinearOpWithSolve::extract_approxFwdOpSrc()
00262 {
00263   RCP<const LinearOpSourceBase<double> >
00264     _approxFwdOpSrc = approxFwdOpSrc_;
00265   approxFwdOpSrc_ = Teuchos::null;
00266   return _approxFwdOpSrc;
00267 }
00268 
00269 
00270 void AztecOOLinearOpWithSolve::uninitialize(
00271   RCP<const LinearOpBase<double> > *fwdOp,
00272   RCP<const LinearOpSourceBase<double> > *fwdOpSrc,
00273   RCP<const PreconditionerBase<double> > *prec,
00274   bool *isExternalPrec_inout,
00275   RCP<const LinearOpSourceBase<double> > *approxFwdOpSrc,
00276   RCP<AztecOO> *aztecFwdSolver,
00277   bool *allowInexactFwdSolve,
00278   RCP<AztecOO> *aztecAdjSolver,
00279   bool *allowInexactAdjSolve,
00280   double *aztecSolverScalar
00281   )
00282 {
00283   if (fwdOp) *fwdOp = fwdOp_;
00284   if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00285   if (prec) *prec = prec_;
00286   if (isExternalPrec_inout) *isExternalPrec_inout = isExternalPrec_;
00287   if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
00288   if (aztecFwdSolver) *aztecFwdSolver = aztecFwdSolver_;
00289   if (allowInexactFwdSolve) *allowInexactFwdSolve = allowInexactFwdSolve_;
00290   if (aztecAdjSolver) *aztecAdjSolver = aztecAdjSolver_;
00291   if (allowInexactAdjSolve) *allowInexactAdjSolve = allowInexactAdjSolve_;
00292   if (aztecSolverScalar) *aztecSolverScalar = aztecSolverScalar_;
00293 
00294   fwdOp_ = Teuchos::null;
00295   fwdOpSrc_ = Teuchos::null;
00296   prec_ = Teuchos::null;
00297   isExternalPrec_ = false; // Just to make unique
00298   approxFwdOpSrc_ = Teuchos::null;
00299   aztecFwdSolver_ = Teuchos::null;
00300   allowInexactFwdSolve_ = false;
00301   aztecAdjSolver_ = Teuchos::null;
00302   allowInexactAdjSolve_ = false;
00303   aztecSolverScalar_ = 0.0;
00304 }
00305 
00306 
00307 // Overridden from LinearOpBase
00308 
00309 
00310 RCP< const VectorSpaceBase<double> >
00311 AztecOOLinearOpWithSolve::range() const
00312 {
00313   return ( fwdOp_.get() ? fwdOp_->range() : Teuchos::null );
00314 }
00315 
00316 
00317 RCP< const VectorSpaceBase<double> >
00318 AztecOOLinearOpWithSolve::domain() const
00319 {
00320   return ( fwdOp_.get() ? fwdOp_->domain() : Teuchos::null );
00321 }
00322 
00323 
00324 RCP<const LinearOpBase<double> >
00325 AztecOOLinearOpWithSolve::clone() const
00326 {
00327   return Teuchos::null; // Not supported yet but could be
00328 }
00329 
00330 
00331 // Overridden from Teuchos::Describable
00332 
00333 
00334 std::string AztecOOLinearOpWithSolve::description() const
00335 {
00336   std::ostringstream oss;
00337   oss << Teuchos::Describable::description();
00338   if (fwdOp_.get()) {
00339     oss << "{";
00340     oss << "fwdOp="<<fwdOp_->description()<<"";
00341     oss << "}";
00342   }
00343   return oss.str();
00344 }
00345 
00346 
00347 void AztecOOLinearOpWithSolve::describe(
00348   Teuchos::FancyOStream &out,
00349   const Teuchos::EVerbosityLevel verbLevel
00350   ) const
00351 {
00352   using Teuchos::OSTab;
00353   using Teuchos::typeName;
00354   using Teuchos::describe;
00355   switch(verbLevel) {
00356     case Teuchos::VERB_DEFAULT:
00357     case Teuchos::VERB_LOW:
00358       out << this->description() << std::endl;
00359       break;
00360     case Teuchos::VERB_MEDIUM:
00361     case Teuchos::VERB_HIGH:
00362     case Teuchos::VERB_EXTREME:
00363     {
00364       out
00365         << Teuchos::Describable::description() << "{"
00366         << "rangeDim=" << this->range()->dim()
00367         << ",domainDim="<< this->domain()->dim() << "}\n";
00368       OSTab tab(out);
00369       if (!is_null(fwdOp_)) {
00370         out << "fwdOp = " << describe(*fwdOp_,verbLevel);
00371       }
00372       if (!is_null(prec_)) {
00373         out << "prec = " << describe(*prec_,verbLevel);
00374       }
00375       if (!is_null(aztecFwdSolver_)) {
00376         if (aztecFwdSolver_->GetUserOperator())
00377           out
00378             << "Aztec Fwd Op = "
00379             << typeName(*aztecFwdSolver_->GetUserOperator()) << "\n";
00380         if (aztecFwdSolver_->GetUserMatrix())
00381           out
00382             << "Aztec Fwd Mat = "
00383             << typeName(*aztecFwdSolver_->GetUserMatrix()) << "\n";
00384         if (aztecFwdSolver_->GetPrecOperator())
00385           out
00386             << "Aztec Fwd Prec Op = "
00387             << typeName(*aztecFwdSolver_->GetPrecOperator()) << "\n";
00388         if (aztecFwdSolver_->GetPrecMatrix())
00389           out
00390             << "Aztec Fwd Prec Mat = "
00391             << typeName(*aztecFwdSolver_->GetPrecMatrix()) << "\n";
00392       }
00393       if (!is_null(aztecAdjSolver_)) {
00394         if (aztecAdjSolver_->GetUserOperator())
00395           out
00396             << "Aztec Adj Op = "
00397             << typeName(*aztecAdjSolver_->GetUserOperator()) << "\n";
00398         if (aztecAdjSolver_->GetUserMatrix())
00399           out
00400             << "Aztec Adj Mat = "
00401             << typeName(*aztecAdjSolver_->GetUserMatrix()) << "\n";
00402         if (aztecAdjSolver_->GetPrecOperator())
00403           out
00404             << "Aztec Adj Prec Op = "
00405             << typeName(*aztecAdjSolver_->GetPrecOperator()) << "\n";
00406         if (aztecAdjSolver_->GetPrecMatrix())
00407           out
00408             << "Aztec Adj Prec Mat = "
00409             << typeName(*aztecAdjSolver_->GetPrecMatrix()) << "\n";
00410       }
00411       break;
00412     }
00413     default:
00414       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00415   }
00416 }
00417 
00418 
00419 // protected
00420 
00421 
00422 // Overridden from LinearOpBase
00423 
00424 
00425 bool AztecOOLinearOpWithSolve::opSupportedImpl(EOpTransp M_trans) const
00426 {
00427   return ::Thyra::opSupported(*fwdOp_,M_trans);
00428 }
00429 
00430 
00431 void AztecOOLinearOpWithSolve::applyImpl(
00432   const EOpTransp M_trans,
00433   const MultiVectorBase<double> &X,
00434   const Ptr<MultiVectorBase<double> > &Y,
00435   const double alpha,
00436   const double beta
00437   ) const
00438 {
00439   Thyra::apply( *fwdOp_, M_trans, X, Y, alpha, beta );
00440 }
00441 
00442 
00443 // Overridden from LinearOpWithSolveBase
00444 
00445 
00446 bool
00447 AztecOOLinearOpWithSolve::solveSupportsImpl(EOpTransp M_trans) const
00448 {
00449   if (real_trans(M_trans)==NOTRANS) return true;
00450   return (nonnull(aztecAdjSolver_));
00451 }
00452 
00453 
00454 bool
00455 AztecOOLinearOpWithSolve::solveSupportsSolveMeasureTypeImpl(
00456   EOpTransp M_trans, const SolveMeasureType& solveMeasureType
00457   ) const
00458 {
00459   if (real_trans(M_trans)==NOTRANS) {
00460     if (solveMeasureType.useDefault())
00461     {
00462       return true;
00463     }
00464     else if (
00465       solveMeasureType(
00466         SOLVE_MEASURE_NORM_RESIDUAL,
00467         SOLVE_MEASURE_NORM_RHS
00468         )
00469       &&
00470       allowInexactFwdSolve_
00471       )
00472     {
00473       return true;
00474     }
00475     else if (
00476       solveMeasureType(
00477         SOLVE_MEASURE_NORM_RESIDUAL,
00478         SOLVE_MEASURE_NORM_INIT_RESIDUAL
00479         )
00480       &&
00481       allowInexactFwdSolve_
00482       )
00483     {
00484       return true;
00485     }
00486   }
00487   else {
00488     // TRANS
00489     if (aztecAdjSolver_.get()==NULL)
00490     {
00491       return false;
00492     }
00493     else if (solveMeasureType.useDefault())
00494     {
00495       return true;
00496     }
00497     else if (
00498       solveMeasureType(
00499         SOLVE_MEASURE_NORM_RESIDUAL,
00500         SOLVE_MEASURE_NORM_RHS
00501         )
00502       &&
00503       allowInexactFwdSolve_
00504       )
00505     {
00506       return true;
00507     }
00508     else if (
00509       solveMeasureType(
00510         SOLVE_MEASURE_NORM_RESIDUAL,
00511         SOLVE_MEASURE_NORM_INIT_RESIDUAL
00512         )
00513       &&
00514       allowInexactFwdSolve_
00515       )
00516     {
00517       return true;
00518     }
00519   }
00520   // If you get here then we don't support the solve measure type!
00521   return false;
00522 }
00523 
00524 
00525 // Overridden from LinearOpWithSolveBase
00526 
00527 
00528 SolveStatus<double>
00529 AztecOOLinearOpWithSolve::solveImpl(
00530   const EOpTransp M_trans,
00531   const MultiVectorBase<double> &B,
00532   const Ptr<MultiVectorBase<double> > &X,
00533   const Ptr<const SolveCriteria<double> > solveCriteria
00534   ) const
00535 {
00536 
00537   using Teuchos::rcp;
00538   using Teuchos::rcpFromRef;
00539   using Teuchos::rcpFromPtr;
00540   using Teuchos::OSTab;
00541   typedef SolveCriteria<double> SC;
00542   typedef SolveStatus<double> SS;
00543 
00544   THYRA_FUNC_TIME_MONITOR("Stratimikos: AztecOOLOWS");
00545   Teuchos::Time totalTimer(""), timer("");
00546   totalTimer.start(true);
00547 
00548   RCP<Teuchos::FancyOStream> out = this->getOStream();
00549   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00550   OSTab tab = this->getOSTab();
00551   if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00552     *out << "\nSolving block system using AztecOO ...\n\n";
00553 
00554   //
00555   // Validate input
00556   //
00557   TEUCHOS_ASSERT(this->solveSupportsImpl(M_trans));
00558   SolveMeasureType solveMeasureType;
00559   if (nonnull(solveCriteria)) {
00560     solveMeasureType = solveCriteria->solveMeasureType;
00561     assertSupportsSolveMeasureType(*this, M_trans, solveMeasureType);
00562   }
00563 
00564   //
00565   // Get the transpose argument
00566   //
00567   const EOpTransp aztecOpTransp = real_trans(M_trans);
00568 
00569   //
00570   // Get the solver, operator, and preconditioner that we will use
00571   //
00572   RCP<AztecOO>
00573     aztecSolver = ( aztecOpTransp == NOTRANS ? aztecFwdSolver_  : aztecAdjSolver_ );
00574   const Epetra_Operator
00575     *aztecOp = aztecSolver->GetUserOperator();
00576 
00577   //
00578   // Get the op(...) range and domain maps
00579   //
00580   const Epetra_Map
00581     &opRangeMap = aztecOp->OperatorRangeMap(),
00582     &opDomainMap = aztecOp->OperatorDomainMap();
00583 
00584   //
00585   // Get the convergence criteria
00586   //
00587   double tol = ( aztecOpTransp==NOTRANS ? fwdDefaultTol() : adjDefaultTol() );
00588   int maxIterations = ( aztecOpTransp==NOTRANS
00589     ? fwdDefaultMaxIterations() : adjDefaultMaxIterations() );
00590   bool isDefaultSolveCriteria = true;
00591   if (nonnull(solveCriteria)) {
00592     if ( solveCriteria->requestedTol != SC::unspecifiedTolerance() ) {
00593       tol = solveCriteria->requestedTol;
00594       isDefaultSolveCriteria = false;
00595     }
00596     if (nonnull(solveCriteria->extraParameters)) {
00597       maxIterations = solveCriteria->extraParameters->get("Maximum Iterations",maxIterations);
00598     }
00599   }
00600 
00601   //
00602   // Get Epetra_MultiVector views of B and X
00603   //
00604 
00605   RCP<const Epetra_MultiVector> epetra_B;
00606   RCP<Epetra_MultiVector> epetra_X;
00607 
00608   const EpetraOperatorWrapper* opWrapper =
00609     dynamic_cast<const EpetraOperatorWrapper*>(aztecOp);
00610 
00611   if (opWrapper == 0) {
00612     epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
00613     epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
00614   }
00615 
00616   //
00617   // Use AztecOO to solve each RHS one at a time (which is all that I can do anyway)
00618   //
00619 
00620   int totalIterations = 0;
00621   SolveStatus<double> solveStatus;
00622   solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00623   solveStatus.achievedTol = -1.0;
00624 
00625   /* Get the number of columns in the multivector. We use Thyra
00626    * functions rather than Epetra functions to do this, as we
00627    * might not yet have created an Epetra multivector. - KL */
00628   //const int m = epetra_B->NumVectors();
00629   const int m = B.domain()->dim();
00630 
00631   for( int j = 0; j < m; ++j ) {
00632 
00633     THYRA_FUNC_TIME_MONITOR_DIFF("Stratimikos: AztecOOLOWS:SingleSolve", SingleSolve);
00634 
00635     //
00636     // Get Epetra_Vector views of B(:,j) and X(:,j)
00637     // How this is done will depend on whether we have a true Epetra operator
00638     // or we are wrapping a general Thyra operator in an Epetra operator.
00639     //
00640 
00641     // We need to declare epetra_x_j as non-const because when we have a phony
00642     // Epetra operator we'll have to copy a thyra vector into it.
00643     RCP<Epetra_Vector> epetra_b_j;
00644     RCP<Epetra_Vector> epetra_x_j;
00645 
00646     if (opWrapper == 0) {
00647       epetra_b_j = rcpFromRef(*const_cast<Epetra_Vector*>((*epetra_B)(j)));
00648       epetra_x_j = rcpFromRef(*(*epetra_X)(j));
00649     }
00650     else {
00651       if (is_null(epetra_b_j)) {
00652         epetra_b_j = rcp(new Epetra_Vector(opRangeMap));
00653         epetra_x_j = rcp(new Epetra_Vector(opDomainMap));
00654       }
00655       opWrapper->copyThyraIntoEpetra(*B.col(j), *epetra_b_j);
00656       opWrapper->copyThyraIntoEpetra(*X->col(j), *epetra_x_j);
00657     }
00658 
00659     //
00660     // Set the RHS and LHS
00661     //
00662 
00663     aztecSolver->SetRHS(&*epetra_b_j);
00664     aztecSolver->SetLHS(&*epetra_x_j);
00665 
00666     //
00667     // Solve the linear system
00668     //
00669     timer.start(true);
00670     {
00671       SetAztecSolveState
00672         setAztecSolveState(aztecSolver,out,verbLevel,solveMeasureType);
00673       aztecSolver->Iterate( maxIterations, tol );
00674       // NOTE: We ignore the returned status but get it below
00675     }
00676     timer.stop();
00677 
00678     //
00679     // Scale the solution 
00680     // (Originally, this was at the end of the loop after all columns had been
00681     // processed. It's moved here because we need to do it before copying the
00682     // solution back into a Thyra vector. - KL
00683     //
00684     if (aztecSolverScalar_ != 1.0)
00685       epetra_x_j->Scale(1.0/aztecSolverScalar_);
00686 
00687     //
00688     // If necessary, convert the solution back to a non-epetra vector
00689     //
00690     if (opWrapper != 0) {
00691       opWrapper->copyEpetraIntoThyra(*epetra_x_j, X->col(j).ptr());
00692     }
00693 
00694     //
00695     // Set the return solve status
00696     //
00697 
00698     const int iterations = aztecSolver->NumIters();
00699     const double achievedTol = aztecSolver->ScaledResidual();
00700     const double *AZ_status = aztecSolver->GetAztecStatus();
00701     std::ostringstream oss;
00702     bool converged = false;
00703     if (AZ_status[AZ_why]==AZ_normal) { oss << "Aztec returned AZ_normal."; converged = true; }
00704     else if (AZ_status[AZ_why]==AZ_param) oss << "Aztec returned AZ_param.";
00705     else if (AZ_status[AZ_why]==AZ_breakdown) oss << "Aztec returned AZ_breakdown.";
00706     else if (AZ_status[AZ_why]==AZ_loss) oss << "Aztec returned AZ_loss.";
00707     else if (AZ_status[AZ_why]==AZ_ill_cond) oss << "Aztec returned AZ_ill_cond.";
00708     else if (AZ_status[AZ_why]==AZ_maxits) oss << "Aztec returned AZ_maxits.";
00709     else oss << "Aztec returned an unknown status?";
00710     oss << "  Iterations = " << iterations << ".";
00711     oss << "  Achieved Tolerance = " << achievedTol << ".";
00712     oss << "  Total time = " << timer.totalElapsedTime() << " sec.";
00713     if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE) && outputEveryRhs())
00714       Teuchos::OSTab(out).o() << "j="<<j<<": " << oss.str() << "\n";
00715     
00716     solveStatus.achievedTol = TEUCHOS_MAX(solveStatus.achievedTol, achievedTol);
00717     // Note, achieveTol may actually be greater than tol due to ill conditioning and roundoff!
00718 
00719     totalIterations += iterations;
00720 
00721     solveStatus.message = oss.str();
00722     if ( isDefaultSolveCriteria ) {
00723       switch(solveStatus.solveStatus) {
00724         case SOLVE_STATUS_UNKNOWN:
00725           // Leave overall unknown!
00726           break;
00727         case SOLVE_STATUS_CONVERGED:
00728           solveStatus.solveStatus = ( converged ? SOLVE_STATUS_CONVERGED : SOLVE_STATUS_UNCONVERGED );
00729           break;
00730         case SOLVE_STATUS_UNCONVERGED:
00731           // Leave overall unconverged!
00732           break;
00733         default:
00734           TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
00735       }
00736     }
00737   }
00738 
00739   aztecSolver->UnsetLHSRHS();
00740   
00741   //
00742   // Release the Epetra_MultiVector views of X and B
00743   //
00744   epetra_X = Teuchos::null;
00745   epetra_B = Teuchos::null;
00746 
00747   //
00748   // Update the overall solve criteria
00749   //
00750   totalTimer.stop();
00751   SolveStatus<double> overallSolveStatus;
00752   if (isDefaultSolveCriteria) {
00753     overallSolveStatus.solveStatus = SOLVE_STATUS_UNKNOWN;
00754     overallSolveStatus.achievedTol = SS::unknownTolerance();
00755   }
00756   else {
00757     overallSolveStatus.solveStatus = solveStatus.solveStatus;
00758     overallSolveStatus.achievedTol = solveStatus.achievedTol;
00759   }
00760   std::ostringstream oss;
00761   oss
00762     << "AztecOO solver "
00763     << ( overallSolveStatus.solveStatus==SOLVE_STATUS_CONVERGED ? "converged" : "unconverged" )
00764     << " on m = "<<m<<" RHSs using " << totalIterations << " cumulative iterations"
00765     << " for an average of " << (totalIterations/m) << " iterations/RHS and"
00766     << " total CPU time of "<<totalTimer.totalElapsedTime()<<" sec.";
00767   overallSolveStatus.message = oss.str();
00768 
00769   // Added these statistics following what was done for Belos
00770   if (overallSolveStatus.extraParameters.is_null()) {
00771     overallSolveStatus.extraParameters = Teuchos::parameterList ();
00772   }
00773   overallSolveStatus.extraParameters->set ("AztecOO/Iteration Count",
00774                                             totalIterations);
00775   // package independent version of the same
00776   overallSolveStatus.extraParameters->set ("Iteration Count",
00777                                             totalIterations);
00778   overallSolveStatus.extraParameters->set ("AztecOO/Achieved Tolerance",
00779                                             overallSolveStatus.achievedTol);
00780 
00781   //
00782   // Report the overall time
00783   //
00784   if (out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00785     *out
00786       << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00787 
00788   return overallSolveStatus;
00789 
00790 }
00791 
00792 
00793 } // end namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines