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