Changeset 551


Ignore:
Timestamp:
Oct 26, 2005 8:31:28 PM (14 years ago)
Author:
andreasw
Message:
  • improvement to the CG-penatly implememtation
Location:
branches/dev-penalty/Algorithm
Files:
10 edited

Legend:

Unmodified
Added
Removed
  • branches/dev-penalty/Algorithm/IpBacktrackingLineSearch.cpp

    r549 r551  
    530530  }
    531531
    532   bool BacktrackingLineSearch::DoBacktrackingLineSearch(bool skip_first_trial_point,
    533       Number& alpha_primal,
    534       bool& corr_taken,
    535       bool& soc_taken,
    536       Index& n_steps,
    537       bool& evaluation_error,
    538       SmartPtr<IteratesVector>& actual_delta)
    539   {
     532  bool BacktrackingLineSearch::DoBacktrackingLineSearch(
     533    bool skip_first_trial_point,
     534    Number& alpha_primal,
     535    bool& corr_taken,
     536    bool& soc_taken,
     537    Index& n_steps,
     538    bool& evaluation_error,
     539    SmartPtr<IteratesVector>& actual_delta)
     540  {
     541    DBG_START_METH("BacktrackingLineSearch::DoBacktrackingLineSearch",
     542                   dbg_verbosity);
     543
    540544    evaluation_error = false;
    541545    bool accept = false;
    542 
    543     DBG_START_METH("BacktrackingLineSearch::DoBacktrackingLineSearch",
    544                    dbg_verbosity);
    545546
    546547    // Compute primal fraction-to-the-boundary value
     
    561562    alpha_primal = alpha_primal_max;
    562563
    563     // Step size used in ftype and armijo tests
     564    // Step size used in Armijo tests
    564565    Number alpha_primal_test = alpha_primal;
    565566    if (in_watchdog_) {
     
    746747  }
    747748
    748   void BacktrackingLineSearch::PerformDualStep(Number alpha_primal,
    749       Number alpha_dual,
    750       SmartPtr<IteratesVector>& delta)
     749  void BacktrackingLineSearch::PerformDualStep(
     750    Number alpha_primal,
     751    Number alpha_dual,
     752    SmartPtr<IteratesVector>& delta)
    751753  {
    752754    DBG_START_FUN("BacktrackingLineSearch::PerformDualStep", dbg_verbosity);
     
    770772      break;
    771773      case FULL_STEP_FOR_Y:
    772       alpha_y = 1;
     774      alpha_y = 1.;
    773775      break;
    774776      case MIN_DUAL_INFEAS_ALPHA_FOR_Y:
  • branches/dev-penalty/Algorithm/IpCGPenaltyLSAcceptor.cpp

    r550 r551  
    5151      "Relaxation factor in the Armijo condition for the penalty function.",
    5252      0.0, true, 0.5, true, 1e-8);
     53    roptions->AddLowerBoundedNumberOption(
     54      "penalty_update_infeasibility_tol",
     55      "Threshold for infeasibility in penalty parameter update test.",
     56      0.0, true, 1e-9,
     57      "If the new constraint violation is smaller than this tolerance, the "
     58      "penalty parameter is not increased.");
     59    roptions->AddLowerBoundedNumberOption(
     60      "eta_min",
     61      "LIFENG WRITES THIS.",
     62      0.0, true, 1e-2,
     63      "");
     64    roptions->AddLowerBoundedNumberOption(
     65      "penalty_update_compl_tol",
     66      "LIFENG WRITES THIS.",
     67      0.0, true, 0.2,
     68      "");
     69    roptions->AddLowerBoundedNumberOption(
     70      "chi_hat",
     71      "LIFENG WRITES THIS.",
     72      0.0, true, 2.,
     73      "");
     74    roptions->AddLowerBoundedNumberOption(
     75      "chi_tilde",
     76      "LIFENG WRITES THIS.",
     77      0.0, true, 5.,
     78      "");
     79    roptions->AddLowerBoundedNumberOption(
     80      "chi_cup",
     81      "LIFENG WRITES THIS.",
     82      0.0, true, 1.5,
     83      "");
     84    roptions->AddLowerBoundedNumberOption(
     85      "gamma_hat",
     86      "LIFENG WRITES THIS.",
     87      0.0, true, 0.04,
     88      "");
     89    roptions->AddLowerBoundedNumberOption(
     90      "gamma_tilde",
     91      "LIFENG WRITES THIS.",
     92      0.0, true, 4.,
     93      "");
     94    roptions->AddLowerBoundedNumberOption(
     95      "penalty_max",
     96      "LIFENG WRITES THIS.",
     97      0.0, true, 1e20,
     98      "");
     99    roptions->AddLowerBoundedNumberOption(
     100      "epsilon_c",
     101      "LIFENG WRITES THIS.",
     102      0.0, true, 1e-2,
     103      "");
    53104  }
    54105
     
    57108  {
    58109    options.GetNumericValue("eta_penalty", eta_penalty_, prefix);
    59 
    60     // The following could become regular options
    61     penalty_update_infeasibility_tol_ = 1e-9;
    62     eta_min_ = 1e-2;
    63     penalty_update_compl_tol_ = 0.2;
    64     chi_hat_ = 2.;
    65     chi_tilde_ = 5.;
    66     chi_cup_ = 1.5;
    67     gamma_hat_ = 0.04;
    68     gamma_tilde_ = 4.;
    69     penalty_max_ = 1e20;
    70     epsilon_c_ = 1e-2;
     110    options.GetNumericValue("penalty_update_infeasibility_tol",
     111                            penalty_update_infeasibility_tol_, prefix);
     112    options.GetNumericValue("eta_min", eta_min_, prefix);
     113    options.GetNumericValue("penalty_update_compl_tol",
     114                            penalty_update_compl_tol_, prefix);
     115    options.GetNumericValue("chi_hat", chi_hat_, prefix);
     116    options.GetNumericValue("chi_tilde", chi_tilde_, prefix);
     117    options.GetNumericValue("chi_cup", chi_cup_, prefix);
     118    options.GetNumericValue("gamma_hat", gamma_hat_, prefix);
     119    options.GetNumericValue("gamma_tilde", gamma_tilde_, prefix);
     120    options.GetNumericValue("penalty_max", penalty_max_, prefix);
     121    options.GetNumericValue("epsilon_c", epsilon_c_, prefix);
     122    // The following two option is registered by FilterLSAcceptor
     123    options.GetIntegerValue("max_soc", max_soc_, prefix);
     124    if (max_soc_>0) {
     125      ASSERT_EXCEPTION(IsValid(pd_solver_), OPTION_INVALID,
     126                       "Option \"max_soc\": This option is non-negative, but no linear solver for computing the SOC given to FilterLSAcceptor object.");
     127    }
     128    options.GetNumericValue("kappa_soc", kappa_soc_, prefix);
    71129
    72130    Reset();
     
    87145    if (!in_watchdog) {
    88146      reference_penalty_function_ = IpCq().curr_penalty_function();
    89       // ToDo we need to make this cleaner!
    90       if (IpData().HaveCgPenDeltas()) {
     147      if (IpData().HaveCgFastDeltas()) {
    91148        // use the fast step
    92149        reference_direct_deriv_penalty_function_ =
    93150          IpCq().curr_fast_direct_deriv_penalty_function();
    94         // Overwrite the original step (hmm.....)
    95         SmartPtr<const IteratesVector> delta_cgpen = IpData().delta_cgpen();
    96         IpData().set_delta(delta_cgpen);
    97151      }
    98152      else {
     
    165219    watchdog_direct_deriv_penalty_function_ =
    166220      IpCq().curr_direct_deriv_penalty_function();
     221    watchdog_delta_cgpen_ = IpData().delta_cgpen();
    167222  }
    168223
     
    174229    reference_direct_deriv_penalty_function_ =
    175230      watchdog_direct_deriv_penalty_function_;
     231    IpData().set_delta_cgpen(watchdog_delta_cgpen_);
     232    watchdog_delta_cgpen_ = NULL;
    176233  }
    177234
     
    195252                   dbg_verbosity);
    196253
    197     return false;
     254    if (max_soc_==0) {
     255      return false;
     256    }
     257
     258    bool accept = false;
     259    Index count_soc = 0;
     260
     261    Number theta_soc_old = 0.;
     262    Number theta_trial = IpCq().trial_constraint_violation();
     263    Number alpha_primal_soc = alpha_primal;
     264
     265    // delta_y_c and delta_y_d are the steps used in the right hand
     266    // side for the SOC step
     267    SmartPtr<const Vector> delta_y_c = IpData().delta()->y_c();
     268    SmartPtr<const Vector> delta_y_d = IpData().delta()->y_d();
     269
     270    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNewCopy();
     271    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNewCopy();
     272
     273    while (count_soc<max_soc_ && !accept &&
     274           (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
     275      theta_soc_old = theta_trial;
     276
     277      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
     278                     "Trying second order correction number %d\n",
     279                     count_soc+1);
     280
     281      // Compute SOC constraint violation
     282      Number c_over_r = IpCq().curr_cg_pert_fact();
     283      c_soc->AddTwoVectors(1.0, *IpCq().trial_c(),
     284                           -c_over_r, *delta_y_c,
     285                           alpha_primal_soc);
     286      dms_soc->AddTwoVectors(1.0, *IpCq().trial_d_minus_s(),
     287                             -c_over_r, *delta_y_d,
     288                             alpha_primal_soc);
     289
     290      // Compute the SOC search direction
     291      SmartPtr<IteratesVector> delta_soc =
     292        actual_delta->MakeNewIteratesVector(true);
     293      SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
     294      rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
     295      rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
     296      rhs->Set_y_c(*c_soc);
     297      rhs->Set_y_d(*dms_soc);
     298      rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
     299      rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
     300      rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
     301      rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
     302      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
     303
     304      // Update the delta_y_c and delta_y_d vectors in case we do
     305      // additional SOC steps
     306      delta_y_c = ConstPtr(delta_soc->y_c());
     307      delta_y_d = ConstPtr(delta_soc->y_d());
     308
     309      // Compute step size
     310      alpha_primal_soc =
     311        IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
     312                                        *delta_soc->x(),
     313                                        *delta_soc->s());
     314
     315      // Check if trial point is acceptable
     316      try {
     317        // Compute the primal trial point
     318        IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());
     319
     320        // in acceptance tests, use original step size!
     321        accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
     322      }
     323      catch(IpoptNLP::Eval_Error& e) {
     324        e.ReportException(Jnlst(), J_DETAILED);
     325        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
     326        IpData().Append_info_string("e");
     327        accept = false;
     328        // There is no point in continuing SOC procedure
     329        break;
     330      }
     331
     332      if (accept) {
     333        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
     334        // Accept all SOC quantities
     335        alpha_primal = alpha_primal_soc;
     336        actual_delta = delta_soc;
     337      }
     338      else {
     339        count_soc++;
     340        theta_trial = IpCq().trial_constraint_violation();
     341      }
     342    }
     343    return accept;
    198344  }
    199345
     
    216362    if (curr_eta_<0.) {
    217363      // We need to initialize the eta tolerance
    218       // Lifeng, at the moment I'm using the scaled 1-norm for E_0
    219364      curr_eta_ = Max(eta_min_, Min(gamma_tilde_,
    220                                     gamma_hat_*IpCq().curr_primal_dual_system_error(0.)));
     365                                    gamma_hat_*IpCq().curr_nlp_error()));
    221366    }
    222367
    223368    // Check if the penalty parameter is to be increased
    224 
    225369    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    226370                   "Starting tests for penalty parameter update:\n");
    227     // Lifeng: Should we use the new infeasibility here?
    228     Number curr_inf = IpCq().curr_primal_infeasibility(NORM_2);
     371
     372    // We use the new infeasibility here...
     373    Number trial_inf = IpCq().trial_primal_infeasibility(NORM_2);
    229374    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    230                    "current infeasibility = %8.2\n", curr_inf);
    231     bool increase = (curr_inf >= penalty_update_infeasibility_tol_);
     375                   "trial infeasibility = %8.2\n", trial_inf);
     376    bool increase = (trial_inf >= penalty_update_infeasibility_tol_);
    232377    if (!increase) {
    233378      info_alpha_primal_char='i';
     
    235380
    236381    if (increase) {
    237       Number max_step = Max(IpData().delta()->x()->Amax(),
    238                             IpData().delta()->s()->Amax());
     382      Number max_step = Max(IpData().delta_cgpen()->x()->Amax(),
     383                            IpData().delta_cgpen()->s()->Amax());
    239384      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    240385                     "Max norm of step = %8.2\n", max_step);
     
    286431    if (increase) {
    287432      SmartPtr<Vector> vec = IpData().curr()->y_c()->MakeNewCopy();
    288       vec->AddTwoVectors(1., *IpData().delta()->y_c(),
     433      vec->AddTwoVectors(1., *IpData().delta_cgpen()->y_c(),
    289434                         -1./IpCq().curr_cg_pert_fact(), *IpCq().curr_c(),
    290435                         1.);
     
    320465      SmartPtr<Vector> vec = IpData().curr()->y_c()->MakeNew();
    321466      vec->AddTwoVectors(1., *IpData().curr()->y_c(),
    322                          1., *IpData().delta()->y_c(), 0.);
     467                         1., *IpData().delta_cgpen()->y_c(), 0.);
    323468      y_full_step_max = vec->Amax();
    324469      vec = IpData().curr()->y_d()->MakeNew();
    325470      vec->AddTwoVectors(1., *IpData().curr()->y_d(),
    326                          1., *IpData().delta()->y_d(), 0.);
     471                         1., *IpData().delta_cgpen()->y_d(), 0.);
    327472      y_full_step_max = Max(y_full_step_max, vec->Amax());
    328473      if (IpCq().curr_primal_infeasibility(NORM_2) >= epsilon_c_) {
  • branches/dev-penalty/Algorithm/IpCGPenaltyLSAcceptor.hpp

    r549 r551  
    153153    Number penalty_max_;
    154154    Number epsilon_c_;
     155    /** Maximal number of second order correction steps */
     156    Index max_soc_;
     157    /** Required reduction in constraint violation before trying
     158     *  multiple second order correction steps \f$ \kappa_{soc}\f$.
     159     */
     160    Number kappa_soc_;
    155161    //@}
    156162
     
    175181     *  respect to which progress is to be made (at watchdog point) */
    176182    Number watchdog_direct_deriv_penalty_function_;
     183    /** Backup for the Chen-Goldfarb search direction (needed in the
     184     *  update rule for the penalty parameter */
     185    SmartPtr<const IteratesVector> watchdog_delta_cgpen_;
    177186    //@}
    178187
  • branches/dev-penalty/Algorithm/IpCGSearchDirCalc.cpp

    r550 r551  
    3737      "penalty_init_min",
    3838      "Minimal value for the intial penalty parameter (for Chen-Goldfarb line search).",
    39       0, true, 1e-3,
     39      0, true, 1.,
    4040      "");
    4141    roptions->AddLowerBoundedNumberOption(
    4242      "penalty_init_max",
    4343      "Maximal value for the intial penalty parameter (for Chen-Goldfarb line search).",
    44       0, true, 1.,
     44      0, true, 10.,
     45      "");
     46    roptions->AddStringOption2(
     47      "never_use_fact_cgpen_direction",
     48      "Toggle to switch off the fast Chen-Goldfarb direction",
     49      "no",
     50      "no", "always compute the fast direction",
     51      "yes", "never compute the fast direction",
    4552      "");
    4653  }
     
    5158    options.GetNumericValue("penalty_init_min", penalty_init_min_, prefix);
    5259    options.GetNumericValue("penalty_init_max", penalty_init_max_, prefix);
     60    options.GetBoolValue("never_use_fact_cgpen_direction",
     61                         never_use_fact_cgpen_direction_, prefix);
    5362
    5463    return pd_solver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
     
    7685
    7786    if (!IpData().PenaltyInitialized()) {
    78       Number penalty_init = Max(penalty_init_max_,
    79                                 Min(Max(IpData().curr()->y_c()->Amax(),
    80                                         IpData().curr()->y_d()->Amax()),
    81                                     penalty_init_min_));
     87      Number y_max = Max(IpData().curr()->y_c()->Amax(),
     88                         IpData().curr()->y_d()->Amax());
     89      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
     90                     "Initializing penalty parameter...\n");
     91      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
     92                     "Max(||y_c||_inf,||y_d||_inf = %8.2e\n",
     93                     y_max);
     94      Number penalty_init = Max(penalty_init_min_,
     95                                Min(y_max, penalty_init_max_));
    8296      IpData().Set_penalty(penalty_init);
    8397      char spen[40];
    8498      sprintf(spen, " penaltyinit=%8.2e", penalty_init);
    8599      IpData().Append_info_string(spen);
     100      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
     101                     "Initial value of the penalty parameter = %8.2e\n",
     102                     penalty_init);
    86103    }
    87104
     
    99116
    100117    // Get space for the search direction
    101     SmartPtr<IteratesVector> delta =
     118    SmartPtr<IteratesVector> delta_cgpen =
    102119      IpData().curr()->MakeNewIteratesVector(true);
    103120
    104121    bool allow_inexact = false;
    105     pd_solver_->Solve(-1.0, 0.0, *rhs, *delta, allow_inexact,
     122    pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_cgpen, allow_inexact,
    106123                      improve_solution);
    107124
    108125    // Store the original search direction in the IpData object
     126    IpData().set_delta_cgpen(delta_cgpen);
     127    IpData().SetHaveCgPenDeltas(true);
     128
     129    bool keep_fast_delta = true;
     130    if (never_use_fact_cgpen_direction_) {
     131      IpData().SetHaveCgFastDeltas(false);
     132      keep_fast_delta = false;
     133    }
     134    else {
     135      // Now solve it again but for the "fast" search direction
     136      rhs->Set_y_c(*IpCq().curr_c());
     137      rhs->Set_y_d(*IpCq().curr_d_minus_s());
     138
     139      // Get space for the search direction
     140      SmartPtr<IteratesVector> delta_fast =
     141        IpData().curr()->MakeNewIteratesVector(true);
     142
     143      allow_inexact = false;
     144      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_fast, allow_inexact,
     145                        improve_solution);
     146
     147      // Store the fast search direction in the IpData object
     148      IpData().set_delta_cgfast(delta_fast);
     149      IpData().SetHaveCgFastDeltas(true);
     150
     151      // Now we check whether the fast direction is good compatible with
     152      // the merit function
     153
     154      // do the || tilde y - hat y ||_2 <= k_dis ||hat y||_2 test
     155      SmartPtr<const Vector> y_c = IpData().curr()->y_c();
     156      SmartPtr<const Vector> y_d = IpData().curr()->y_d();
     157      SmartPtr<const Vector> delta_fast_y_c = IpData().delta_cgfast()->y_c();
     158      SmartPtr<const Vector> delta_fast_y_d = IpData().delta_cgfast()->y_d();
     159      SmartPtr<const Vector> delta_y_c = IpData().delta_cgpen()->y_c();
     160      SmartPtr<const Vector> delta_y_d = IpData().delta_cgpen()->y_d();
     161
     162      Number hat_y_nrm = sqrt(pow(y_c->Nrm2(), 2.)
     163                              + pow(y_d->Nrm2(), 2.)
     164                              + 2.*y_c->Dot(*delta_y_c)
     165                              + 2.*y_d->Dot(*delta_y_d)
     166                              + pow(delta_y_c->Nrm2(), 2.)
     167                              + pow(delta_y_d->Nrm2(), 2.));
     168      Number diff_y_nrm = sqrt(pow(delta_fast_y_c->Nrm2(), 2.)
     169                               + pow(delta_fast_y_d->Nrm2(), 2.)
     170                               - 2.*delta_y_c->Dot(*delta_fast_y_c)
     171                               - 2.*delta_y_d->Dot(*delta_fast_y_d)
     172                               + pow(delta_y_c->Nrm2(), 2.)
     173                               + pow(delta_y_d->Nrm2(), 2.));
     174
     175      Number kappa_dis_ = 1e1;
     176      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
     177                     "Testing if fast direction can be used.\n"
     178                     "  diff_y_nrm = %8.2e hat_y_norm = %8.2e\n",
     179                     diff_y_nrm, hat_y_nrm);
     180      if (diff_y_nrm > kappa_dis_*hat_y_nrm) {
     181        keep_fast_delta = false;
     182        IpData().Append_info_string("G");
     183      }
     184
     185      if (keep_fast_delta) {
     186        // For now, I just check if the directional derivative for the
     187        // penalty functions are not too much off
     188        Number direct_deriv = IpCq().curr_direct_deriv_penalty_function();
     189        Number fast_direct_deriv = IpCq().curr_fast_direct_deriv_penalty_function();
     190        Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
     191                       "direct_deriv = %23.15e  fast_direct_deriv = %23.15e\n",
     192                       direct_deriv, fast_direct_deriv);
     193        Number need_name_ = 1e-1;
     194        if (fast_direct_deriv > need_name_*direct_deriv) {
     195          // Discard the fast direction
     196          //delta_fast = NULL;
     197          //IpData().set_delta_cgpen(delta_fast);
     198          keep_fast_delta = false;
     199          IpData().Append_info_string("g");
     200        }
     201      }
     202    }
     203
     204    SmartPtr<const IteratesVector> delta;
     205    if (!keep_fast_delta) {
     206      IpData().SetHaveCgFastDeltas(false);
     207      delta = IpData().delta_cgpen();
     208    }
     209    else {
     210      delta = IpData().delta_cgfast();
     211    }
    109212    IpData().set_delta(delta);
    110213
    111     // Now solve it again but for the "fast" search direction
    112     rhs->Set_y_c(*IpCq().curr_c());
    113     rhs->Set_y_d(*IpCq().curr_d_minus_s());
    114 
    115     // Get space for the search direction
    116     SmartPtr<IteratesVector> delta_fast =
    117       IpData().curr()->MakeNewIteratesVector(true);
    118 
    119     allow_inexact = false;
    120     pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_fast, allow_inexact,
    121                       improve_solution);
    122 
    123     // Store the fast search direction in the IpData object
    124     IpData().set_delta_cgpen(delta_fast);
    125     IpData().SetHaveCgPenDeltas(true);
    126 
    127     // Now we check whether the fast direction is good compatible with
    128     // the merit function
    129 
    130     bool keep_fast_delta = true;
    131     // do the || tilde y - hat y ||_2 <= k_dis ||hat y||_2 test
    132     SmartPtr<const Vector> y_c = IpData().curr()->y_c();
    133     SmartPtr<const Vector> y_d = IpData().curr()->y_d();
    134     SmartPtr<const Vector> delta_fast_y_c = IpData().delta_cgpen()->y_c();
    135     SmartPtr<const Vector> delta_fast_y_d = IpData().delta_cgpen()->y_d();
    136     SmartPtr<const Vector> delta_y_c = IpData().delta()->y_c();
    137     SmartPtr<const Vector> delta_y_d = IpData().delta()->y_d();
    138 
    139     Number hat_y_nrm = sqrt(pow(y_c->Nrm2(), 2.) +
    140                             pow(y_d->Nrm2(), 2.) +
    141                             2.*y_c->Dot(*delta_y_c) +
    142                             2.*y_d->Dot(*delta_y_d) +
    143                             pow(delta_y_c->Nrm2(), 2.) +
    144                             pow(delta_y_d->Nrm2(), 2.));
    145     Number diff_y = sqrt(pow(delta_fast_y_c->Nrm2(), 2.) +
    146                          pow(delta_fast_y_d->Nrm2(), 2.) -
    147                          2.*delta_y_c->Dot(*delta_fast_y_c) -
    148                          2.*delta_y_d->Dot(*delta_fast_y_d) +
    149                          pow(delta_y_c->Nrm2(), 2.) +
    150                          pow(delta_y_d->Nrm2(), 2.));
    151 
    152     Number kappa_dis_ = 1e1;
    153     if (diff_y > kappa_dis_*hat_y_nrm) {
    154       keep_fast_delta = false;
    155       IpData().Append_info_string("G");
    156     }
    157 
    158     if (keep_fast_delta) {
    159       // For now, I just check if the directional derivative for the
    160       // penalty functions are not too much off
    161       Number direct_deriv = IpCq().curr_direct_deriv_penalty_function();
    162       Number fast_direct_deriv = IpCq().curr_fast_direct_deriv_penalty_function();
    163       Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
    164                      "direct_deriv = %23.16e  fast_direct_deriv = %23.15e\n",
    165                      direct_deriv, fast_direct_deriv);
    166       Number need_name_ = 1e-1;
    167       if (fast_direct_deriv > need_name_*direct_deriv) {
    168         // Discard the fast direction
    169         //delta_fast = NULL;
    170         //IpData().set_delta_cgpen(delta_fast);
    171         keep_fast_delta = false;
    172         IpData().Append_info_string("g");
    173       }
    174     }
    175 
    176     if (!keep_fast_delta) {
    177       IpData().SetHaveCgPenDeltas(false);
    178     }
    179 
    180214  }
    181215
  • branches/dev-penalty/Algorithm/IpCGSearchDirCalc.hpp

    r549 r551  
    7474    /** Maximal value for initial penalty parameter. */
    7575    Number penalty_init_max_;
     76    /** Flag indicating whether the fast Chen-Goldfarb direction
     77     *  should never be used */
     78    bool never_use_fact_cgpen_direction_;
    7679    //@}
    7780
  • branches/dev-penalty/Algorithm/IpFilterLSAcceptor.cpp

    r549 r551  
    424424    Number alpha_primal_soc = alpha_primal;
    425425
    426     SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNew();
    427     SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNew();
    428     c_soc->Copy(*IpCq().curr_c());
    429     dms_soc->Copy(*IpCq().curr_d_minus_s());
     426    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNewCopy();
     427    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNewCopy();
    430428    while (count_soc<max_soc_ && !accept &&
    431429           (count_soc==0 || theta_trial<=kappa_soc_*theta_soc_old) ) {
  • branches/dev-penalty/Algorithm/IpIpoptCalculatedQuantities.cpp

    r549 r551  
    33373337    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
    33383338    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
    3339     SmartPtr<const Vector> dy_c = ip_data_->delta()->y_c();
    3340     SmartPtr<const Vector> dy_d = ip_data_->delta()->y_d();
    3341     std::vector<const TaggedObject*> tdeps(6);
     3339    SmartPtr<const Vector> dy_c = ip_data_->delta_cgpen()->y_c();
     3340    SmartPtr<const Vector> dy_d = ip_data_->delta_cgpen()->y_d();
     3341    SmartPtr<const Vector> dx = ip_data_->delta_cgpen()->x();
     3342    SmartPtr<const Vector> ds = ip_data_->delta_cgpen()->s();
     3343    std::vector<const TaggedObject*> tdeps(8);
    33423344    tdeps[0] = GetRawPtr(x);
    33433345    tdeps[1] = GetRawPtr(s);
     
    33463348    tdeps[4] = GetRawPtr(dy_c);
    33473349    tdeps[5] = GetRawPtr(dy_d);
     3350    tdeps[6] = GetRawPtr(dx);
     3351    tdeps[7] = GetRawPtr(ds);
    33483352    Number mu = ip_data_->curr_mu();
    33493353    Number penalty = ip_data_->curr_penalty();
     
    33553359      SmartPtr<const Vector> c = curr_c();
    33563360      SmartPtr<const Vector> d_minus_s = curr_d_minus_s();
    3357       result = curr_gradBarrTDelta();
     3361      result = curr_grad_barrier_obj_x()->Dot(*dx) +
     3362               curr_grad_barrier_obj_s()->Dot(*ds);
    33583363      result -= penalty*curr_primal_infeasibility(NORM_2);
    33593364      result += c->Dot(*y_c);
     
    33763381    SmartPtr<const Vector> s = ip_data_->curr()->s();
    33773382    DBG_ASSERT(ip_data_->HaveCgPenDeltas());
    3378     SmartPtr<const Vector> dy_c = ip_data_->delta_cgpen()->y_c();
    3379     SmartPtr<const Vector> dy_d = ip_data_->delta_cgpen()->y_d();
    3380     std::vector<const TaggedObject*> tdeps(4);
     3383    SmartPtr<const Vector> dy_c = ip_data_->delta_cgfast()->y_c();
     3384    SmartPtr<const Vector> dy_d = ip_data_->delta_cgfast()->y_d();
     3385    SmartPtr<const Vector> dx = ip_data_->delta_cgfast()->x();
     3386    SmartPtr<const Vector> ds = ip_data_->delta_cgfast()->s();
     3387    std::vector<const TaggedObject*> tdeps(6);
    33813388    tdeps[0] = GetRawPtr(x);
    33823389    tdeps[1] = GetRawPtr(s);
    33833390    tdeps[2] = GetRawPtr(dy_c);
    33843391    tdeps[3] = GetRawPtr(dy_d);
     3392    tdeps[4] = GetRawPtr(dx);
     3393    tdeps[5] = GetRawPtr(ds);
    33853394    Number mu = ip_data_->curr_mu();
    33863395    Number penalty = ip_data_->curr_penalty();
     
    33923401      SmartPtr<const Vector> c = curr_c();
    33933402      SmartPtr<const Vector> d_minus_s = curr_d_minus_s();
    3394       result = curr_gradBarrTDelta();
     3403      result = curr_grad_barrier_obj_x()->Dot(*dx) +
     3404               curr_grad_barrier_obj_s()->Dot(*ds);
    33953405      result -= penalty*curr_primal_infeasibility(NORM_2);
    33963406      result += c->Dot(*dy_c);
  • branches/dev-penalty/Algorithm/IpIpoptData.cpp

    r549 r551  
    5757    have_affine_deltas_ = false;
    5858    have_cgpen_deltas_ = false;
     59    have_cgfast_deltas_ = false;
    5960
    6061    free_mu_mode_ = false;
     
    129130    debug_delta_aff_tag_sum_ = 0;
    130131    debug_delta_cgpen_tag_ = 0;
     132    debug_delta_cgfast_tag_ = 0;
    131133    debug_delta_cgpen_tag_sum_ = 0;
     134    debug_delta_cgfast_tag_sum_ = 0;
    132135#endif
    133136
     
    138141    delta_aff_ = NULL;
    139142    delta_cgpen_ = NULL;
     143    delta_cgfast_ = NULL;
    140144
    141145    have_prototypes_ = true;
     
    143147    have_affine_deltas_ = false;
    144148    have_cgpen_deltas_ = false;
     149    have_cgfast_deltas_ = false;
    145150
    146151    return true;
     
    231236    // Free the memory for the Chen-Goldfarb step
    232237    delta_cgpen_ = NULL;
     238    delta_cgfast_ = NULL;
    233239
    234240    have_deltas_ = false;
    235241    have_affine_deltas_ = false;
    236242    have_cgpen_deltas_ = false;
     243    have_cgfast_deltas_ = false;
    237244  }
    238245
  • branches/dev-penalty/Algorithm/IpIpoptData.hpp

    r550 r551  
    116116    SmartPtr<const IteratesVector> delta() const;
    117117
    118     /** Set the current delta - like the trial point, this method copies
    119      *  the pointer for efficiency (no copy and to keep cache tags the
    120      *  same) so after you call set, you cannot modify the data
     118    /** Set the current delta - like the trial point, this method
     119     *  copies the pointer for efficiency (no copy and to keep cache
     120     *  tags the same) so after you call set, you cannot modify the
     121     *  data.
    121122     */
    122123    void set_delta(SmartPtr<IteratesVector>& delta);
     
    126127     *  tags the same) so after you call set, you cannot modify the
    127128     *  data.  This is the version that is happy with a pointer to
    128      *  const IteratesVector
     129     *  const IteratesVector.
    129130     */
    130131    void set_delta(SmartPtr<const IteratesVector>& delta);
     
    135136    /** Set the affine delta - like the trial point, this method copies
    136137     *  the pointer for efficiency (no copy and to keep cache tags the
    137      *  same) so after you call set, you cannot modify the data
     138     *  same) so after you call set, you cannot modify the data.
    138139     */
    139140    void set_delta_aff(SmartPtr<IteratesVector>& delta_aff);
    140141
    141     /** Delta for the fast Chen-Goldfarb search direction */
     142    /** Delta for the Chen-Goldfarb search direction */
    142143    SmartPtr<const IteratesVector> delta_cgpen() const;
    143144
    144145    /** Set the delta_cgpen - like the trial point, this method copies
    145146     *  the pointer for efficiency (no copy and to keep cache tags the
    146      *  same) so after you call set, you cannot modify the data
     147     *  same) so after you call set, you cannot modify the data.
    147148     */
    148149    void set_delta_cgpen(SmartPtr<IteratesVector>& delta_pen);
    149150
    150     /** Hessian or Hessian approximation (do not hold on to it, it might be changed) */
     151    /** Set the delta_cgpen - like the trial point, this method copies
     152     *  the pointer for efficiency (no copy and to keep cache tags the
     153     *  same) so after you call set, you cannot modify the data.  This
     154     *  is the version that is happy with a pointer to const
     155     *  IteratesVector.
     156     */
     157    void set_delta_cgpen(SmartPtr<const IteratesVector>& delta_pen);
     158
     159    /** Delta for the fast Chen-Goldfarb search direction */
     160    SmartPtr<const IteratesVector> delta_cgfast() const;
     161
     162    /** Set the delta_cgpen - like the trial point, this method copies
     163     *  the pointer for efficiency (no copy and to keep cache tags the
     164     *  same) so after you call set, you cannot modify the data.
     165     */
     166    void set_delta_cgfast(SmartPtr<IteratesVector>& delta_fast);
     167
     168    /** Hessian or Hessian approximation (do not hold on to it, it
     169     *  might be changed) */
    151170    SmartPtr<const SymMatrix> W()
    152171    {
     
    231250    }
    232251
     252    bool HaveCgFastDeltas() const
     253    {
     254      return have_cgfast_deltas_;
     255    }
     256    void SetHaveCgFastDeltas(bool have_cgfast_deltas)
     257    {
     258      have_cgfast_deltas_ = have_cgfast_deltas;
     259    }
     260
    233261
    234262    /** @name Public Methods for updating iterates */
     
    473501     */
    474502    bool have_cgpen_deltas_;
     503    //@}
     504
     505    /** @name Fast Chen-Goldfarb step for the penatly function.  This
     506     *  used to transfer the information about the step from the
     507     *  computation of the overall search direction to the line
     508     *  search. */
     509    //@{
     510    SmartPtr<const IteratesVector> delta_cgfast_;
     511    /** The following flag is set to true, if some other part of the
     512     *  algorithm has already computed the fast Chen-Goldfarb step.
     513     *  This flag is reset when the AcceptTrialPoint method is called.
     514     *  * ToDo: we could cue off of a null delta_cgfast_;
     515     */
     516    bool have_cgfast_deltas_;
    475517    //@}
    476518
     
    572614    TaggedObject::Tag debug_delta_aff_tag_;
    573615    TaggedObject::Tag debug_delta_cgpen_tag_;
     616    TaggedObject::Tag debug_delta_cgfast_tag_;
    574617    TaggedObject::Tag debug_curr_tag_sum_;
    575618    TaggedObject::Tag debug_trial_tag_sum_;
     
    577620    TaggedObject::Tag debug_delta_aff_tag_sum_;
    578621    TaggedObject::Tag debug_delta_cgpen_tag_sum_;
     622    TaggedObject::Tag debug_delta_cgfast_tag_sum_;
    579623    //@}
    580624#endif
     
    630674
    631675    return delta_cgpen_;
     676  }
     677
     678  inline
     679  SmartPtr<const IteratesVector> IpoptData::delta_cgfast() const
     680  {
     681#   ifdef IP_DEBUG
     682    DBG_ASSERT(IsNull(delta_cgfast_) || (delta_cgfast_->GetTag() == debug_delta_cgfast_tag_ && delta_cgfast_->GetTagSum() == debug_delta_cgfast_tag_sum_) );
     683#   endif
     684
     685    return delta_cgfast_;
    632686  }
    633687
     
    747801  }
    748802
     803  inline
     804  void IpoptData::set_delta_cgpen(SmartPtr<const IteratesVector>& delta_cgpen)
     805  {
     806    delta_cgpen_ = delta_cgpen;
     807#ifdef IP_DEBUG
     808
     809    if (IsValid(delta_cgpen)) {
     810      debug_delta_cgpen_tag_ = delta_cgpen->GetTag();
     811      debug_delta_cgpen_tag_sum_ = delta_cgpen->GetTagSum();
     812    }
     813    else {
     814      debug_delta_cgpen_tag_ = 0;
     815      debug_delta_cgpen_tag_sum_ = delta_cgpen->GetTagSum();
     816    }
     817#endif
     818
     819    delta_cgpen = NULL;
     820  }
     821
     822  inline
     823  void IpoptData::set_delta_cgfast(SmartPtr<IteratesVector>& delta_cgfast)
     824  {
     825    delta_cgfast_ = ConstPtr(delta_cgfast);
     826#ifdef IP_DEBUG
     827
     828    if (IsValid(delta_cgfast)) {
     829      debug_delta_cgfast_tag_ = delta_cgfast->GetTag();
     830      debug_delta_cgfast_tag_sum_ = delta_cgfast->GetTagSum();
     831    }
     832    else {
     833      debug_delta_cgfast_tag_ = 0;
     834      debug_delta_cgfast_tag_sum_ = delta_cgfast->GetTagSum();
     835    }
     836#endif
     837
     838    delta_cgfast = NULL;
     839  }
     840
    749841} // namespace Ipopt
    750842
  • branches/dev-penalty/Algorithm/IpStdAugSystemSolver.cpp

    r517 r551  
    130130                            *J_c, D_c, delta_c, *J_d, D_d, delta_d,
    131131                            rhs_x, rhs_s, rhs_c, rhs_d);
     132      Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
     133                     "Augmented system has changed, need to refactorize.\n");
    132134    }
    133135
     
    179181                               check_NegEVals, numberOfNegEVals);
    180182    if (retval==SYMSOLVER_SUCCESS) {
    181       Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Factorization successful.\n");
     183      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Solution successful.\n");
    182184      augmented_sol->Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "SOL");
    183185    }
     
    186188    }
    187189    else {
    188       Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Factorization failed with retval = %d\n", retval);
     190      Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Solution failed with retval = %d\n", retval);
    189191    }
    190192
Note: See TracChangeset for help on using the changeset viewer.