Thyra_AztecOOLinearOpWithSolve.cpp

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

Generated on Wed Jul 22 12:56:30 2009 for Stratimikos by  doxygen 1.5.8