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