Thyra_BelosLinearOpWithSolve.hpp

00001 
00002 #ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
00003 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
00004 
00005 #include "Thyra_BelosLinearOpWithSolveDecl.hpp"
00006 #include "Thyra_LinearOpWithSolveHelpers.hpp"
00007 #include "Teuchos_TimeMonitor.hpp"
00008 
00009 namespace Thyra {
00010 
00011 // Constructors/initializers/accessors
00012 
00013 template<class Scalar>
00014 BelosLinearOpWithSolve<Scalar>::BelosLinearOpWithSolve()
00015   :isExternalPrec_(false)
00016   ,supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED)
00017   ,defaultTol_(-1.0)
00018 {}
00019 
00020 template<class Scalar>
00021 BelosLinearOpWithSolve<Scalar>::BelosLinearOpWithSolve(
00022   const Teuchos::RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> >         &lp
00023   ,const Teuchos::RCP<Teuchos::ParameterList>                         &solverPL
00024   ,const Teuchos::RCP<Belos::SolverManager<Scalar,MV_t,LO_t> >        &iterativeSolver
00025   ,const Teuchos::RCP<const LinearOpSourceBase<Scalar> >              &fwdOpSrc
00026   ,const Teuchos::RCP<const PreconditionerBase<Scalar> >              &prec
00027   ,const bool                                                         isExternalPrec
00028   ,const Teuchos::RCP<const LinearOpSourceBase<Scalar> >              &approxFwdOpSrc
00029   ,const ESupportSolveUse                                             &supportSolveUse
00030   )
00031 {
00032   initialize(
00033     lp,solverPL,iterativeSolver
00034     ,fwdOpSrc,prec,isExternalPrec,approxFwdOpSrc,supportSolveUse
00035     );
00036 }
00037 
00038 template<class Scalar>
00039 void BelosLinearOpWithSolve<Scalar>::initialize(
00040   const Teuchos::RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> >         &lp
00041   ,const Teuchos::RCP<Teuchos::ParameterList>                         &solverPL
00042   ,const Teuchos::RCP<Belos::SolverManager<Scalar,MV_t,LO_t> >        &iterativeSolver
00043   ,const Teuchos::RCP<const LinearOpSourceBase<Scalar> >              &fwdOpSrc
00044   ,const Teuchos::RCP<const PreconditionerBase<Scalar> >              &prec
00045   ,const bool                                                         isExternalPrec
00046   ,const Teuchos::RCP<const LinearOpSourceBase<Scalar> >              &approxFwdOpSrc
00047   ,const ESupportSolveUse                                             &supportSolveUse
00048   )
00049 {
00050   this->setLinePrefix("BELOS/T");
00051   // ToDo: Validate input
00052   lp_ = lp;
00053   solverPL_ = solverPL;
00054   iterativeSolver_ = iterativeSolver;
00055   fwdOpSrc_ = fwdOpSrc;
00056   prec_ = prec;
00057   isExternalPrec_ = isExternalPrec;
00058   approxFwdOpSrc_ = approxFwdOpSrc;
00059   supportSolveUse_ = supportSolveUse;
00060   // Check if "Convergence Tolerance" is in the solver parameter list.  If not, use the default from the solver.
00061   if ( !is_null(solverPL_) ) {
00062     if (solverPL_->isParameter("Convergence Tolerance")) {
00063       defaultTol_ = solverPL_->get<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>("Convergence Tolerance");
00064     }
00065   }
00066   else {
00067     Teuchos::RCP<const Teuchos::ParameterList> defaultPL = iterativeSolver->getValidParameters();
00068     defaultTol_ = defaultPL->get<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>("Convergence Tolerance");
00069   }
00070 }
00071 
00072 template<class Scalar>
00073 Teuchos::RCP<const LinearOpSourceBase<Scalar> >
00074 BelosLinearOpWithSolve<Scalar>::extract_fwdOpSrc()
00075 {
00076   Teuchos::RCP<const LinearOpSourceBase<Scalar> >
00077     _fwdOpSrc = fwdOpSrc_;
00078   fwdOpSrc_ = Teuchos::null;
00079   return _fwdOpSrc;
00080 }
00081 
00082 template<class Scalar>
00083 Teuchos::RCP<const PreconditionerBase<Scalar> >
00084 BelosLinearOpWithSolve<Scalar>::extract_prec()
00085 {
00086   Teuchos::RCP<const PreconditionerBase<Scalar> >
00087     _prec = prec_;
00088   prec_ = Teuchos::null;
00089   return _prec;
00090 }
00091 
00092 template<class Scalar>
00093 bool BelosLinearOpWithSolve<Scalar>::isExternalPrec() const
00094 {
00095   return isExternalPrec_;
00096 }
00097 
00098 template<class Scalar>
00099 Teuchos::RCP<const LinearOpSourceBase<Scalar> >
00100 BelosLinearOpWithSolve<Scalar>::extract_approxFwdOpSrc()
00101 {
00102   Teuchos::RCP<const LinearOpSourceBase<Scalar> >
00103     _approxFwdOpSrc = approxFwdOpSrc_;
00104   approxFwdOpSrc_ = Teuchos::null;
00105   return _approxFwdOpSrc;
00106 }
00107 
00108 template<class Scalar>
00109 ESupportSolveUse BelosLinearOpWithSolve<Scalar>::supportSolveUse() const
00110 {
00111   return supportSolveUse_;
00112 }
00113 
00114 template<class Scalar>
00115 void BelosLinearOpWithSolve<Scalar>::uninitialize(
00116   Teuchos::RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> >         *lp
00117   ,Teuchos::RCP<Teuchos::ParameterList>                         *solverPL
00118   ,Teuchos::RCP<Belos::SolverManager<Scalar,MV_t,LO_t> >        *iterativeSolver
00119   ,Teuchos::RCP<const LinearOpSourceBase<Scalar> >              *fwdOpSrc
00120   ,Teuchos::RCP<const PreconditionerBase<Scalar> >              *prec
00121   ,bool                                                         *isExternalPrec
00122   ,Teuchos::RCP<const LinearOpSourceBase<Scalar> >              *approxFwdOpSrc
00123   ,ESupportSolveUse                                             *supportSolveUse
00124   )
00125 {
00126   if(lp) *lp = lp_;
00127   if(solverPL) *solverPL = solverPL_;
00128   if(iterativeSolver) *iterativeSolver = iterativeSolver_;
00129   if(fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
00130   if(prec) *prec = prec_;
00131   if(isExternalPrec) *isExternalPrec = isExternalPrec_;
00132   if(approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
00133   if(supportSolveUse) *supportSolveUse = supportSolveUse_;
00134 
00135   lp_ = Teuchos::null;
00136   solverPL_ = Teuchos::null;
00137   iterativeSolver_ = Teuchos::null;
00138   fwdOpSrc_ = Teuchos::null;
00139   prec_ = Teuchos::null;
00140   isExternalPrec_ = false;
00141   approxFwdOpSrc_ = Teuchos::null;
00142   supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
00143 }
00144 
00145 // Overridden from LinearOpBase
00146 
00147 template<class Scalar>
00148 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00149 BelosLinearOpWithSolve<Scalar>::range() const
00150 {
00151   if (!is_null(lp_))
00152     return lp_->getOperator()->range();
00153   return Teuchos::null;
00154 }
00155 
00156 template<class Scalar>
00157 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00158 BelosLinearOpWithSolve<Scalar>::domain() const
00159 {
00160   if (!is_null(lp_))
00161     return lp_->getOperator()->domain();
00162   return Teuchos::null;
00163 }
00164 
00165 template<class Scalar>
00166 Teuchos::RCP<const LinearOpBase<Scalar> >
00167 BelosLinearOpWithSolve<Scalar>::clone() const
00168 {
00169   return Teuchos::null; // Not supported yet but could be
00170 }
00171 
00172 // Overridden from Teuchos::Describable
00173 
00174 template<class Scalar>
00175 std::string BelosLinearOpWithSolve<Scalar>::description() const
00176 {
00177   std::ostringstream oss;
00178   oss << Teuchos::Describable::description();
00179   if( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
00180     oss << "{";
00181     oss << "iterativeSolver=\'"<<iterativeSolver_->description()<<"\'";
00182     oss << ",fwdOp=\'"<<lp_->getOperator()->description()<<"\'";
00183     if(lp_->getLeftPrec().get())
00184       oss << ",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<"\'";
00185     if(lp_->getRightPrec().get())
00186       oss << ",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<"\'";
00187     oss << "}";
00188   }
00189   // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
00190   // that we can get better information.
00191   return oss.str();
00192 }
00193 
00194 template<class Scalar>
00195 void BelosLinearOpWithSolve<Scalar>::describe(
00196   Teuchos::FancyOStream                &out_arg
00197   ,const Teuchos::EVerbosityLevel      verbLevel
00198   ) const
00199 {
00200   typedef Teuchos::ScalarTraits<Scalar>  ST;
00201   using Teuchos::RCP;
00202   using Teuchos::FancyOStream;
00203   using Teuchos::OSTab;
00204   using Teuchos::describe;
00205   RCP<FancyOStream> out = rcp(&out_arg,false);
00206   OSTab tab(out);
00207   switch(verbLevel) {
00208     case Teuchos::VERB_DEFAULT:
00209     case Teuchos::VERB_LOW:
00210       *out << this->description() << std::endl;
00211       break;
00212     case Teuchos::VERB_MEDIUM:
00213     case Teuchos::VERB_HIGH:
00214     case Teuchos::VERB_EXTREME:
00215     {
00216       *out
00217         << Teuchos::Describable::description()<< "{"
00218         << "rangeDim=" << this->range()->dim()
00219         << ",domainDim=" << this->domain()->dim() << "}\n";
00220       if(lp_->getOperator().get()) {
00221         OSTab tab(out);
00222         *out
00223           << "iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
00224           << "fwdOp = " << describe(*lp_->getOperator(),verbLevel);
00225         if(lp_->getLeftPrec().get())
00226           *out << "leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
00227         if(lp_->getRightPrec().get())
00228           *out << "rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
00229       }
00230       break;
00231     }
00232     default:
00233       TEST_FOR_EXCEPT(true); // Should never get here!
00234   }
00235 }
00236 
00237 // protected
00238 
00239 // Overridden from SingleScalarLinearOpBase
00240 
00241 template<class Scalar>
00242 bool BelosLinearOpWithSolve<Scalar>::opSupported(EOpTransp M_trans) const
00243 {
00244   return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
00245 }
00246 
00247 template<class Scalar>
00248 void BelosLinearOpWithSolve<Scalar>::apply(
00249   const EOpTransp                     M_trans
00250   ,const MultiVectorBase<Scalar>    &X
00251   ,MultiVectorBase<Scalar>          *Y
00252   ,const Scalar                     alpha
00253   ,const Scalar                     beta
00254   ) const
00255 {
00256   ::Thyra::apply(*lp_->getOperator(),M_trans,X,Y,alpha,beta);
00257 }
00258 
00259 // Overridden from SingleScalarLinearOpWithSolveBase
00260 
00261 template<class Scalar>
00262 bool BelosLinearOpWithSolve<Scalar>::solveSupportsTrans(EOpTransp M_trans) const
00263 {
00264   if(real_trans(M_trans)==NOTRANS) return true;
00265   return false; // ToDo: Support adjoint solves!
00266 }
00267 
00268 template<class Scalar>
00269 bool BelosLinearOpWithSolve<Scalar>::solveSupportsSolveMeasureType(EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
00270 {
00271   if(real_trans(M_trans)==NOTRANS) {
00272     return (
00273       solveMeasureType.useDefault()
00274       ||
00275       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
00276       ||
00277       solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
00278       );
00279   }
00280   // TRANS
00281   return false; // ToDo: Support adjoint solves!
00282 }
00283 
00284 template<class Scalar>
00285 void BelosLinearOpWithSolve<Scalar>::solve(
00286   const EOpTransp                         M_trans
00287   ,const MultiVectorBase<Scalar>        &B
00288   ,MultiVectorBase<Scalar>              *X
00289   ,const int                            numBlocks
00290   ,const BlockSolveCriteria<Scalar>     blockSolveCriteria[]
00291   ,SolveStatus<Scalar>                  blockSolveStatus[]
00292   ) const
00293 {
00294 
00295   TEUCHOS_FUNC_TIME_MONITOR("BelosLOWS");
00296 
00297   using Teuchos::RCP;
00298   using Teuchos::rcp;
00299   using Teuchos::FancyOStream;
00300   using Teuchos::OSTab;
00301   using Teuchos::describe;
00302   typedef Teuchos::ScalarTraits<Scalar> ST;
00303   typedef typename ST::magnitudeType ScalarMag;
00304   Teuchos::Time totalTimer(""), timer("");
00305   totalTimer.start(true);
00306   
00307   TEST_FOR_EXCEPT(numBlocks > 1); // ToDo: Deal with multiple solve criteria later if needed
00308   TEST_FOR_EXCEPT(!this->solveSupportsTrans(M_trans)); // ToDo: Support adjoint solves later!
00309 
00310   const int numRhs = B.domain()->dim();
00311   const int numEquations = B.range()->dim();
00312 
00313   Teuchos::RCP<Teuchos::FancyOStream>  out = this->getOStream();
00314   Teuchos::EVerbosityLevel                     verbLevel = this->getVerbLevel();
00315   OSTab tab = this->getOSTab();
00316   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE)) {
00317     *out << "\nStarting iterations with Belos:\n";
00318     OSTab tab(out);
00319     *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
00320     *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
00321     *out << "With #Eqns="<<numEquations<<", #RHSs="<<numRhs<<" ...\n";
00322   }
00323   //
00324   // Set RHS and LHS
00325   //
00326   bool ret = lp_->setProblem(Teuchos::rcp(X,false),Teuchos::rcp(&B,false));
00327   TEST_FOR_EXCEPTION(
00328     ret == false, CatastrophicSolveFailure
00329     ,"Error, the Belos::LinearProblem could not be set for the current solve!"
00330     );
00331   //
00332   // Set the solution criteria
00333   //
00334   Teuchos::RCP<Teuchos::ParameterList> tmpPL = Teuchos::rcp( new Teuchos::ParameterList() );
00335 
00336   SolveMeasureType solveMeasureType;
00337   if(numBlocks==1) {
00338     solveMeasureType = blockSolveCriteria[0].solveCriteria.solveMeasureType;
00339     const ScalarMag requestedTol = blockSolveCriteria[0].solveCriteria.requestedTol;
00340     assertSupportsSolveMeasureType(*this,M_trans,solveMeasureType);
00341     if( solveMeasureType.useDefault() ) {
00342       tmpPL->set("Convergence Tolerance", defaultTol_);
00343     }
00344     else if( solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS) ) {
00345       if( requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance() )
00346         tmpPL->set("Convergence Tolerance", requestedTol);
00347       else
00348         tmpPL->set("Convergence Tolerance", defaultTol_);
00349       tmpPL->set("Explicit Residual Scaling", "Norm of RHS");
00350     }
00351     else if( solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL) ) {
00352       if( requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance() )
00353         tmpPL->set("Convergence Tolerance", requestedTol);
00354       else
00355         tmpPL->set("Convergence Tolerance", defaultTol_);
00356       tmpPL->set("Explicit Residual Scaling", "Norm of Initial Residual");
00357     }
00358     else {
00359       TEST_FOR_EXCEPT(true); // Should never get there.
00360     }
00361   }
00362   else {
00363     tmpPL->set("Convergence Tolerance", defaultTol_);
00364   }
00365   //
00366   // Reset the blocksize if we adding more vectors than half the number of equations,
00367   // orthogonalization will fail on the first iteration!
00368   //
00369   Teuchos::RCP<const Teuchos::ParameterList> solverParams = iterativeSolver_->getCurrentParameters();
00370   const int currBlockSize = Teuchos::getParameter<int>(*solverParams, "Block Size");
00371   bool isNumBlocks = false;
00372   int currNumBlocks = 0;
00373   if (Teuchos::isParameterType<int>(*solverParams, "Num Blocks")) {
00374     currNumBlocks = Teuchos::getParameter<int>(*solverParams, "Num Blocks");
00375     isNumBlocks = true;
00376   }
00377   const int newBlockSize = TEUCHOS_MIN(currBlockSize,numEquations/2);
00378   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE) && newBlockSize != currBlockSize)
00379     *out << "\nAdjusted block size = " << newBlockSize << "\n";
00380   //
00381   tmpPL->set("Block Size",newBlockSize);
00382   //
00383   // Set the number of Krylov blocks if we are using a GMRES solver, or a solver
00384   // that recognizes "Num Blocks". Otherwise the solver will throw an error!
00385   if (isNumBlocks) {
00386     const int Krylov_length = (currNumBlocks*currBlockSize)/newBlockSize;
00387     tmpPL->set("Num Blocks",Krylov_length);
00388   
00389     if(newBlockSize != currBlockSize) {
00390       if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00391         *out
00392           << "\nAdjusted max number of Krylov basis blocks = " << Krylov_length << "\n";
00393     }
00394   }
00395   //
00396   // Solve the linear system
00397   //
00398   Belos::ReturnType belosSolveStatus;
00399   {
00400     RCP<std::ostream>
00401       outUsed =
00402       ( static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE)
00403         ? out
00404         : rcp(new FancyOStream(rcp(new Teuchos::oblackholestream())))
00405         );
00406     Teuchos::OSTab tab(outUsed,1,"BELOS");
00407     tmpPL->set("Output Stream", outUsed);
00408     iterativeSolver_->setParameters( tmpPL );
00409     belosSolveStatus = iterativeSolver_->solve();
00410   }
00411   //
00412   // Report the solve status
00413   //
00414   totalTimer.stop();
00415   //
00416   ESolveStatus solveStatus = SOLVE_STATUS_UNKNOWN;
00417   ScalarMag    achievedTol = SolveStatus<Scalar>::unknownTolerance();
00418   switch(belosSolveStatus) {
00419     case Belos::Unconverged:
00420       solveStatus = SOLVE_STATUS_UNCONVERGED;
00421       break;
00422     case Belos::Converged:
00423       solveStatus = SOLVE_STATUS_CONVERGED;
00424       achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
00425       break;
00426     default:
00427       TEST_FOR_EXCEPT(true); // Should never get here!
00428   }
00429   //const int iterations = iterativeSolver_->getNumIters();
00430 
00431   std::ostringstream ossmessage;
00432   ossmessage
00433     << "The Belos solver of type \""<<iterativeSolver_->description()<<"\" returned a solve status of \""
00434     << toString(solveStatus) << 
00435     /* "\" in " << iterations << " iterations and achieved an approximate max tolerance of "
00436        << belosAchievedTol << 
00437     */
00438     " with total CPU time of " << totalTimer.totalElapsedTime() << " sec" ;
00439   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
00440     *out << "\n" << ossmessage.str() << "\n";
00441   //
00442   if( numBlocks && blockSolveStatus ) {
00443     blockSolveStatus[0].solveStatus = solveStatus;
00444     blockSolveStatus[0].achievedTol = achievedTol;
00445     blockSolveStatus[0].message     = ossmessage.str();
00446   }
00447 
00448   if(out.get() && static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW))
00449     *out
00450       << "\nTotal solve time = "<<totalTimer.totalElapsedTime()<<" sec\n";
00451 }
00452 
00453 } // end namespace Thyra
00454 
00455 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP

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