OptiPack_NonlinearCG_def.hpp

00001 /*
00002 // @HEADER
00003 // ***********************************************************************
00004 // 
00005 //    OptiPack: Collection of simple Thyra-based Optimization ANAs
00006 //                 Copyright (2009) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 // 
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 // 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00026 // 
00027 // ***********************************************************************
00028 // @HEADER
00029 */
00030 
00031 #ifndef OPTIPACK_NONLINEAR_CG_DEF_HPP
00032 #define OPTIPACK_NONLINEAR_CG_DEF_HPP
00033 
00034 
00035 #include "OptiPack_NonlinearCG_decl.hpp"
00036 #include "OptiPack_DefaultPolyLineSearchPointEvaluator.hpp"
00037 #include "OptiPack_UnconstrainedOptMeritFunc1D.hpp"
00038 #include "Thyra_ModelEvaluatorHelpers.hpp"
00039 #include "Thyra_VectorStdOps.hpp"
00040 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00041 #include "Teuchos_StandardParameterEntryValidators.hpp"
00042 #include "Teuchos_Tuple.hpp"
00043 
00044 
00045 namespace OptiPack {
00046 
00047 
00048 template<typename Scalar>
00049 RCP<Teuchos::ParameterEntryValidator>
00050 NonlinearCG<Scalar>::solverType_validator_ = Teuchos::null;
00051 
00052 
00053 // Constructor/Initializers/Accessors
00054 
00055 
00056 template<typename Scalar>
00057 NonlinearCG<Scalar>::NonlinearCG()
00058   : paramIndex_(-1),
00059     responseIndex_(-1),
00060     solverType_(NonlinearCGUtils::solverType_default_integral_val),
00061     alpha_init_(NonlinearCGUtils::alpha_init_default),
00062     alpha_reinit_(NonlinearCGUtils::alpha_reinit_default),
00063     and_conv_tests_(NonlinearCGUtils::and_conv_tests_default),
00064     minIters_(NonlinearCGUtils::minIters_default),
00065     maxIters_(NonlinearCGUtils::maxIters_default),
00066     g_reduct_tol_(NonlinearCGUtils::g_reduct_tol_default),
00067     g_grad_tol_(NonlinearCGUtils::g_grad_tol_default),
00068     g_mag_(NonlinearCGUtils::g_mag_default),
00069     numIters_(0)
00070 {}
00071 
00072 
00073 template<typename Scalar>
00074 void NonlinearCG<Scalar>::initialize(
00075   const RCP<const Thyra::ModelEvaluator<Scalar> > &model,
00076   const int paramIndex,
00077   const int responseIndex,
00078   const RCP<GlobiPack::LineSearchBase<Scalar> > &linesearch
00079   )
00080 {
00081   // ToDo: Validate input objects!
00082   model_ = model.assert_not_null();
00083   paramIndex_ = paramIndex;
00084   responseIndex_ = responseIndex;
00085   linesearch_ = linesearch.assert_not_null();
00086 }
00087 
00088 
00089 template<typename Scalar>
00090 NonlinearCGUtils::ESolverTypes NonlinearCG<Scalar>::get_solverType() const
00091 {
00092   return solverType_;
00093 }
00094 
00095 
00096 template<typename Scalar>
00097 typename NonlinearCG<Scalar>::ScalarMag
00098 NonlinearCG<Scalar>::get_alpha_init() const
00099 {
00100   return alpha_init_;
00101 }
00102 
00103 
00104 template<typename Scalar>
00105 bool NonlinearCG<Scalar>::get_alpha_reinit() const
00106 {
00107   return alpha_reinit_;
00108 }
00109 
00110 
00111 template<typename Scalar>
00112 bool NonlinearCG<Scalar>::get_and_conv_tests() const
00113 {
00114   return and_conv_tests_;
00115 }
00116 
00117 
00118 template<typename Scalar>
00119 int NonlinearCG<Scalar>::get_minIters() const
00120 {
00121   return minIters_;
00122 }
00123 
00124 
00125 template<typename Scalar>
00126 int NonlinearCG<Scalar>::get_maxIters() const
00127 {
00128   return maxIters_;
00129 }
00130 
00131 
00132 template<typename Scalar>
00133 typename NonlinearCG<Scalar>::ScalarMag
00134 NonlinearCG<Scalar>::get_g_reduct_tol() const
00135 {
00136   return g_reduct_tol_;
00137 }
00138 
00139 
00140 template<typename Scalar>
00141 typename NonlinearCG<Scalar>::ScalarMag
00142 NonlinearCG<Scalar>::get_g_grad_tol() const
00143 {
00144   return g_grad_tol_;
00145 }
00146 
00147 
00148 template<typename Scalar>
00149 typename NonlinearCG<Scalar>::ScalarMag
00150 NonlinearCG<Scalar>::get_g_mag() const
00151 {
00152   return g_mag_;
00153 }
00154 
00155 
00156 // Overridden from ParameterListAcceptor (simple forwarding functions)
00157 
00158 
00159 template<typename Scalar>
00160 void NonlinearCG<Scalar>::setParameterList(RCP<ParameterList> const& paramList)
00161 {
00162   typedef ScalarTraits<Scalar> ST;
00163   typedef ScalarTraits<ScalarMag> SMT;
00164   namespace NCGU = NonlinearCGUtils;
00165   using Teuchos::getParameter;
00166   using Teuchos::getIntegralValue;
00167   paramList->validateParametersAndSetDefaults(*this->getValidParameters());
00168   solverType_ = getIntegralValue<NCGU::ESolverTypes>(*paramList, NCGU::solverType_name);
00169   alpha_init_ = getParameter<double>(*paramList, NCGU::alpha_init_name);
00170   alpha_reinit_ = getParameter<bool>(*paramList, NCGU::alpha_reinit_name);
00171   and_conv_tests_ = getParameter<bool>(*paramList, NCGU::and_conv_tests_name);
00172   minIters_ = getParameter<int>(*paramList, NCGU::minIters_name);
00173   maxIters_ = getParameter<int>(*paramList, NCGU::maxIters_name);
00174   g_reduct_tol_ = getParameter<double>(*paramList, NCGU::g_reduct_tol_name);
00175   g_grad_tol_ = getParameter<double>(*paramList, NCGU::g_grad_tol_name);
00176   g_mag_ = getParameter<double>(*paramList, NCGU::g_mag_name);
00177   TEUCHOS_ASSERT_INEQUALITY( alpha_init_, >, SMT::zero() );
00178   TEUCHOS_ASSERT_INEQUALITY( minIters_, >=, 0 );
00179   TEUCHOS_ASSERT_INEQUALITY( minIters_, <, maxIters_ );
00180   TEUCHOS_ASSERT_INEQUALITY( g_reduct_tol_, >=, SMT::zero() );
00181   TEUCHOS_ASSERT_INEQUALITY( g_grad_tol_, >=, SMT::zero() );
00182   TEUCHOS_ASSERT_INEQUALITY( g_mag_, >, SMT::zero() );
00183   Teuchos::readVerboseObjectSublist(&*paramList, this);
00184   setMyParamList(paramList);
00185 }
00186 
00187 
00188 template<typename Scalar>
00189 RCP<const ParameterList>
00190 NonlinearCG<Scalar>::getValidParameters() const
00191 {
00192   using Teuchos::tuple;
00193   namespace NCGU = NonlinearCGUtils;
00194   static RCP<const ParameterList> validPL;
00195   if (is_null(validPL)) {
00196     RCP<Teuchos::ParameterList>
00197       pl = Teuchos::rcp(new Teuchos::ParameterList());
00198     solverType_validator_ = 
00199       Teuchos::stringToIntegralParameterEntryValidator<NCGU::ESolverTypes>(
00200         tuple<std::string>(
00201           "FR",
00202           "PR+",
00203           "FR-PR",
00204           "HS"
00205           ),
00206         tuple<std::string>(
00207           "Fletcher-Reeves Method",
00208           "Polak-Ribiere Method",
00209           "Fletcher-Reeves Polak-Ribiere Hybrid Method",
00210           "Hestenes-Stiefel Method"
00211           ),
00212         tuple<NCGU::ESolverTypes>(
00213           NCGU::NONLINEAR_CG_FR,
00214           NCGU::NONLINEAR_CG_PR_PLUS,
00215           NCGU::NONLINEAR_CG_FR_PR,
00216           NCGU::NONLINEAR_CG_HS
00217           ),
00218         NCGU::solverType_name
00219         );
00220     pl->set( NCGU::solverType_name, NCGU::solverType_default,
00221       "Set the type of nonlinear CG solver algorithm to use.",
00222       solverType_validator_ );
00223     pl->set( NCGU::alpha_init_name, NCGU::alpha_init_default );
00224     pl->set( NCGU::alpha_reinit_name, NCGU::alpha_reinit_default );
00225     pl->set( NCGU::and_conv_tests_name, NCGU::and_conv_tests_default );
00226     pl->set( NCGU::minIters_name, NCGU::minIters_default );
00227     pl->set( NCGU::maxIters_name, NCGU::maxIters_default );
00228     pl->set( NCGU::g_reduct_tol_name, NCGU::g_reduct_tol_default );
00229     pl->set( NCGU::g_grad_tol_name, NCGU::g_grad_tol_default );
00230     pl->set( NCGU::g_mag_name, NCGU::g_mag_default );
00231     Teuchos::setupVerboseObjectSublist(&*pl);
00232     validPL = pl;
00233     // ToDo: Add documentation for these parameters
00234   }
00235   return validPL;
00236 }
00237 
00238 
00239 // Solve
00240 
00241 
00242 template<typename Scalar>
00243 NonlinearCGUtils::ESolveReturn
00244 NonlinearCG<Scalar>::doSolve(
00245   const Ptr<Thyra::VectorBase<Scalar> > &p_inout,
00246   const Ptr<ScalarMag> &g_opt_out,
00247   const Ptr<const ScalarMag> &g_reduct_tol_in,
00248   const Ptr<const ScalarMag> &g_grad_tol_in,
00249   const Ptr<const ScalarMag> &alpha_init_in,
00250   const Ptr<int> &numIters_out
00251   )
00252 {
00253 
00254   typedef ScalarTraits<Scalar> ST;
00255   typedef ScalarTraits<ScalarMag> SMT;
00256   
00257   using Teuchos::null;
00258   using Teuchos::as;
00259   using Teuchos::tuple;
00260   using Teuchos::rcpFromPtr;
00261   using Teuchos::optInArg;
00262   using Teuchos::inOutArg;
00263   using GlobiPack::computeValue;
00264   using GlobiPack::PointEval1D;
00265   using Thyra::VectorSpaceBase;
00266   using Thyra::VectorBase;
00267   using Thyra::MultiVectorBase;
00268   using Thyra::scalarProd;
00269   using Thyra::createMember;
00270   using Thyra::createMembers;
00271   using Thyra::get_ele;
00272   using Thyra::norm;
00273   using Thyra::V_StV;
00274   using Thyra::Vt_S;
00275   using Thyra::eval_g_DgDp;
00276   typedef Thyra::Ordinal Ordinal;
00277   typedef Thyra::ModelEvaluatorBase MEB;
00278   namespace NCGU = NonlinearCGUtils;
00279   using std::max;
00280 
00281   // Validate input
00282 
00283   g_opt_out.assert_not_null();
00284 
00285   // Set streams
00286 
00287   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00288   linesearch_->setOStream(out);
00289 
00290   // Determine what step constants will be computed
00291 
00292   const bool compute_beta_PR =
00293     (
00294       solverType_ == NCGU::NONLINEAR_CG_PR_PLUS
00295       ||
00296       solverType_ == NCGU::NONLINEAR_CG_FR_PR
00297       );
00298 
00299   const bool compute_beta_HS = (solverType_ == NCGU::NONLINEAR_CG_HS);
00300 
00301   //
00302   // A) Set up the storage for the algorithm
00303   //
00304   
00305   const RCP<DefaultPolyLineSearchPointEvaluator<Scalar> >
00306     pointEvaluator = defaultPolyLineSearchPointEvaluator<Scalar>();
00307 
00308   const RCP<UnconstrainedOptMeritFunc1D<Scalar> >
00309     meritFunc = unconstrainedOptMeritFunc1D<Scalar>(
00310       model_, paramIndex_, responseIndex_ );
00311 
00312   const RCP<const VectorSpaceBase<Scalar> >
00313     p_space = model_->get_p_space(paramIndex_),
00314     g_space = model_->get_g_space(responseIndex_);
00315 
00316   // Stoarge for current iteration
00317   RCP<VectorBase<Scalar> >
00318     p_k = rcpFromPtr(p_inout),        // Current solution for p
00319     p_kp1 = createMember(p_space),    // Trial point for p (in line search)
00320     g_vec = createMember(g_space),    // Vector (size 1) form of objective g(p) 
00321     g_grad_k = createMember(p_space), // Gradient of g DgDp^T
00322     d_k = createMember(p_space),      // Search direction
00323     g_grad_k_diff_km1 = null;         // g_grad_k - g_grad_km1 (if needed)
00324 
00325   // Storage for previous iteration
00326   RCP<VectorBase<Scalar> >
00327     g_grad_km1 = null, // Will allocate if we need it!
00328     d_km1 = null; // Will allocate if we need it!
00329   ScalarMag
00330     alpha_km1 = SMT::zero(),
00331     g_km1 = SMT::zero(),
00332     g_grad_km1_inner_g_grad_km1 = SMT::zero(),
00333     g_grad_km1_inner_d_km1 = SMT::zero();
00334   
00335   if (compute_beta_PR || compute_beta_HS) {
00336     g_grad_km1 = createMember(p_space);
00337     g_grad_k_diff_km1 = createMember(p_space);
00338   }
00339   
00340   if (compute_beta_HS) {
00341     d_km1 = createMember(p_space);
00342   }
00343 
00344   //
00345   // B) Do the nonlinear CG iterations
00346   //
00347 
00348   *out << "\nStarting nonlinear CG iterations ...\n";
00349 
00350   if (and_conv_tests_) {
00351     *out << "\nNOTE: Using AND of convergence tests!\n";
00352   }
00353   else {
00354     *out << "\nNOTE: Using OR of convergence tests!\n";
00355   }
00356 
00357   const Scalar alpha_init =
00358     ( !is_null(alpha_init_in) ? *alpha_init_in : alpha_init_ );
00359   const Scalar g_reduct_tol =
00360     ( !is_null(g_reduct_tol_in) ? *g_reduct_tol_in : g_reduct_tol_ );
00361   const Scalar g_grad_tol =
00362     ( !is_null(g_grad_tol_in) ? *g_grad_tol_in : g_grad_tol_ );
00363 
00364   const Ordinal globalDim = p_space->dim();
00365 
00366   bool foundSolution = false;
00367   bool fatalLinesearchFailure = false;
00368   bool restart = true;
00369   int numConsecutiveLineSearchFailures = 0;
00370 
00371   int numConsecutiveIters = 0;
00372 
00373   for (numIters_ = 0; numIters_ < maxIters_; ++numIters_, ++numConsecutiveIters) {
00374 
00375     Teuchos::OSTab tab(out);
00376 
00377     *out << "\nNonlinear CG Iteration k = " << numIters_ << "\n";
00378 
00379     Teuchos::OSTab tab2(out);
00380 
00381     //
00382     // B.1) Evaluate the point (on first iteration)
00383     //
00384     
00385     eval_g_DgDp(
00386       *model_, paramIndex_, *p_k, responseIndex_,
00387       numIters_ == 0 ? g_vec.ptr() : null, // Only on first iteration
00388       MEB::Derivative<Scalar>(g_grad_k, MEB::DERIV_MV_GRADIENT_FORM) );
00389 
00390     const ScalarMag g_k = get_ele(*g_vec, 0);
00391     // Above: If numIters_ > 0, then g_vec was updated in meritFunc->eval(...).
00392 
00393     //
00394     // B.2) Check for convergence
00395     //
00396 
00397     // B.2.a) ||g_k - g_km1|| |g_k + g_mag| <= g_reduct_tol
00398 
00399     bool g_reduct_converged = false;
00400 
00401     if (numIters_ > 0) {
00402 
00403       const ScalarMag g_reduct = g_k - g_km1;
00404       
00405       *out << "\ng_k - g_km1 = "<<g_reduct<<"\n";
00406       
00407       const ScalarMag g_reduct_err =
00408         SMT::magnitude(g_reduct / SMT::magnitude(g_k + g_mag_));
00409       
00410       g_reduct_converged = (g_reduct_err <= g_reduct_tol);
00411       
00412       *out << "\nCheck convergence: |g_k - g_km1| / |g_k + g_mag| = "<<g_reduct_err
00413            << (g_reduct_converged ? " <= " : " > ")
00414            << "g_reduct_tol = "<<g_reduct_tol<<"\n";
00415       
00416     }
00417 
00418     // B.2.b) ||g_grad_k|| g_mag <= g_grad_tol
00419 
00420     const Scalar g_grad_k_inner_g_grad_k = scalarProd<Scalar>(*g_grad_k, *g_grad_k);
00421     const ScalarMag norm_g_grad_k = ST::magnitude(ST::squareroot(g_grad_k_inner_g_grad_k));
00422 
00423     *out << "\n||g_grad_k|| = "<<norm_g_grad_k << "\n";
00424 
00425     const ScalarMag g_grad_err = norm_g_grad_k / g_mag_;
00426 
00427     const bool g_grad_converged = (g_grad_err <= g_grad_tol);
00428 
00429     *out << "\nCheck convergence: ||g_grad_k|| / g_mag = "<<g_grad_err
00430          << (g_grad_converged ? " <= " : " > ")
00431          << "g_grad_tol = "<<g_grad_tol<<"\n";
00432 
00433     // B.2.c) Convergence status
00434     
00435     bool isConverged = false;
00436     if (and_conv_tests_) {
00437       isConverged = g_reduct_converged && g_grad_converged;
00438     }
00439     else {
00440       isConverged = g_reduct_converged || g_grad_converged;
00441     }
00442 
00443     if (isConverged) {
00444       if (numIters_ < minIters_) {
00445         *out << "\nnumIters="<<numIters_<<" < minIters="<<minIters_
00446              << ", continuing on!\n";
00447       }
00448       else {
00449         *out << "\nFound solution, existing algorithm!\n";
00450         foundSolution = true;
00451       }
00452     }
00453     else {
00454       *out << "\nNot converged!\n";
00455     }
00456     
00457     if (foundSolution) {
00458       break;
00459     }
00460 
00461     //
00462     // B.3) Compute the search direction d_k
00463     //
00464 
00465     if (numConsecutiveIters == globalDim) {
00466 
00467       *out << "\nThe number of consecutive iterations exceeds the"
00468            << " global dimension so restarting!\n";
00469 
00470       restart = true;
00471 
00472     }
00473 
00474     if (restart) {
00475 
00476       *out << "\nResetting search direction back to steppest descent!\n";
00477 
00478       // d_k = -g_grad_k
00479       V_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
00480 
00481       restart = false;
00482 
00483     }
00484     else {
00485       
00486       // g_grad_k - g_grad_km1
00487       if (!is_null(g_grad_k_diff_km1)) {
00488         V_VmV( g_grad_k_diff_km1.ptr(), *g_grad_k, *g_grad_km1 );
00489       }
00490 
00491       // beta_FR = inner(g_grad_k, g_grad_k) / inner(g_grad_km1, g_grad_km1)
00492       const Scalar beta_FR =
00493         g_grad_k_inner_g_grad_k / g_grad_km1_inner_g_grad_km1;
00494       *out << "\nbeta_FR = " << beta_FR << "\n";
00495       // NOTE: Computing beta_FR is free so we might as well just do it!
00496 
00497       // beta_PR = inner(g_grad_k, g_grad_k - g_grad_km1) /
00498       //    inner(g_grad_km1, g_grad_km1)
00499       Scalar beta_PR = ST::zero();
00500       if (compute_beta_PR) {
00501         beta_PR =
00502           inner(*g_grad_k, *g_grad_k_diff_km1) / g_grad_km1_inner_g_grad_km1;
00503         *out << "\nbeta_PR = " << beta_PR << "\n";
00504       }
00505 
00506       // beta_HS = inner(g_grad_k, g_grad_k - g_grad_km1) /
00507       //    inner(g_grad_k - g_grad_km1, d_km1)
00508       Scalar beta_HS = ST::zero();
00509       if (compute_beta_HS) {
00510         beta_HS =
00511           inner(*g_grad_k, *g_grad_k_diff_km1) / inner(*g_grad_k_diff_km1, *d_km1);
00512         *out << "\nbeta_HS = " << beta_HS << "\n";
00513       }
00514       
00515       Scalar beta_k = ST::zero();
00516       switch(solverType_) {
00517         case NCGU::NONLINEAR_CG_FR: {
00518           beta_k = beta_FR;
00519           break;
00520         }
00521         case NCGU::NONLINEAR_CG_PR_PLUS: {
00522           beta_k = max(beta_PR, ST::zero());
00523           break;
00524         }
00525         case NCGU::NONLINEAR_CG_FR_PR: {
00526           // NOTE: This does not seem to be working :-(
00527           if (numConsecutiveIters < 2) {
00528             beta_k = beta_PR;
00529           }
00530           else if (beta_PR < -beta_FR)
00531             beta_k = -beta_FR;
00532           else if (ST::magnitude(beta_PR) <= beta_FR)
00533             beta_k = beta_PR;
00534           else // beta_PR > beta_FR
00535             beta_k = beta_FR;
00536         }
00537         case NCGU::NONLINEAR_CG_HS: {
00538           beta_k = beta_HS;
00539           break;
00540         }
00541         default:
00542           TEST_FOR_EXCEPT(true);
00543       }
00544       *out << "\nbeta_k = " << beta_k << "\n";
00545 
00546       // d_k = beta_k * d_last + -g_grad_k
00547       if (!is_null(d_km1))
00548         V_StV( d_k.ptr(), beta_k, *d_km1 );
00549       else
00550         Vt_S( d_k.ptr(), beta_k );
00551       Vp_StV( d_k.ptr(), as<Scalar>(-1.0), *g_grad_k );
00552 
00553     }
00554     
00555     //
00556     // B.4) Perform the line search
00557     //
00558 
00559     // B.4.a) Compute the initial step length
00560 
00561     Scalar alpha_k = as<Scalar>(-1.0);
00562 
00563     if (numIters_ == 0) {
00564       alpha_k = alpha_init;
00565     }
00566     else {
00567       if (alpha_reinit_) {
00568         alpha_k = alpha_init;
00569       }
00570       else {
00571         alpha_k = alpha_km1;
00572         // ToDo: Implement better logic from Nocedal and Wright for selecting
00573         // this step length after first iteration!
00574       }
00575     }
00576 
00577     // B.4.b) Perform the linesearch (computing updated quantities in process)
00578 
00579     pointEvaluator->initialize(tuple<RCP<const VectorBase<Scalar> > >(p_k, d_k)());
00580 
00581     ScalarMag g_grad_k_inner_d_k = ST::zero();
00582 
00583     // Set up the merit function to only compute the value
00584     meritFunc->setEvaluationQuantities(pointEvaluator, p_kp1, g_vec, null);
00585 
00586     PointEval1D<ScalarMag> point_k(ST::zero(), g_k);
00587     if (linesearch_->requiresBaseDeriv()) {
00588       g_grad_k_inner_d_k = scalarProd(*g_grad_k, *d_k);
00589       point_k.Dphi = g_grad_k_inner_d_k;
00590     }
00591 
00592     ScalarMag g_kp1 = computeValue(*meritFunc, alpha_k);
00593     // NOTE: The above call updates p_kp1 and g_vec as well!
00594 
00595     PointEval1D<ScalarMag> point_kp1(alpha_k, g_kp1);
00596 
00597     const bool linesearchResult = linesearch_->doLineSearch(
00598       *meritFunc, point_k, inOutArg(point_kp1), null );
00599 
00600     alpha_k = point_kp1.alpha;
00601     g_kp1 = point_kp1.phi;
00602 
00603     if (linesearchResult) {
00604       numConsecutiveLineSearchFailures = 0;
00605     }
00606     else {
00607       if (numConsecutiveLineSearchFailures==0) {
00608         *out << "\nLine search failure, resetting the search direction!\n";
00609         restart = true;
00610       }
00611       if (numConsecutiveLineSearchFailures==1) {
00612         *out << "\nLine search failure on last iteration also, terminating algorithm!\n";
00613         fatalLinesearchFailure = true;
00614       }
00615       ++numConsecutiveLineSearchFailures;
00616     }
00617 
00618     if (fatalLinesearchFailure) {
00619       break;
00620     }
00621 
00622     //
00623     // B.5) Transition to the next iteration
00624     //
00625     
00626     alpha_km1 = alpha_k;
00627     g_km1 = g_k;
00628     g_grad_km1_inner_g_grad_km1 = g_grad_k_inner_g_grad_k;
00629     g_grad_km1_inner_d_km1 = g_grad_k_inner_d_k;
00630     std::swap(p_k, p_kp1);
00631     if (!is_null(g_grad_km1))
00632       std::swap(g_grad_km1, g_grad_k);
00633     if (!is_null(d_km1))
00634       std::swap(d_k, d_km1);
00635     
00636 #ifdef TEUCHOS_DEBUG
00637     // Make sure we compute these correctly before they are used!
00638     V_S(g_grad_k.ptr(), ST::nan());
00639     V_S(p_kp1.ptr(), ST::nan());
00640 #endif
00641 
00642   }
00643 
00644   //
00645   // C) Final clean up
00646   //
00647   
00648   // Get the most current value of g(p)
00649   *g_opt_out = get_ele(*g_vec, 0);
00650 
00651   // Make sure that the final value for p has been copied in!
00652   V_V( p_inout, *p_k );
00653 
00654   if (!is_null(numIters_out)) {
00655     *numIters_out = numIters_;
00656   }
00657 
00658   if (numIters_ == maxIters_) {
00659     *out << "\nMax nonlinear CG iterations exceeded!\n";
00660   }
00661   
00662   if (foundSolution) {
00663     return NonlinearCGUtils::SOLVE_SOLUTION_FOUND;
00664   }
00665   else if(fatalLinesearchFailure) {
00666     return NonlinearCGUtils::SOLVE_LINSEARCH_FAILURE;
00667   }
00668 
00669   // Else, the max number of iterations was exceeded
00670   return NonlinearCGUtils::SOLVE_MAX_ITERS_EXCEEDED;
00671 
00672 }
00673 
00674 
00675 } // namespace OptiPack
00676 
00677 
00678 #endif // OPTIPACK_NONLINEAR_CG_DEF_HPP

Generated on Tue Jul 13 09:37:12 2010 for OptiPack by  doxygen 1.4.7