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 
00039 namespace Thyra {
00040 
00059 template <class Scalar>
00060 class DampenedNewtonNonlinearSolver : public NonlinearSolverBase<Scalar> {
00061 public:
00062 
00064   typedef Teuchos::ScalarTraits<Scalar> ST;
00066   typedef typename ST::magnitudeType ScalarMag;
00068   typedef Teuchos::ScalarTraits<ScalarMag> SMT;
00069 
00071    STANDARD_MEMBER_COMPOSITION_MEMBERS( ScalarMag, defaultTol )
00072 
00073   
00074   STANDARD_MEMBER_COMPOSITION_MEMBERS( int, defaultMaxNewtonIterations )
00075   
00077   STANDARD_MEMBER_COMPOSITION_MEMBERS( Scalar, armijoConstant )
00078   
00080   STANDARD_MEMBER_COMPOSITION_MEMBERS( int, maxLineSearchIterations )
00081 
00083   DampenedNewtonNonlinearSolver(
00084     const ScalarMag           defaultTol                   = 1e-2
00085     ,const int                defaultMaxNewtonIterations   = 1000
00086     ,const Scalar             armijoConstant               = 1e-4
00087     ,const int                maxLineSearchIterations      = 20
00088     );
00089 
00091   static Teuchos::RefCountPtr<const Teuchos::ParameterList> getValidSolveCriteriaExtraParameters();
00092 
00095 
00097   void setModel(
00098     const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > &model
00099     );
00101   Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > getModel() const;
00103   SolveStatus<Scalar> solve(
00104     VectorBase<Scalar>              *x
00105     ,const SolveCriteria<Scalar>    *solveCriteria
00106     ,VectorBase<Scalar>             *delta = NULL
00107     );
00109   Teuchos::RefCountPtr<const VectorBase<Scalar> > get_current_x() const;
00111   bool is_W_current() const;
00113   Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> > get_nonconst_W();
00115   Teuchos::RefCountPtr<const LinearOpWithSolveBase<Scalar> > get_W() const;
00117   void set_W_is_current(bool W_is_current);
00118 
00120 
00121 private:
00122 
00123   Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >    model_;
00124   Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> >   J_;
00125   Teuchos::RefCountPtr<VectorBase<Scalar> >              current_x_;
00126   bool                                                   J_is_current_;
00127 
00128 };
00129 
00130 // ////////////////////////
00131 // Defintions
00132 
00133 template <class Scalar>
00134 DampenedNewtonNonlinearSolver<Scalar>::DampenedNewtonNonlinearSolver(
00135   const ScalarMag           defaultTol
00136   ,const int                defaultMaxNewtonIterations
00137   ,const Scalar             armijoConstant
00138   ,const int                maxLineSearchIterations
00139   )
00140   :defaultTol_(defaultTol)
00141   ,defaultMaxNewtonIterations_(defaultMaxNewtonIterations)
00142   ,armijoConstant_(armijoConstant)
00143   ,maxLineSearchIterations_(maxLineSearchIterations)
00144   ,J_is_current_(false)
00145 {}
00146 
00147 template <class Scalar>
00148 Teuchos::RefCountPtr<const Teuchos::ParameterList>
00149 DampenedNewtonNonlinearSolver<Scalar>::getValidSolveCriteriaExtraParameters()
00150 {
00151   static Teuchos::RefCountPtr<const Teuchos::ParameterList> validSolveCriteriaExtraParameters;
00152   if(!validSolveCriteriaExtraParameters.get()) {
00153     Teuchos::RefCountPtr<Teuchos::ParameterList>
00154       paramList = Teuchos::rcp(new Teuchos::ParameterList);
00155     paramList->set("Max Iters",int(1000));
00156     validSolveCriteriaExtraParameters = paramList;
00157   }
00158   return validSolveCriteriaExtraParameters;
00159 }
00160 
00161 // Overridden from NonlinearSolverBase
00162 
00163 template <class Scalar>
00164 void DampenedNewtonNonlinearSolver<Scalar>::setModel(
00165   const Teuchos::RefCountPtr<const ModelEvaluator<Scalar> > &model
00166   )
00167 {
00168   TEST_FOR_EXCEPT(model.get()==NULL);
00169   model_ = model;
00170   J_ = Teuchos::null;
00171   current_x_ = Teuchos::null;
00172   J_is_current_ = false;
00173 }
00174 
00175 template <class Scalar>
00176 Teuchos::RefCountPtr<const ModelEvaluator<Scalar> >
00177 DampenedNewtonNonlinearSolver<Scalar>::getModel() const
00178 {
00179   return model_;
00180 }
00181 
00182 template <class Scalar>
00183 SolveStatus<Scalar>
00184 DampenedNewtonNonlinearSolver<Scalar>::solve(
00185   VectorBase<Scalar>                    *x_inout
00186   ,const SolveCriteria<Scalar>          *solveCriteria
00187   ,VectorBase<Scalar>                   *delta
00188   ) 
00189 {
00190   using std::endl;
00191   // Validate input
00192   THYRA_ASSERT_VEC_SPACES(
00193     "DampenedNewtonNonlinearSolver<Scalar>::solve(...)",*x_inout->space(),*model_->get_x_space());
00194   // Get the output stream and verbosity level
00195   const Teuchos::RefCountPtr<Teuchos::FancyOStream> out = this->getOStream();
00196   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00197   const bool showNewtonIters = (verbLevel==Teuchos::VERB_LOW);
00198   const bool showLineSearchIters = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
00199   const bool showNewtonDetails = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_HIGH));
00200   const bool dumpAll = (static_cast<int>(verbLevel) == static_cast<int>(Teuchos::VERB_EXTREME)); 
00201   TEUCHOS_OSTAB;
00202   if(out.get() && showNewtonIters) *out
00203     << "\nBeginning dampended Newton solve of model = " << model_->description() << "\n\n";
00204   // Initialize storage for algorithm
00205   if(!J_.get())                                        J_    = model_->create_W();
00206   Teuchos::RefCountPtr<VectorBase<Scalar> >            f     = createMember(model_->get_f_space());
00207   Teuchos::RefCountPtr<VectorBase<Scalar> >            x     = Teuchos::rcp(x_inout,false);
00208   Teuchos::RefCountPtr<VectorBase<Scalar> >            dx    = createMember(model_->get_x_space());
00209   Teuchos::RefCountPtr<VectorBase<Scalar> >            x_new = createMember(model_->get_x_space());
00210   Teuchos::RefCountPtr<VectorBase<Scalar> >            ee    = createMember(model_->get_x_space());
00211   V_S(&*ee,ST::zero());
00212   // Get convergence criteria
00213   ScalarMag tol = this->defaultTol();
00214   int maxIters = this->defaultMaxNewtonIterations();
00215   if(solveCriteria && !solveCriteria->solveMeasureType.useDefault()) {
00216     TEST_FOR_EXCEPTION(
00217       !solveCriteria->solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS), CatastrophicSolveFailure
00218       ,"DampenedNewtonNonlinearSolver<Scalar>::solve(...): Error, can only support resudual-based"
00219       " convergence criteria!");
00220     tol = solveCriteria->requestedTol;
00221     if(solveCriteria->extraParameters.get()) {
00222       solveCriteria->extraParameters->validateParameters(*getValidSolveCriteriaExtraParameters());
00223       maxIters = solveCriteria->extraParameters->get("Max Iters",int(maxIters));
00224     }
00225   }
00226   if(out.get() && showNewtonDetails)
00227     *out << "\nCompute the initial starting point ...\n";
00228   eval_f_W( *model_, *x, &*f, &*J_ );
00229   if(out.get() && dumpAll) {
00230     *out << "\nInitial starting point:\n";
00231     *out << "\nx =\n" << *x;
00232     *out << "\nf =\n" << *f;
00233     *out << "\nJ =\n" << *J_;
00234   }
00235   // Peform the Newton iterations
00236   int newtonIter, num_residual_evals = 1;
00237   SolveStatus<Scalar> solveStatus;
00238   solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
00239   for( newtonIter = 1; newtonIter <= maxIters; ++newtonIter ) {
00240     if(out.get() && showNewtonDetails) *out << "\n*** newtonIter = " << newtonIter << endl;
00241     // Check convergence
00242     if(out.get() && showNewtonDetails) *out << "\nChecking for convergence ... : ";
00243     const Scalar phi = scalarProd(*f,*f), sqrt_phi = ST::squareroot(phi); // merit function: phi(f) = <f,f>
00244     solveStatus.achievedTol = sqrt_phi;
00245     const bool isConverged = sqrt_phi <= tol;
00246     if(out.get() && showNewtonDetails) *out
00247       << "sqrt(phi) = sqrt(<f,f>) = ||f|| = " << sqrt_phi << ( isConverged ? " <= " : " > " ) << "tol = " << tol << endl;
00248     if(out.get() && showNewtonIters) *out
00249       << "newton_iter="<<newtonIter<<": Check convergence: ||f|| = "
00250       << sqrt_phi << ( isConverged ? " <= " : " > " ) << "tol = " << tol << ( isConverged ? ", Converged!!!" : "" ) << endl;
00251     if(isConverged) {
00252       if(x_inout != x.get()) assign( x_inout, *x );  // Assign the solution if we have to
00253       if(out.get() && showNewtonDetails) {
00254         *out << "\nWe have converged :-)\n"
00255              << "\n||x||inf = " << norm_inf(*x) << endl;
00256         if(dumpAll) *out << "\nx =\n" << *x;
00257         *out << "\nExiting SimpleNewtonSolver::solve(...)\n";
00258       }
00259       std::ostringstream oss;
00260       oss << "Converged! ||f|| = " << sqrt_phi << ", num_newton_iters="<<newtonIter<<", num_residual_evals="<<num_residual_evals<<".";
00261       solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
00262       solveStatus.message = oss.str();
00263       break;
00264     }
00265     if(out.get() && showNewtonDetails) *out << "\nWe have to keep going :-(\n";
00266     // Compute the Jacobian if we have not already
00267     if(newtonIter > 1) {
00268       if(out.get() && showNewtonDetails) *out << "\nComputing the Jacobian J_ at current point ...\n";
00269       eval_f_W<Scalar>( *model_, *x, NULL, &*J_ );
00270       if(out.get() && dumpAll) *out << "\nJ =\n" << *J_;
00271     }
00272     // Compute the newton step: dx = -inv(J)*f
00273     if(out.get() && showNewtonDetails) *out << "\nComputing the Newton step: dx = - inv(J)*f ...\n";
00274     if(out.get() && showNewtonIters) *out << "newton_iter="<<newtonIter<<": Computing Newton step ...\n";
00275     assign( &*dx, ST::zero() );       // Initial guess for the linear solve
00276     Thyra::solve(*J_,NOTRANS,*f,&*dx); // Solve: J*dx = f
00277     Vt_S( &*dx, Scalar(-ST::one()) ); // dx *= -1.0
00278     Vp_V( &*ee, *dx);                 // ee += dx
00279     if(out.get() && showNewtonDetails) *out << "\n||dx||inf = " << norm_inf(*dx) << endl;
00280     if(out.get() && dumpAll) *out << "\ndy =\n" << *dx;
00281     // Perform backtracking armijo line search
00282     if(out.get() && showNewtonDetails) *out << "\nStarting backtracking line search iterations ...\n";
00283     if(out.get() && showNewtonIters) *out << "newton_iter="<<newtonIter<<": Starting backtracking line search ...\n";
00284     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
00285     Scalar alpha = 1.0; // Try a full step initially since it will eventually be accepted near solution
00286     int lineSearchIter;
00287     ++num_residual_evals;
00288     for( lineSearchIter = 1; lineSearchIter <= maxLineSearchIterations(); ++lineSearchIter, ++num_residual_evals ) {
00289       TEUCHOS_OSTAB;
00290       if(out.get() && showNewtonDetails) *out << "\n*** lineSearchIter = " << lineSearchIter << endl;
00291       // x_new = x + alpha*dx
00292       assign( &*x_new, *x ); Vp_StV( &*x_new, alpha, *dx );
00293       if(out.get() && showNewtonDetails) *out << "\n||x_new||inf = " << norm_inf(*x_new) << endl;
00294       if(out.get() && dumpAll) *out << "\nx_new =\n" << *x_new;
00295       // Compute the residual at the updated point
00296       eval_f(*model_,*x_new,&*f);
00297       if(out.get() && dumpAll) *out << "\nf_new =\n" << *f;
00298       const Scalar phi_new = scalarProd(*f,*f), phi_frac = phi + alpha * armijoConstant() * Dphi;
00299       if(out.get() && showNewtonDetails) *out << "\nphi_new = <f_new,f_new> = " << phi_new << endl;
00300       if( Teuchos::ScalarTraits<Scalar>::isnaninf(phi_new) ) {
00301         if(out.get() && showNewtonDetails) *out << "\nphi_new is not a valid number, backtracking (alpha = 0.1*alpha) ...\n";
00302         alpha *= 0.1;
00303         continue;
00304       }
00305       const bool acceptPoint = (phi_new <= phi_frac);
00306       if(out.get() && showNewtonDetails) *out
00307         << "\nphi_new = " << phi_new << ( acceptPoint ? " <= " : " > " )
00308         << "phi + alpha * eta * Dphi = " << phi << " + " << alpha << " * " << armijoConstant() << " * " << Dphi
00309         << " = " << phi_frac << endl;
00310       if(out.get() && (showLineSearchIters || (showNewtonIters && acceptPoint))) *out
00311         << "newton_iter="<<newtonIter<<", ls_iter="<<lineSearchIter<<" : "
00312         << "phi(alpha="<<alpha<<") = "<<phi_new<<(acceptPoint ? " <=" : " >")<<" armijo_cord = " << phi_frac << endl;
00313       if( acceptPoint ) {
00314         if(out.get() && showNewtonDetails) *out << "\nAccepting the current step with step length alpha = " << alpha << "!\n";
00315         break;
00316       }
00317       if(out.get() && showNewtonDetails) *out << "\nBacktracking (alpha = 0.5*alpha) ...\n";
00318       alpha *= 0.5;
00319     }
00320     // Check for line search failure
00321     if( lineSearchIter > maxLineSearchIterations() ) {
00322       std::ostringstream oss;
00323       oss
00324         << "lineSearchIter = " << lineSearchIter << " > maxLineSearchIterations = " << maxLineSearchIterations()
00325         << ": Linear search failure!  Algorithm terminated!";
00326       solveStatus.message = oss.str();
00327       if(out.get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
00328       goto exit;
00329     }
00330     // Take the Newton step
00331     std::swap<Teuchos::RefCountPtr<VectorBase<Scalar> > >( x_new, x ); // Now x is current point!
00332   }
00333 exit:
00334   if(out.get() && showNewtonIters) *out
00335     << "\n[Final] newton_iters="<<newtonIter<<", num_residual_evals="<<num_residual_evals<<"\n";
00336   if(newtonIter > maxIters) {
00337     std::ostringstream oss;
00338     oss
00339       << "newton_iter = " << newtonIter << " > maxIters = " << maxIters
00340       << ": Newton algorithm terminated!";
00341     solveStatus.message = oss.str();
00342     if( out.get() && (showNewtonIters || showNewtonDetails)) *out << endl << oss.str() << endl;
00343   }
00344   if(x_inout != x.get()) assign( x_inout, *x ); // Assign the final point
00345   if(delta != NULL) assign( delta, *ee );
00346   current_x_ = x_inout->clone_v(); // Remember the final point
00347   J_is_current_ = newtonIter==1; // J is only current with x if initial point was converged!
00348   if(out.get() && showNewtonDetails) *out
00349     << "\n*** Ending dampended Newton solve." << endl; 
00350   return solveStatus;
00351 }
00352 
00353 template <class Scalar>
00354 Teuchos::RefCountPtr<const VectorBase<Scalar> >
00355 DampenedNewtonNonlinearSolver<Scalar>::get_current_x() const
00356 {
00357   return current_x_;
00358 }
00359 
00360 template <class Scalar>
00361 bool DampenedNewtonNonlinearSolver<Scalar>::is_W_current() const
00362 {
00363   return J_is_current_;
00364 }
00365 
00366 template <class Scalar>
00367 Teuchos::RefCountPtr<LinearOpWithSolveBase<Scalar> >
00368 DampenedNewtonNonlinearSolver<Scalar>::get_nonconst_W()
00369 {
00370   return J_;
00371 }
00372 
00373 template <class Scalar>
00374 Teuchos::RefCountPtr<const LinearOpWithSolveBase<Scalar> >
00375 DampenedNewtonNonlinearSolver<Scalar>::get_W() const
00376 {
00377   return J_;
00378 }
00379 
00380 template <class Scalar>
00381 void DampenedNewtonNonlinearSolver<Scalar>::set_W_is_current(bool W_is_current)
00382 {
00383   J_is_current_ = W_is_current;
00384 }
00385 
00386 } // namespace Thyra
00387 
00388 #endif // THYRA_DAMPENED_NEWTON_NONLINEAR_SOLVER_HPP

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1