source: branches/dev/Algorithm/IpIpoptAlg.cpp @ 416

Last change on this file since 416 was 416, checked in by andreasw, 14 years ago
  • revised handling of "acceptable level of accuracy" (now in ConvCheck?)
  • fixed uncaught evaluation error exceptions
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.4 KB
Line 
1// Copyright (C) 2004, International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpIpoptAlg.cpp 416 2005-07-29 19:11:41Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpIpoptAlg.hpp"
10#include "IpJournalist.hpp"
11#include "IpRestoPhase.hpp"
12
13namespace Ipopt
14{
15  DBG_SET_VERBOSITY(0);
16
17  DefineIpoptType(IpoptAlgorithm);
18
19  IpoptAlgorithm::IpoptAlgorithm(const SmartPtr<PDSystemSolver>& pd_solver,
20                                 const SmartPtr<LineSearch>& line_search,
21                                 const SmartPtr<MuUpdate>& mu_update,
22                                 const SmartPtr<ConvergenceCheck>& conv_check,
23                                 const SmartPtr<IterateInitializer>& iterate_initializer,
24                                 const SmartPtr<IterationOutput>& iter_output)
25      :
26      pd_solver_(pd_solver),
27      line_search_(line_search),
28      mu_update_(mu_update),
29      conv_check_(conv_check),
30      iterate_initializer_(iterate_initializer),
31      iter_output_(iter_output)
32  {
33    DBG_START_METH("IpoptAlgorithm::IpoptAlgorithm",
34                   dbg_verbosity);
35    DBG_ASSERT(IsValid(pd_solver_));
36    DBG_ASSERT(IsValid(line_search_));
37    DBG_ASSERT(IsValid(mu_update_));
38    DBG_ASSERT(IsValid(conv_check_));
39    DBG_ASSERT(IsValid(iterate_initializer_));
40    DBG_ASSERT(IsValid(iter_output_));
41  }
42
43  IpoptAlgorithm::~IpoptAlgorithm()
44  {
45    DBG_START_METH("IpoptAlgorithm::~IpoptAlgorithm()",
46                   dbg_verbosity);
47  }
48
49  void IpoptAlgorithm::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
50  {
51    roptions->AddLowerBoundedNumberOption(
52      "kappa_sigma",
53      "Factor limiting deviation of dual variables from primal estimates.",
54      0, true, 1e10,
55      "(see Eqn. (16) in the implementation paper.) Setting to value less than "
56      "one disables correction.");
57  }
58
59  bool IpoptAlgorithm::InitializeImpl(const OptionsList& options,
60                                      const std::string& prefix)
61  {
62    DBG_START_METH("IpoptAlgorithm::InitializeImpl",
63                   dbg_verbosity);
64
65    // Read the IpoptAlgorithm options
66    // Initialize the Data object
67    bool retvalue = IpData().Initialize(Jnlst(),
68                                        options, prefix);
69    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
70                     "the IpIpoptData object failed to initialize.");
71
72    // Initialize the CQ object
73    retvalue = IpCq().Initialize(Jnlst(),
74                                 options, prefix);
75    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
76                     "the IpIpoptCalculatedQuantities object failed to initialize.");
77
78    // Initialize the CQ object
79    retvalue = IpNLP().Initialize(Jnlst(),
80                                  options, prefix);
81    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
82                     "the IpIpoptNLP object failed to initialize.");
83
84    // Initialize all the strategies
85    retvalue = iterate_initializer_->Initialize(Jnlst(), IpNLP(), IpData(),
86               IpCq(), options, prefix);
87    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
88                     "the iterate_initializer strategy failed to initialize.");
89
90    retvalue = mu_update_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
91                                      options, prefix);
92    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
93                     "the mu_update strategy failed to initialize.");
94
95    retvalue = pd_solver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
96                                      options, prefix);
97    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
98                     "the pd_solver strategy failed to initialize.");
99
100    retvalue = line_search_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
101                                        options,prefix);
102    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
103                     "the line_search strategy failed to initialize.");
104
105    retvalue = conv_check_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
106                                       options, prefix);
107    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
108                     "the conv_check strategy failed to initialize.");
109
110    retvalue = iter_output_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
111                                        options, prefix);
112    ASSERT_EXCEPTION(retvalue, FAILED_INITIALIZATION,
113                     "the iter_output strategy failed to initialize.");
114
115    options.GetNumericValue("kappa_sigma", kappa_sigma_, prefix);
116
117    if (prefix=="resto.") {
118      skip_print_problem_stats_ = true;
119    }
120    else {
121      skip_print_problem_stats_ = false;
122    }
123
124    return true;
125  }
126
127  SolverReturn IpoptAlgorithm::Optimize()
128  {
129    DBG_START_METH("IpoptAlgorithm::Optimize", dbg_verbosity);
130
131    // Initialize the iterates
132    InitializeIterates();
133
134    if (!skip_print_problem_stats_) {
135      PrintProblemStatistics();
136    }
137
138    ConvergenceCheck::ConvergenceStatus conv_status
139    = conv_check_->CheckConvergence();
140
141    try {
142      // main loop
143      while (conv_status == ConvergenceCheck::CONTINUE) {
144        // Set the Hessian Matrix
145        ActualizeHessian();
146
147        // do all the output for this iteration
148        OutputIteration();
149        IpData().ResetInfo();
150
151        // update the barrier parameter if necessary
152        UpdateBarrierParameter();
153
154        // solve the primal-dual system to get the full step
155        ComputeSearchDirection();
156
157        // Compute the new iterate
158        ComputeAcceptableTrialPoint();
159        //ApplyFractionToBoundary();
160
161        // Accept the new iterate
162        AcceptTrialPoint();
163
164        IpData().Set_iter_count(IpData().iter_count()+1);
165
166        conv_status  = conv_check_->CheckConvergence();
167      }
168
169      OutputIteration();
170
171      if (conv_status == ConvergenceCheck::CONVERGED) {
172        return SUCCESS;
173      }
174      if (conv_status == ConvergenceCheck::CONVERGED_TO_ACCEPTABLE_POINT) {
175        return STOP_AT_ACCEPTABLE_POINT;
176      }
177      else if (conv_status == ConvergenceCheck::MAXITER_EXCEEDED) {
178        return MAXITER_EXCEEDED;
179      }
180    }
181    catch(TINY_STEP_DETECTED& exc) {
182      exc.ReportException(Jnlst());
183      return STOP_AT_TINY_STEP;
184    }
185    catch(ACCEPTABLE_POINT_REACHED& exc) {
186      exc.ReportException(Jnlst());
187      return STOP_AT_ACCEPTABLE_POINT;
188    }
189    catch(LOCALLY_INFEASIBLE& exc) {
190      exc.ReportException(Jnlst());
191      return LOCAL_INFEASIBILITY;
192    }
193    catch(RESTORATION_FAILED& exc) {
194      exc.ReportException(Jnlst());
195      return RESTORATION_FAILURE;
196    }
197    catch(FEASIBILITY_PROBLEM_SOLVED& exc) {
198      exc.ReportException(Jnlst());
199      return SUCCESS;
200    }
201
202    DBG_ASSERT(false && "Unknown return code in the algorithm");
203  }
204
205  void IpoptAlgorithm::ActualizeHessian()
206  {
207    // At this point, just compute the exact Hessian
208    IpData().Set_W(IpCq().curr_exact_hessian());
209  }
210
211
212  void IpoptAlgorithm::UpdateBarrierParameter()
213  {
214    Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n");
215    Jnlst().Printf(J_DETAILED, J_MAIN, "*** Update Barrier Parameter for Iteration %d:", IpData().iter_count());
216    Jnlst().Printf(J_DETAILED, J_MAIN, "\n**************************************************\n\n");
217    mu_update_->UpdateBarrierParameter();
218
219    Jnlst().Printf(J_DETAILED, J_MAIN, "Barrier Parameter: %e\n", IpData().curr_mu());
220
221  }
222
223  void IpoptAlgorithm::ComputeSearchDirection()
224  {
225    DBG_START_METH("IpoptAlgorithm::ComputeSearchDirection", dbg_verbosity);
226
227    Jnlst().Printf(J_DETAILED, J_MAIN,
228                   "\n**************************************************\n");
229    Jnlst().Printf(J_DETAILED, J_MAIN,
230                   "*** Solving the Primal Dual System for Iteration %d:",
231                   IpData().iter_count());
232    Jnlst().Printf(J_DETAILED, J_MAIN,
233                   "\n**************************************************\n\n");
234
235    if (IpData().HaveDeltas()) {
236      Jnlst().Printf(J_DETAILED, J_MAIN,
237                     "No need to compute search direction - it has already been computed.\n");
238      return;
239    }
240
241    SmartPtr<IteratesVector> rhs = IpData().curr()->MakeNewContainer();
242    rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
243    rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
244    rhs->Set_y_c(*IpCq().curr_c());
245    rhs->Set_y_d(*IpCq().curr_d_minus_s());
246    rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
247    rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
248    rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
249    rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
250
251    DBG_PRINT_VECTOR(2, "rhs", *rhs);
252
253    //ToDo: allow us to delete entries in IpData to save memory?
254    // To save memory, delete the old search directions
255    //     IpData().SetFromPtr_delta_x(NULL);
256    //     IpData().SetFromPtr_delta_s(NULL);
257    //     IpData().SetFromPtr_delta_y_c(NULL);
258    //     IpData().SetFromPtr_delta_y_d(NULL);
259    //     IpData().SetFromPtr_delta_z_L(NULL);
260    //     IpData().SetFromPtr_delta_z_U(NULL);
261    //     IpData().SetFromPtr_delta_v_L(NULL);
262    //     IpData().SetFromPtr_delta_v_U(NULL);
263
264    // Get space for the search direction
265    SmartPtr<IteratesVector> delta = IpData().curr()->MakeNewIteratesVector(true);
266
267    pd_solver_->Solve(-1.0, 0.0, *rhs, *delta);
268
269    // Store the search directions in the IpData object
270    IpData().set_delta(delta);
271
272    Jnlst().Printf(J_MOREVECTOR, J_MAIN,
273                   "*** Step Calculated for Iteration: %d\n",
274                   IpData().iter_count());
275    Jnlst().PrintVector(J_MOREVECTOR, J_MAIN, "delta", *IpData().delta());
276  }
277
278  void IpoptAlgorithm::ComputeAcceptableTrialPoint()
279  {
280    Jnlst().Printf(J_DETAILED, J_MAIN,
281                   "\n**************************************************\n");
282    Jnlst().Printf(J_DETAILED, J_MAIN,
283                   "*** Finding Acceptable Trial Point for Iteration %d:",
284                   IpData().iter_count());
285    Jnlst().Printf(J_DETAILED, J_MAIN,
286                   "\n**************************************************\n\n");
287    line_search_->FindAcceptableTrialPoint();
288  }
289
290  void IpoptAlgorithm::OutputIteration()
291  {
292    iter_output_->WriteOutput();
293  }
294
295  void IpoptAlgorithm::InitializeIterates()
296  {
297    DBG_START_METH("IpoptAlgorithm::InitializeIterates", dbg_verbosity);
298
299    iterate_initializer_->SetInitialIterates();
300  }
301
302  void IpoptAlgorithm::AcceptTrialPoint()
303  {
304    DBG_START_METH("IpoptAlgorithm::AcceptTrialPoint", dbg_verbosity);
305    // If the line search didn't determine a new acceptable trial
306    // point, do not accept a new iterate
307    if (line_search_->CheckSkippedLineSearch()) {
308      Jnlst().Printf(J_SUMMARY, J_MAIN,
309                     "Line search didn't find acceptable trial point.\n");
310      return;
311    }
312
313    // Adjust the bounds if necessary
314    Index adjusted_slacks = IpCq().AdjustedTrialSlacks();
315    DBG_PRINT((1, "adjusted_slacks = %d\n", adjusted_slacks));
316    if (adjusted_slacks>0) {
317      IpCq().ResetAdjustedTrialSlacks();
318      if (adjusted_slacks==1) {
319        Jnlst().Printf(J_SUMMARY, J_MAIN,
320                       "%d Slack too small, adjusting variable bound\n",
321                       adjusted_slacks);
322      }
323      else {
324        Jnlst().Printf(J_SUMMARY, J_MAIN,
325                       "%d Slacks too small, adjusting variable bounds\n",
326                       adjusted_slacks);
327      }
328      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
329        Jnlst().PrintVector(J_VECTOR, J_MAIN, "old_x_L", *IpNLP().x_L());
330        Jnlst().PrintVector(J_VECTOR, J_MAIN, "old_x_U", *IpNLP().x_U());
331        Jnlst().PrintVector(J_VECTOR, J_MAIN, "old_d_L", *IpNLP().d_L());
332        Jnlst().PrintVector(J_VECTOR, J_MAIN, "old_d_U", *IpNLP().d_U());
333      }
334
335      SmartPtr<Vector> new_x_l = IpNLP().x_L()->MakeNew();
336      IpNLP().Px_L()->TransMultVector(1.0, *IpData().trial()->x(),
337                                      0.0, *new_x_l);
338      new_x_l->Axpy(-1.0, *IpCq().trial_slack_x_L());
339
340      SmartPtr<Vector> new_x_u = IpNLP().x_U()->MakeNew();
341      IpNLP().Px_U()->TransMultVector(1.0, *IpData().trial()->x(),
342                                      0.0, *new_x_u);
343      new_x_u->Axpy(1.0, *IpCq().trial_slack_x_U());
344
345      SmartPtr<Vector> new_d_l = IpNLP().d_L()->MakeNew();
346      IpNLP().Pd_L()->TransMultVector(1.0, *IpData().trial()->s(),
347                                      0.0, *new_d_l);
348      new_d_l->Axpy(-1.0, *IpCq().trial_slack_s_L());
349
350      SmartPtr<Vector> new_d_u = IpNLP().d_U()->MakeNew();
351      IpNLP().Pd_U()->TransMultVector(1.0, *IpData().trial()->s(),
352                                      0.0, *new_d_u);
353      new_d_u->Axpy(1.0, *IpCq().trial_slack_s_U());
354
355      IpNLP().AdjustVariableBounds(*new_x_l, *new_x_u, *new_d_l, *new_d_u);
356
357      if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
358        Jnlst().PrintVector(J_VECTOR, J_MAIN, "new_x_L", *IpNLP().x_L());
359        Jnlst().PrintVector(J_VECTOR, J_MAIN, "new_x_U", *IpNLP().x_U());
360        Jnlst().PrintVector(J_VECTOR, J_MAIN, "new_d_L", *IpNLP().d_L());
361        Jnlst().PrintVector(J_VECTOR, J_MAIN, "new_d_U", *IpNLP().d_U());
362      }
363
364    }
365
366    // Make sure that bound multipliers are not too far from \mu * S^{-1}
367    // (see kappa_sigma in paper)
368    bool corrected = false;
369    Number max_correction;
370    SmartPtr<const Vector> new_z_L;
371    max_correction = correct_bound_multiplier(
372                       *IpData().trial()->z_L(),
373                       *IpCq().trial_slack_x_L(),
374                       *IpCq().trial_compl_x_L(),
375                       new_z_L);
376    if (max_correction>0.) {
377      Jnlst().Printf(J_DETAILED, J_MAIN,
378                     "Some value in z_L becomes too large - maximal correction = %8.2e\n",
379                     max_correction);
380      corrected = true;
381    }
382    SmartPtr<const Vector> new_z_U;
383    max_correction = correct_bound_multiplier(
384                       *IpData().trial()->z_U(),
385                       *IpCq().trial_slack_x_U(),
386                       *IpCq().trial_compl_x_U(),
387                       new_z_U);
388    if (max_correction>0.) {
389      Jnlst().Printf(J_DETAILED, J_MAIN,
390                     "Some value in z_U becomes too large - maximal correction = %8.2e\n",
391                     max_correction);
392      corrected = true;
393    }
394    SmartPtr<const Vector> new_v_L;
395    max_correction = correct_bound_multiplier(
396                       *IpData().trial()->v_L(),
397                       *IpCq().trial_slack_s_L(),
398                       *IpCq().trial_compl_s_L(),
399                       new_v_L);
400    if (max_correction>0.) {
401      Jnlst().Printf(J_DETAILED, J_MAIN,
402                     "Some value in v_L becomes too large - maximal correction = %8.2e\n",
403                     max_correction);
404      corrected = true;
405    }
406    SmartPtr<const Vector> new_v_U;
407    max_correction = correct_bound_multiplier(
408                       *IpData().trial()->v_U(),
409                       *IpCq().trial_slack_s_U(),
410                       *IpCq().trial_compl_s_U(),
411                       new_v_U);
412    if (max_correction>0.) {
413      Jnlst().Printf(J_DETAILED, J_MAIN,
414                     "Some value in v_U becomes too large - maximal correction = %8.2e\n",
415                     max_correction);
416      corrected = true;
417    }
418    SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
419    trial->Set_bound_mult(*new_z_L, *new_z_U, *new_v_L, *new_v_U);
420    IpData().set_trial(trial);
421
422    if (corrected) {
423      IpData().Append_info_string("z");
424    }
425
426    // Accept the step
427    IpData().AcceptTrialPoint();
428  }
429
430  void IpoptAlgorithm::PrintProblemStatistics()
431  {
432    if (!Jnlst().ProduceOutput(J_SUMMARY, J_STATISTICS)) {
433      // nothing to print
434      return;
435    }
436
437    SmartPtr<const Vector> x = IpData().curr()->x();
438    SmartPtr<const Vector> x_L = IpNLP().x_L();
439    SmartPtr<const Vector> x_U = IpNLP().x_U();
440    SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
441    SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
442
443    Index nx_tot, nx_only_lower, nx_both, nx_only_upper;
444    calc_number_of_bounds(*IpData().curr()->x(), *IpNLP().x_L(), *IpNLP().x_U(),
445                          *IpNLP().Px_L(), *IpNLP().Px_U(),
446                          nx_tot, nx_only_lower, nx_both, nx_only_upper);
447
448    Index ns_tot, ns_only_lower, ns_both, ns_only_upper;
449    calc_number_of_bounds(*IpData().curr()->s(), *IpNLP().d_L(), *IpNLP().d_U(),
450                          *IpNLP().Pd_L(), *IpNLP().Pd_U(),
451                          ns_tot, ns_only_lower, ns_both, ns_only_upper);
452
453    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
454                   "Total number of variables............................: %8d\n",nx_tot);
455    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
456                   "                     variables with only lower bounds: %8d\n",
457                   nx_only_lower);
458    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
459                   "                variables with lower and upper bounds: %8d\n",nx_both);
460    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
461                   "                     variables with only upper bounds: %8d\n",
462                   nx_only_upper);
463    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
464                   "Total number of equality constraints.................: %8d\n",
465                   IpData().curr()->y_c()->Dim());
466    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
467                   "Total number of inequality constraints...............: %8d\n",ns_tot);
468    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
469                   "        inequality constraints with only lower bounds: %8d\n",
470                   ns_only_lower);
471    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
472                   "   inequality constraints with lower and upper bounds: %8d\n",ns_both);
473    Jnlst().Printf(J_SUMMARY, J_STATISTICS,
474                   "        inequality constraints with only upper bounds: %8d\n\n",
475                   ns_only_upper);
476  }
477
478  void IpoptAlgorithm::calc_number_of_bounds(
479    const Vector& x,
480    const Vector& x_L,
481    const Vector& x_U,
482    const Matrix& Px_L,
483    const Matrix& Px_U,
484    Index& n_tot,
485    Index& n_only_lower,
486    Index& n_both,
487    Index& n_only_upper)
488  {
489    DBG_START_METH("IpoptAlgorithm::calc_number_of_bounds",
490                   dbg_verbosity);
491
492    n_tot = x.Dim();
493
494    SmartPtr<Vector> tmpx = x.MakeNew();
495    SmartPtr<Vector> tmpxL = x_L.MakeNew();
496    SmartPtr<Vector> tmpxU = x_U.MakeNew();
497
498    tmpxL->Set(-1.);
499    tmpxU->Set(2.);
500    Px_L.MultVector(1.0, *tmpxL, 0.0, *tmpx);
501    Px_U.MultVector(1.0, *tmpxU, 1.0, *tmpx);
502    // Now, x has elements
503    //  -1 : if component has only lower bound
504    //   0 : if component has no bound
505    //   1 : if component has both lower and upper bound
506    //   2 : if component has only upper bound
507    DBG_PRINT_VECTOR(2, "x-indicator", *tmpx);
508
509    SmartPtr<Vector> tmpx0 = x.MakeNew();
510    tmpx0->Set(0.);
511
512    SmartPtr<Vector> tmpx2 = x.MakeNew();
513    tmpx2->Set(-1.0);
514    tmpx2->Axpy(1.0, *tmpx);
515    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
516    // components with only upper bounds
517    n_only_upper = (Index)tmpx2->Asum();
518
519    tmpx->Axpy(-2., *tmpx2);       // now make all those entries for
520    // only upper bounds zero in tmpx
521
522    tmpx2->Copy(*tmpx);
523    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
524    // components with both bounds
525    n_both = (Index)tmpx2->Asum();
526
527    tmpx->Axpy(-1., *tmpx2);
528    tmpx->ElementWiseMin(*tmpx);   // tmpx is now -1 in those with only
529    // lower bounds
530    n_only_lower = (Index)tmpx->Asum();
531
532  }
533
534  Number IpoptAlgorithm::correct_bound_multiplier(
535    const Vector& trial_z,
536    const Vector& trial_slack,
537    const Vector& trial_compl,
538    SmartPtr<const Vector>& new_trial_z)
539  {
540    DBG_START_METH("IpoptAlgorithm::CorrectBoundMultiplier",
541                   dbg_verbosity);
542
543    if (kappa_sigma_<1. || trial_z.Dim()==0) {
544      new_trial_z = &trial_z;
545      return 0.;
546    }
547
548    // We choose as barrier parameter to be used either the current
549    // algorithmic barrier parameter (if we are not in the free mode),
550    // or the average complementarity (at the trial point)
551    Number mu;
552    if (IpData().FreeMuMode()) {
553      mu = IpCq().trial_avrg_compl();
554      mu = Min(mu, 1e3);
555    }
556    else {
557      mu = IpData().curr_mu();
558    }
559    DBG_PRINT((1,"mu = %8.2e\n", mu));
560    DBG_PRINT_VECTOR(2, "trial_z", trial_z);
561
562    // First check quickly if anything need to be corrected, using the
563    // trial complementarity directly.  Here, Amax is the same as Max
564    // (and we use Amax because that can be used later)
565    if (trial_compl.Amax() <= kappa_sigma_*mu &&
566        trial_compl.Min() >= 1./kappa_sigma_*mu) {
567      new_trial_z = &trial_z;
568      return 0.;
569    }
570
571    SmartPtr<Vector> one_over_s = trial_z.MakeNew();
572    one_over_s->Copy(trial_slack);
573    one_over_s->ElementWiseReciprocal();
574
575    SmartPtr<Vector> step_z = trial_z.MakeNew();
576    step_z->AddTwoVectors(kappa_sigma_*mu, *one_over_s, -1., trial_z, 0.);
577
578    DBG_PRINT_VECTOR(2, "step_z", *step_z);
579
580    Number max_correction_up = Max(0., -step_z->Min());
581    if (max_correction_up>0.) {
582      SmartPtr<Vector> tmp = trial_z.MakeNew();
583      tmp->Set(0.);
584      step_z->ElementWiseMin(*tmp);
585      tmp->AddTwoVectors(1., trial_z, 1., *step_z, 0.);
586      new_trial_z = GetRawPtr(tmp);
587    }
588    else {
589      new_trial_z = &trial_z;
590    }
591
592    step_z->AddTwoVectors(1./kappa_sigma_*mu, *one_over_s, -1., *new_trial_z, 0.);
593
594    Number max_correction_low = Max(0., step_z->Max());
595    if (max_correction_low>0.) {
596      SmartPtr<Vector> tmp = trial_z.MakeNew();
597      tmp->Set(0.);
598      step_z->ElementWiseMax(*tmp);
599      tmp->AddTwoVectors(1., *new_trial_z, 1., *step_z, 0.);
600      new_trial_z = GetRawPtr(tmp);
601    }
602
603    DBG_PRINT_VECTOR(2, "new_trial_z", *new_trial_z);
604
605    return Max(max_correction_up, max_correction_low);
606  }
607
608} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.