Thyra_DampenedNewtonNonlinearSolver.hpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) 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 #ifndef THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
00030 #define THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP
00031 
00032 #include "Thyra_NonlinearSolverBase.hpp"
00033 #include "Thyra_ModelEvaluatorHelpers.hpp"
00034 #include "Thyra_TestingTools.hpp"
00035 #include "Teuchos_StandardMemberCompositionMacros.hpp"
00036 #include "Teuchos_StandardCompositionMacros.hpp"
00037 #include "Teuchos_VerboseObject.hpp"
00038 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00039 #include "Teuchos_StandardParameterEntryValidators.hpp"
00040 #include "Teuchos_as.hpp"
00041 
00042 namespace Thyra {
00043 
00062 template <class Scalar>
00063 class DampenedNewtonNonlinearSolver : public NonlinearSolverBase<Scalar> {
00064 public:
00065 
00067   typedef Teuchos::ScalarTraits<Scalar> ST;
00069   typedef typename ST::magnitudeType ScalarMag;
00071   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00072 
00074   STANDARD_MEMBER_COMPOSITION_MEMBERS( ScalarMag, defaultTol );
00075 
00077   STANDARD_MEMBER_COMPOSITION_MEMBERS( int, defaultMaxNewtonIterations );
00078 
00080   STANDARD_MEMBER_COMPOSITION_MEMBERS( bool, useDampenedLineSearch  );
00081   
00083   STANDARD_MEMBER_COMPOSITION_MEMBERS( Scalar, armijoConstant );
00084   
00086   STANDARD_MEMBER_COMPOSITION_MEMBERS( int, maxLineSearchIterations );
00087 
00089   DampenedNewtonNonlinearSolver(
00090     const ScalarMag defaultTol = 1e-2
00091     ,const int defaultMaxNewtonIterations = 1000
00092     ,const bool useDampenedLineSearch = true
00093     ,const Scalar armijoConstant = 1e-4
00094     ,const int maxLineSearchIterations = 20
00095     );
00096 
00098   static RCP<const Teuchos::ParameterList>
00099   getValidSolveCriteriaExtraParameters();
00100 
00103 
00105   void setParameterList(RCP<Teuchos::ParameterList> const& paramList);
00107   RCP<Teuchos::ParameterList> getNonconstParameterList();
00109   RCP<Teuchos::ParameterList> unsetParameterList();
00111   RCP<const Teuchos::ParameterList> getParameterList() const;
00113   RCP<const Teuchos::ParameterList> getValidParameters() const;
00114 
00116 
00119 
00121   void setModel(
00122     const RCP<const ModelEvaluator<Scalar> > &model
00123     );
00125   RCP<const ModelEvaluator<Scalar> > getModel() const;
00127   SolveStatus<Scalar> solve(
00128     VectorBase<Scalar> *x,
00129     const SolveCriteria<Scalar> *solveCriteria,
00130     VectorBase<Scalar> *delta
00131     );
00133   RCP<const VectorBase<Scalar> > get_current_x() const;
00135   bool is_W_current() const;
00137   RCP<LinearOpWithSolveBase<Scalar> > get_nonconst_W(const bool forceUpToDate);
00139   RCP<const LinearOpWithSolveBase<Scalar> > get_W() const;
00141   void set_W_is_current(bool W_is_current);
00142 
00144 
00145 private:
00146 
00147   RCP<Teuchos::ParameterList> paramList_;
00148   RCP<const ModelEvaluator<Scalar> > model_;
00149   RCP<LinearOpWithSolveBase<Scalar> > J_;
00150   RCP<VectorBase<Scalar> > current_x_;
00151   bool J_is_current_;
00152 
00153 };
00154 
00155 // ////////////////////////
00156 // Defintions
00157 
00158 template <class Scalar>
00159 DampenedNewtonNonlinearSolver<Scalar>::DampenedNewtonNonlinearSolver(
00160   const ScalarMag defaultTol
00161   ,const int defaultMaxNewtonIterations
00162   ,const bool useDampenedLineSearch
00163   ,const Scalar armijoConstant
00164   ,const int maxLineSearchIterations
00165   )
00166   :defaultTol_(defaultTol)
00167   ,defaultMaxNewtonIterations_(defaultMaxNewtonIterations)
00168   ,useDampenedLineSearch_(useDampenedLineSearch)
00169   ,armijoConstant_(armijoConstant)
00170   ,maxLineSearchIterations_(maxLineSearchIterations)
00171   ,J_is_current_(false)
00172 {}
00173 
00174 template <class Scalar>
00175 RCP<const Teuchos::ParameterList>
00176 DampenedNewtonNonlinearSolver<Scalar>::getValidSolveCriteriaExtraParameters()
00177 {
00178   static RCP<const Teuchos::ParameterList> validSolveCriteriaExtraParameters;
00179   if(!validSolveCriteriaExtraParameters.get()) {
00180     RCP<Teuchos::ParameterList>
00181       paramList = Teuchos::rcp(new Teuchos::ParameterList);
00182     paramList->set("Max Iters",int(1000));
00183     validSolveCriteriaExtraParameters = paramList;
00184   }
00185   return validSolveCriteriaExtraParameters;
00186 }
00187 
00188 // Overridden from Teuchos::ParameterListAcceptor
00189 
00190 template<class Scalar>
00191 void DampenedNewtonNonlinearSolver<Scalar>::setParameterList(
00192   RCP<Teuchos::ParameterList> const& paramList
00193   )
00194 {
00195   using Teuchos::get;
00196   TEST_FOR_EXCEPT(is_null(paramList));
00197   paramList->validateParametersAndSetDefaults(*getValidParameters(),0);
00198   paramList_ = paramList;
00199   TEST_FOR_EXCEPT("ToDo: Implement!");
00200   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00201 #ifdef TEUCHOS_DEBUG
00202   paramList_->validateParameters(*getValidParameters(),0);
00203 #endif // TEUCHOS_DEBUG
00204 }
00205 
00206 template<class Scalar>
00207 RCP<Teuchos::ParameterList>
00208 DampenedNewtonNonlinearSolver<Scalar>::getNonconstParameterList()
00209 {
00210   return paramList_;
00211 }
00212 
00213 template<class Scalar>
00214 RCP<Teuchos::ParameterList>
00215 DampenedNewtonNonlinearSolver<Scalar>::unsetParameterList()
00216 {
00217   RCP<Teuchos::ParameterList> _paramList = paramList_;
00218   paramList_ = Teuchos::null;
00219   return _paramList;
00220 }
00221 
00222 template<class Scalar>
00223 RCP<const Teuchos::ParameterList>
00224 DampenedNewtonNonlinearSolver<Scalar>::getParameterList() const
00225 {
00226   return paramList_;
00227 }
00228 
00229 template<class Scalar>
00230 RCP<const Teuchos::ParameterList>
00231 DampenedNewtonNonlinearSolver<Scalar>::getValidParameters() const
00232 {
00233   using Teuchos::setDoubleParameter; using Teuchos::setIntParameter;
00234   static RCP<const Teuchos::ParameterList> validPL;
00235   if (is_null(validPL)) {
00236     RCP<Teuchos::ParameterList>
00237       pl = Teuchos::parameterList();
00238     TEST_FOR_EXCEPT("ToDo: Implement!");
00239     Teuchos::setupVerboseObjectSublist(&*pl);
00240     validPL = pl;
00241   }
00242   return validPL;
00243 }
00244 
00245 // Overridden from NonlinearSolverBase
00246 
00247 template <class Scalar>
00248 void DampenedNewtonNonlinearSolver<Scalar>::setModel(
00249   const RCP<const ModelEvaluator<Scalar> > &model
00250   )
00251 {
00252   TEST_FOR_EXCEPT(model.get()==NULL);
00253   model_ = model;
00254   J_ = Teuchos::null;
00255   current_x_ = Teuchos::null;
00256   J_is_current_ = false;
00257 }
00258 
00259 template <class Scalar>
00260 RCP<const ModelEvaluator<Scalar> >
00261 DampenedNewtonNonlinearSolver<Scalar>::getModel() const
00262 {
00263   return model_;
00264 }
00265 
00266 template <class Scalar>
00267 SolveStatus<Scalar>
00268 DampenedNewtonNonlinearSolver<Scalar>::solve(
00269   VectorBase<Scalar> *x_inout
00270   ,const SolveCriteria<Scalar> *solveCriteria
00271   ,VectorBase<Scalar> *delta
00272   ) 
00273 {
00274 
00275   using std::endl;
00276   using Teuchos::as;
00277 
00278   // Validate input
00279 #ifdef TEUCHOS_DEBUG
00280   TEST_FOR_EXCEPT(0==x_inout);
00281   THYRA_ASSERT_VEC_SPACES(
00282     "DampenedNewtonNonlinearSolver<Scalar>::solve(...)",
00283     *x_inout->space(), *model_->get_x_space() );
00284 #endif
00285 
00286   // Get the output stream and verbosity level
00287   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00288   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00289   const bool showNewtonIters = (verbLevel==Teuchos::VERB_LOW);
00290   const bool showLineSearchIters = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_MEDIUM));
00291   const bool showNewtonDetails = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
00292   const bool dumpAll = (as<int>(verbLevel) == as<int>(Teuchos::VERB_EXTREME)); 
00293   TEUCHOS_OSTAB;
00294   if(out.get() && showNewtonIters) {
00295     *out << "\nBeginning dampended Newton solve of model = " << model_->description() << "\n\n";
00296     if (!useDampenedLineSearch())
00297       *out << "\nDoing undampened newton ...\n\n";
00298   }
00299 
00300   // Initialize storage for algorithm
00301   if(!J_.get()) J_ = model_->create_W();
00302   RCP<VectorBase<Scalar> > f = createMember(model_->get_f_space());
00303   RCP<VectorBase<Scalar> > x = Teuchos::rcp(x_inout,false);
00304   RCP<VectorBase<Scalar> > dx = createMember(model_->get_x_space());
00305   RCP<VectorBase<Scalar> > x_new = createMember(model_->get_x_space());
00306   RCP<VectorBase<Scalar> > ee = createMember(model_->get_x_space());
00307   V_S(&*ee,ST::zero());
00308 
00309   // Get convergence criteria
00310   ScalarMag tol = this->defaultTol();
00311   int maxIters = this->defaultMaxNewtonIterations();
00312   if(solveCriteria && !solveCriteria->solveMeasureType.useDefault()) {
00313     TEST_FOR_EXCEPTION(
00314       !solveCriteria->solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS), CatastrophicSolveFailure
00315       ,"DampenedNewtonNonlinearSolver<Scalar>::solve(...): Error, can only support resudual-based"
00316       " convergence criteria!");
00317     tol = solveCriteria->requestedTol;
00318     if(solveCriteria->extraParameters.get()) {
00319       solveCriteria->extraParameters->validateParameters(*getValidSolveCriteriaExtraParameters());
00320       maxIters = solveCriteria->extraParameters->get("Max Iters",int(maxIters));
00321     }
00322   }
00323 
00324   if(out.get() && showNewtonDetails)
00325     *out << "\nCompute the initial starting point ...\n";
00326 
00327   eval_f_W( *model_, *x, &*f, &*J_ );
00328   if(out.get() && dumpAll) {
00329     *out << "\nInitial starting point:\n";
00330     *out << "\nx =\n" << *x;
00331     *out << "\nf =\n" << *f;
00332     *out << "\nJ =\n" << *J_;
00333   }
00334 
00335   // Peform the Newton iterations
00336   int newtonIter, num_residual_evals = 1;
00337   SolveStatus<Scalar> solveStatus;
00338   solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
00339 
00340   for( newtonIter = 1; newtonIter <= maxIters; ++newtonIter ) {
00341 
00342     if(out.get() && showNewtonDetails) *out << "\n*** newtonIter = " << newtonIter << endl;
00343 
00344     // Check convergence
00345     if(out.get() && showNewtonDetails) *out << "\nChecking for convergence ... : ";
00346     const Scalar phi = scalarProd(*f,*f), sqrt_phi = ST::squareroot(phi); // merit function: phi(f) = <f,f>
00347     solveStatus.achievedTol = sqrt_phi;
00348     const bool isConverged = sqrt_phi <= tol;
00349     if(out.get() && showNewtonDetails) *out
00350       << "sqrt(phi) = sqrt(<f,f>) = ||f|| = " << sqrt_phi << ( isConverged ? " <= " : " > " ) << "tol = " << tol << endl;
00351     if(out.get() && showNewtonIters) *out
00352       << "newton_iter="<<newtonIter<<": Check convergence: ||f|| = "
00353       << sqrt_phi << ( isConverged ? " <= " : " > " ) << "tol = " << tol << ( isConverged ? ", Converged!!!" : "" ) << endl;
00354     if(isConverged) {
00355       if(x_inout != x.get()) assign( x_inout, *x ); // Assign the solution if we have to
00356       if(out.get() && showNewtonDetails) {
00357         *out << "\nWe have converged :-)\n"
00358              << "\n||x||inf = " << norm_inf(*x) << endl;
00359         if(dumpAll) *out << "\nx =\n" << *x;
00360         *out << "\nExiting SimpleNewtonSolver::solve(...)\n";
00361       }
00362       std::ostringstream oss;
00363       oss << "Converged! ||f|| = " << sqrt_phi << ", num_newton_iters="<<newtonIter<<", num_residual_evals="<<num_residual_evals<<".";
00364       solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00365       solveStatus.message = oss.str();
00366       break;
00367     }
00368     if(out.get() && showNewtonDetails) *out << "\nWe have to keep going :-(\n";
00369 
00370     // Compute the Jacobian if we have not already
00371     if(newtonIter > 1) {
00372       if(out.get() && showNewtonDetails) *out << "\nComputing the Jacobian J_ at current point ...\n";
00373       eval_f_W<Scalar>( *model_, *x, NULL, &*J_ );
00374       if(out.get() && dumpAll) *out << "\nJ =\n" << *J_;
00375     }
00376 
00377     // Compute the newton step: dx = -inv(J)*f
00378     if(out.get() && showNewtonDetails) *out << "\nComputing the Newton step: dx = - inv(J)*f ...\n";
00379     if(out.get() && showNewtonIters) *out << "newton_iter="<<newtonIter<<": Computing Newton step ...\n";
00380     assign( &*dx, ST::zero() ); // Initial guess for the linear solve
00381     Thyra::solve(*J_,NOTRANS,*f,&*dx); // Solve: J*dx = f
00382     Vt_S( &*dx, Scalar(-ST::one()) ); // dx *= -1.0
00383     Vp_V( &*ee, *dx); // ee += dx
00384     if(out.get() && showNewtonDetails) *out << "\n||dx||inf = " << norm_inf(*dx) << endl;
00385     if(out.get() && dumpAll) *out << "\ndy =\n" << *dx;
00386 
00387     // Perform backtracking armijo line search
00388     if(out.get() && showNewtonDetails) *out << "\nStarting backtracking line search iterations ...\n";
00389     if(out.get() && showNewtonIters) *out << "newton_iter="<<newtonIter<<": Starting backtracking line search ...\n";
00390     const Scalar Dphi = -2.0*phi; // D(phi(x+alpha*dx))/D(alpha) at alpha=0.0 => -2.0*<f,c>: where dx = -inv(J)*f
00391     Scalar alpha = 1.0; // Try a full step initially since it will eventually be accepted near solution
00392     int lineSearchIter;
00393     ++num_residual_evals;
00394     for( lineSearchIter = 1; lineSearchIter <= maxLineSearchIterations(); ++lineSearchIter, ++num_residual_evals ) {
00395       TEUCHOS_OSTAB;
00396       if(out.get() && showNewtonDetails) *out << "\n*** lineSearchIter = " << lineSearchIter << endl;
00397       // x_new = x + alpha*dx
00398       assign( &*x_new, *x ); Vp_StV( &*x_new, alpha, *dx );
00399       if(out.get() && showNewtonDetails) *out << "\n||x_new||inf = " << norm_inf(*x_new) << endl;
00400       if(out.get() && dumpAll) *out << "\nx_new =\n" << *x_new;
00401       // Compute the residual at the updated point
00402       eval_f(*model_,*x_new,&*f);
00403       if(out.get() && dumpAll) *out << "\nf_new =\n" << *f;
00404       const Scalar phi_new = scalarProd(*f,*f), phi_frac = phi + alpha * armijoConstant() * Dphi;
00405       if(out.get() && showNewtonDetails) *out << "\nphi_new = <f_new,f_new> = " << phi_new << endl;
00406       if( Teuchos::ScalarTraits<Scalar>::isnaninf(phi_new) ) {
00407         if(out.get() && showNewtonDetails) *out << "\nphi_new is not a valid number, backtracking (alpha = 0.1*alpha) ...\n";
00408         alpha *= 0.1;
00409         continue;
00410       }
00411       const bool acceptPoint = (phi_new <= phi_frac);
00412       if(out.get() && showNewtonDetails) *out
00413         << "\nphi_new = " << phi_new << ( acceptPoint ? " <= " : " > " )
00414         << "phi + alpha * eta * Dphi = " << phi << " + " << alpha << " * " << armijoConstant() << " * " << Dphi
00415         << " = " << phi_frac << endl;
00416       if(out.get() && (showLineSearchIters || (showNewtonIters && acceptPoint))) *out
00417         << "newton_iter="<<newtonIter<<", ls_iter="<<lineSearchIter<<" : "
00418         << "phi(alpha="<<alpha<<") = "<<phi_new<<(acceptPoint ? " <=" : " >")<<" armijo_cord = " << phi_frac << endl;
00419       if (out.get() && showNewtonDetails && !useDampenedLineSearch())
00420         *out << "\nUndamped newton, always accpeting the point!\n";
00421       if( acceptPoint || !useDampenedLineSearch() ) {
00422         if(out.get() && showNewtonDetails) *out << "\nAccepting the current step with step length alpha = " << alpha << "!\n";
00423         break;
00424       }
00425       if(out.get() && showNewtonDetails) *out << "\nBacktracking (alpha = 0.5*alpha) ...\n";
00426       alpha *= 0.5;
00427     }
00428 
00429     // Check for line search failure
00430     if( lineSearchIter > maxLineSearchIterations() ) {
00431       std::ostringstream oss;
00432       oss
00433         << "lineSearchIter = " << lineSearchIter << " > maxLineSearchIterations = " << maxLineSearchIterations()
00434         << ": Linear search failure! Algorithm terminated!";
00435       solveStatus.message = oss.str();
00436       if(out.get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
00437       goto exit;
00438     }
00439 
00440     // Take the Newton step
00441     std::swap<RCP<VectorBase<Scalar> > >( x_new, x ); // Now x is current point!
00442 
00443   }
00444 
00445 exit:
00446 
00447   if(out.get() && showNewtonIters) *out
00448     << "\n[Final] newton_iters="<<newtonIter<<", num_residual_evals="<<num_residual_evals<<"\n";
00449 
00450   if(newtonIter > maxIters) {
00451     std::ostringstream oss;
00452     oss
00453       << "newton_iter = " << newtonIter << " > maxIters = " << maxIters
00454       << ": Newton algorithm terminated!";
00455     solveStatus.message = oss.str();
00456     if( out.get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
00457   }
00458 
00459   if(x_inout != x.get()) assign( x_inout, *x ); // Assign the final point
00460   if(delta != NULL) assign( delta, *ee );
00461   current_x_ = x_inout->clone_v(); // Remember the final point
00462   J_is_current_ = newtonIter==1; // J is only current with x if initial point was converged!
00463 
00464   if(out.get() && showNewtonDetails) *out
00465     << "\n*** Ending dampended Newton solve." << endl; 
00466 
00467   return solveStatus;
00468 
00469 }
00470 
00471 template <class Scalar>
00472 RCP<const VectorBase<Scalar> >
00473 DampenedNewtonNonlinearSolver<Scalar>::get_current_x() const
00474 {
00475   return current_x_;
00476 }
00477 
00478 template <class Scalar>
00479 bool DampenedNewtonNonlinearSolver<Scalar>::is_W_current() const
00480 {
00481   return J_is_current_;
00482 }
00483 
00484 template <class Scalar>
00485 RCP<LinearOpWithSolveBase<Scalar> >
00486 DampenedNewtonNonlinearSolver<Scalar>::get_nonconst_W(const bool forceUpToDate)
00487 {
00488   if (forceUpToDate) {
00489     TEST_FOR_EXCEPT(forceUpToDate);
00490   }
00491   return J_;
00492 }
00493 
00494 template <class Scalar>
00495 RCP<const LinearOpWithSolveBase<Scalar> >
00496 DampenedNewtonNonlinearSolver<Scalar>::get_W() const
00497 {
00498   return J_;
00499 }
00500 
00501 template <class Scalar>
00502 void DampenedNewtonNonlinearSolver<Scalar>::set_W_is_current(bool W_is_current)
00503 {
00504   J_is_current_ = W_is_current;
00505 }
00506 
00507 } // namespace Thyra
00508 
00509 #endif // THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP

Generated on Wed May 12 21:27:21 2010 for Thyra Nonlinear Solver Support by  doxygen 1.4.7