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       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     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
00183   ,const double fwdDefaultTol
00184   ,const int adjDefaultMaxIterations
00185   ,const double adjDefaultTol
00186   ,const bool outputEveryRhs
00187   )
00188   :fwdDefaultMaxIterations_(fwdDefaultMaxIterations)
00189   ,fwdDefaultTol_(fwdDefaultTol)
00190   ,adjDefaultMaxIterations_(adjDefaultMaxIterations)
00191   ,adjDefaultTol_(adjDefaultTol)
00192   ,outputEveryRhs_(outputEveryRhs)
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
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   TEST_FOR_EXCEPT(fwdOp.get()==NULL);
00215   TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00216   TEST_FOR_EXCEPT(aztecFwdSolver.get()==NULL);
00217 #endif
00218   fwdOp_ = fwdOp;
00219   fwdOpSrc_ = fwdOpSrc;
00220   isExternalPrec_ = isExternalPrec;
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,
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) *isExternalPrec = 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       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   TEUCHOS_FUNC_TIME_MONITOR("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   }
00597 
00598   //
00599   // Get Epetra_MultiVector views of B and X
00600   //
00601 
00602   RCP<const Epetra_MultiVector> epetra_B;
00603   RCP<Epetra_MultiVector> epetra_X;
00604 
00605   const EpetraOperatorWrapper* opWrapper =
00606     dynamic_cast<const EpetraOperatorWrapper*>(aztecOp);
00607 
00608   if (opWrapper == 0) {
00609     epetra_B = get_Epetra_MultiVector(opRangeMap, rcpFromRef(B));
00610     epetra_X = get_Epetra_MultiVector(opDomainMap, rcpFromPtr(X));
00611   }
00612 
00613   //
00614   // Use AztecOO to solve each RHS one at a time (which is all that I can do anyway)
00615   //
00616 
00617   int totalIterations = 0;
00618   SolveStatus<double> solveStatus;
00619   solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00620   solveStatus.achievedTol = -1.0;
00621 
00622   /* Get the number of columns in the multivector. We use Thyra
00623    * functions rather than Epetra functions to do this, as we
00624    * might not yet have created an Epetra multivector. - KL */
00625   //const int m = epetra_B->NumVectors();
00626   const int m = B.domain()->dim();
00627 
00628   for( int j = 0; j < m; ++j ) {
00629 
00630     TEUCHOS_FUNC_TIME_MONITOR_DIFF("AztecOOLOWS:SingleSolve", SingleSolve);
00631 
00632     //
00633     // Get Epetra_Vector views of B(:,j) and X(:,j)
00634     // How this is done will depend on whether we have a true Epetra operator
00635     // or we are wrapping a general Thyra operator in an Epetra operator.
00636     //
00637 
00638     // We need to declare epetra_x_j as non-const because when we have a phony
00639     // Epetra operator we'll have to copy a thyra vector into it.
00640     RCP<Epetra_Vector> epetra_b_j;
00641     RCP<Epetra_Vector> epetra_x_j;
00642 
00643     if (opWrapper == 0) {
00644       epetra_b_j = rcpFromRef(*const_cast<Epetra_Vector*>((*epetra_B)(j)));
00645       epetra_x_j = rcpFromRef(*(*epetra_X)(j));
00646     }
00647     else {
00648       if (is_null(epetra_b_j)) {
00649         epetra_b_j = rcp(new Epetra_Vector(opRangeMap));
00650         epetra_x_j = rcp(new Epetra_Vector(opDomainMap));
00651       }
00652       opWrapper->copyThyraIntoEpetra(*B.col(j), *epetra_b_j);
00653       opWrapper->copyThyraIntoEpetra(*X->col(j), *epetra_x_j);
00654     }
00655 
00656     //
00657     // Set the RHS and LHS
00658     //
00659 
00660     aztecSolver->SetRHS(&*epetra_b_j);
00661     aztecSolver->SetLHS(&*epetra_x_j);
00662 
00663     //
00664     // Solve the linear system
00665     //
00666     timer.start(true);
00667     {
00668       SetAztecSolveState
00669         setAztecSolveState(aztecSolver,out,verbLevel,solveMeasureType);
00670       aztecSolver->Iterate( maxIterations, tol );
00671       // NOTE: We ignore the returned status but get it below
00672     }
00673     timer.stop();
00674 
00675     //
00676     // Scale the solution 
00677     // (Originally, this was at the end of the loop after all columns had been
00678     // processed. It's moved here because we need to do it before copying the
00679     // solution back into a Thyra vector. - KL
00680     //
00681     if (aztecSolverScalar_ != 1.0)
00682       epetra_x_j->Scale(1.0/aztecSolverScalar_);
00683 
00684     //
00685     // If necessary, convert the solution back to a non-epetra vector
00686     //
00687     if (opWrapper != 0) {
00688       opWrapper->copyEpetraIntoThyra(*epetra_x_j, X->col(j).ptr());
00689     }
00690 
00691     //
00692     // Set the return solve status
00693     //
00694 
00695     const int iterations = aztecSolver->NumIters();
00696     const double achievedTol = aztecSolver->ScaledResidual();
00697     const double *AZ_status = aztecSolver->GetAztecStatus();
00698     std::ostringstream oss;
00699     bool converged = false;
00700     if (AZ_status[AZ_why]==AZ_normal) { oss << "Aztec returned AZ_normal."; converged = true; }
00701     else if (AZ_status[AZ_why]==AZ_param) oss << "Aztec returned AZ_param.";
00702     else if (AZ_status[AZ_why]==AZ_breakdown) oss << "Aztec returned AZ_breakdown.";
00703     else if (AZ_status[AZ_why]==AZ_loss) oss << "Aztec returned AZ_loss.";
00704     else if (AZ_status[AZ_why]==AZ_ill_cond) oss << "Aztec returned AZ_ill_cond.";
00705     else if (AZ_status[AZ_why]==AZ_maxits) oss << "Aztec returned AZ_maxits.";
00706     else oss << "Aztec returned an unknown status?";
00707     oss << "  Iterations = " << iterations << ".";
00708     oss << "  Achieved Tolerance = " << achievedTol << ".";
00709     oss << "  Total time = " << timer.totalElapsedTime() << " sec.";
00710     if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE) && outputEveryRhs())
00711       Teuchos::OSTab(out).o() << "j="<<j<<": " << oss.str() << "\n";
00712     
00713     solveStatus.achievedTol = TEUCHOS_MAX(solveStatus.achievedTol, achievedTol);
00714     // Note, achieveTol may actually be greater than tol due to ill conditioning and roundoff!
00715 
00716     solveStatus.message = oss.str();
00717     if ( isDefaultSolveCriteria ) {
00718       switch(solveStatus.solveStatus) {
00719         case SOLVE_STATUS_UNKNOWN:
00720           // Leave overall unknown!
00721           break;
00722         case SOLVE_STATUS_CONVERGED:
00723           solveStatus.solveStatus = ( converged ? SOLVE_STATUS_CONVERGED : SOLVE_STATUS_UNCONVERGED );
00724           break;
00725         case SOLVE_STATUS_UNCONVERGED:
00726           // Leave overall unconverged!
00727           break;
00728         default:
00729           TEST_FOR_EXCEPT(true); // Should never get here!
00730       }
00731     }
00732   }
00733 
00734   aztecSolver->UnsetLHSRHS();
00735   
00736   //
00737   // Release the Epetra_MultiVector views of X and B
00738   //
00739   epetra_X = Teuchos::null;
00740   epetra_B = Teuchos::null;
00741 
00742   //
00743   // Update the overall solve criteria
00744   //
00745   totalTimer.stop();
00746   SolveStatus<double> overallSolveStatus;
00747   std::ostringstream oss;
00748   oss
00749     << "AztecOO solver "
00750     << ( overallSolveStatus.solveStatus==SOLVE_STATUS_CONVERGED ? "converged" : "unconverged" )
00751     << " on m = "<<m<<" RHSs using " << totalIterations << " cumulative iterations"
00752     << " for an average of " << (totalIterations/m) << " iterations/RHS and"
00753     << " total CPU time of "<<totalTimer.totalElapsedTime()<<" sec.";
00754   overallSolveStatus.message = oss.str();
00755   if (isDefaultSolveCriteria) {
00756     overallSolveStatus.solveStatus = SOLVE_STATUS_UNKNOWN;
00757     overallSolveStatus.achievedTol = SS::unknownTolerance();
00758   }
00759   else {
00760     overallSolveStatus.solveStatus = solveStatus.solveStatus;
00761     overallSolveStatus.achievedTol = solveStatus.achievedTol;
00762   }
00763 
00764   //
00765   // Report the overall time
00766   //
00767   if (out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00768     *out
00769       << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00770 
00771   return overallSolveStatus;
00772 
00773 }
00774 
00775 
00776 } // end namespace Thyra
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends
Generated on Wed Apr 13 10:19:51 2011 for Stratimikos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3