source: branches/dev/Algorithm/IpQualityFunctionMuOracle.cpp @ 567

Last change on this file since 567 was 567, checked in by andreasw, 15 years ago

adapted for Portland compilers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 50.7 KB
Line 
1// Copyright (C) 2004, 2005 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: IpQualityFunctionMuOracle.cpp 567 2005-11-01 23:27:11Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter            IBM    2004-11-12
8
9#include "IpQualityFunctionMuOracle.hpp"
10
11#ifdef HAVE_CMATH
12# include <cmath>
13#else
14# ifdef HAVE_MATH_H
15#  include <math.h>
16# else
17#  error "don't have header file for math"
18# endif
19#endif
20
21namespace Ipopt
22{
23#ifdef IP_DEBUG
24  static const Index dbg_verbosity = 0;
25#endif
26
27  QualityFunctionMuOracle::QualityFunctionMuOracle(const SmartPtr<PDSystemSolver>& pd_solver)
28      :
29      MuOracle(),
30      pd_solver_(pd_solver),
31
32      tmp_step_x_L_(NULL),
33      tmp_step_x_U_(NULL),
34      tmp_step_s_L_(NULL),
35      tmp_step_s_U_(NULL),
36      tmp_step_z_L_(NULL),
37      tmp_step_z_U_(NULL),
38      tmp_step_v_L_(NULL),
39      tmp_step_v_U_(NULL),
40
41      tmp_slack_x_L_(NULL),
42      tmp_slack_x_U_(NULL),
43      tmp_slack_s_L_(NULL),
44      tmp_slack_s_U_(NULL),
45      tmp_z_L_(NULL),
46      tmp_z_U_(NULL),
47      tmp_v_L_(NULL),
48      tmp_v_U_(NULL),
49
50      count_qf_evals_(0)
51  {
52    DBG_ASSERT(IsValid(pd_solver_));
53  }
54
55  QualityFunctionMuOracle::~QualityFunctionMuOracle()
56  {}
57
58  void QualityFunctionMuOracle::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
59  {
60    roptions->AddLowerBoundedNumberOption(
61      "sigma_max",
62      "Maximum value of the centering parameter.",
63      0.0, true, 1e2);
64    roptions->AddLowerBoundedNumberOption(
65      "sigma_min",
66      "Minimum value of the centering parameter.",
67      0.0, false, 1e-6);
68    roptions->AddStringOption4(
69      "quality_function_norm_type",
70      "Norm used for components of the quality function.",
71      "2-norm-squared",
72      "1-norm", "use the 1-norm (abs sum)",
73      "2-norm-squared", "use the 2-norm squared (sum of squares)",
74      "max-norm", "use the infinity norm (max)",
75      "2-norm", "use 2-norm");
76    roptions->AddStringOption4(
77      "quality_function_centrality",
78      "The penalty term for centrality that is included in quality function.",
79      "none",
80      "none", "no penalty term is added",
81      "log", "complementarity * the log of the centrality measure",
82      "reciprocal", "complementarity * the reciprocal of the centrality measure",
83      "cubed-reciprocal", "complementarity * the reciprocal of the centrality measure cubed",
84      "This determines whether a term is added to the quality function to "
85      "penalize deviation from centrality with respect to complementarity.  The "
86      "complementarity measure here is the xi in the Loqo update rule.");
87    roptions->AddStringOption2(
88      "quality_function_balancing_term",
89      "The balancing term included in the quality function for centrality.",
90      "none",
91      "none", "no balancing term is added",
92      "cubic", "Max(0,Max(dual_inf,primal_inf)-compl)^3",
93      "This determines whether a term is added to the quality function that "
94      "penalizes situations where the complementarity is much smaller "
95      "than dual and primal infeasibilities.");
96    roptions->AddLowerBoundedIntegerOption(
97      "quality_function_max_section_steps",
98      "Maximum number of search steps during direct search procedure "
99      "determining the optimal centering parameter.",
100      0, 8,
101      "The golden section search is performed for the quality function based "
102      "mu oractle.");
103    roptions->AddBoundedNumberOption(
104      "quality_function_section_sigma_tol",
105      "Tolerance for the section search procedure determining "
106      "the optimal centering parameter (in sigma space).",
107      0.0, false, 1.0, true,
108      1e-2,
109      "The golden section search is performed for the quality function based "
110      "mu oractle.");
111    roptions->AddBoundedNumberOption(
112      "quality_function_section_qf_tol",
113      "Tolerance for the golden section search procedure determining "
114      "the optimal centering parameter (in the function value space).",
115      0.0, false, 1.0, true,
116      0e-2,
117      "The section search is performed for the quality function based mu "
118      "oractle.");
119  }
120
121
122  bool QualityFunctionMuOracle::InitializeImpl(const OptionsList& options,
123      const std::string& prefix)
124  {
125    Index enum_int;
126
127    options.GetNumericValue("sigma_max", sigma_max_, prefix);
128    options.GetNumericValue("sigma_min", sigma_min_, prefix);
129
130    options.GetEnumValue("quality_function_norm_type", enum_int, prefix);
131    quality_function_norm_ = NormEnum(enum_int);
132    options.GetEnumValue("quality_function_centrality", enum_int, prefix);
133    quality_function_centrality_ = CentralityEnum(enum_int);
134    options.GetEnumValue("quality_function_balancing_term", enum_int, prefix);
135    quality_function_balancing_term_ = BalancingTermEnum(enum_int);
136    options.GetIntegerValue("quality_function_max_section_steps",
137                            quality_function_max_section_steps_, prefix);
138    options.GetNumericValue("quality_function_section_sigma_tol",
139                            quality_function_section_sigma_tol_, prefix);
140    options.GetNumericValue("quality_function_section_qf_tol",
141                            quality_function_section_qf_tol_, prefix);
142
143    initialized_ = false;
144
145    return true;
146  }
147
148  Number QualityFunctionMuOracle::CalculateMu(Number mu_min, Number mu_max)
149  {
150    DBG_START_METH("QualityFunctionMuOracle::CalculateMu",
151                   dbg_verbosity);
152
153    ///////////////////////////////////////////////////////////////////////////
154    // Reserve memory for temporary vectors used in CalculateQualityFunction //
155    ///////////////////////////////////////////////////////////////////////////
156
157    tmp_step_x_L_ = IpNLP().x_L()->MakeNew();
158    tmp_step_x_U_ = IpNLP().x_U()->MakeNew();
159    tmp_step_s_L_ = IpNLP().d_L()->MakeNew();
160    tmp_step_s_U_ = IpNLP().d_U()->MakeNew();
161    tmp_step_z_L_ = IpNLP().x_L()->MakeNew();
162    tmp_step_z_U_ = IpNLP().x_U()->MakeNew();
163    tmp_step_v_L_ = IpNLP().d_L()->MakeNew();
164    tmp_step_v_U_ = IpNLP().d_U()->MakeNew();
165
166    tmp_slack_x_L_ = IpNLP().x_L()->MakeNew();
167    tmp_slack_x_U_ = IpNLP().x_U()->MakeNew();
168    tmp_slack_s_L_ = IpNLP().d_L()->MakeNew();
169    tmp_slack_s_U_ = IpNLP().d_U()->MakeNew();
170    tmp_z_L_ = IpNLP().x_L()->MakeNew();
171    tmp_z_U_ = IpNLP().x_U()->MakeNew();
172    tmp_v_L_ = IpNLP().d_L()->MakeNew();
173    tmp_v_U_ = IpNLP().d_U()->MakeNew();
174
175    /////////////////////////////////////
176    // Compute the affine scaling step //
177    /////////////////////////////////////
178
179    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
180                   "Solving the Primal Dual System for the affine step\n");
181    // First get the right hand side
182    SmartPtr<IteratesVector> rhs_aff = IpData().curr()->MakeNewIteratesVector(false);
183    rhs_aff->Set_x(*IpCq().curr_grad_lag_x());
184    rhs_aff->Set_s(*IpCq().curr_grad_lag_s());
185    rhs_aff->Set_y_c(*IpCq().curr_c());
186    rhs_aff->Set_y_d(*IpCq().curr_d_minus_s());
187    rhs_aff->Set_z_L(*IpCq().curr_compl_x_L());
188    rhs_aff->Set_z_U(*IpCq().curr_compl_x_U());
189    rhs_aff->Set_v_L(*IpCq().curr_compl_s_L());
190    rhs_aff->Set_v_U(*IpCq().curr_compl_s_U());
191
192    // Get space for the affine scaling step
193    SmartPtr<IteratesVector> step_aff = IpData().curr()->MakeNewIteratesVector(true);
194
195    // Now solve the primal-dual system to get the step.  We allow a
196    // somewhat inexact solution, iterative refinement will be done
197    // after mu is known
198    bool allow_inexact = true;
199    pd_solver_->Solve(-1.0, 0.0,
200                      *rhs_aff,
201                      *step_aff,
202                      allow_inexact
203                     );
204
205    DBG_PRINT_VECTOR(2, "step_aff", *step_aff);
206
207    /////////////////////////////////////
208    // Compute the pure centering step //
209    /////////////////////////////////////
210
211    Number avrg_compl = IpCq().curr_avrg_compl();
212
213    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
214                   "Solving the Primal Dual System for the centering step\n");
215    // First get the right hand side
216    SmartPtr<IteratesVector> rhs_cen = IpData().curr()->MakeNewIteratesVector(true);
217    rhs_cen->x_NonConst()->AddOneVector(-avrg_compl,
218                                        *IpCq().grad_kappa_times_damping_x(),
219                                        0.);
220    rhs_cen->s_NonConst()->AddOneVector(-avrg_compl,
221                                        *IpCq().grad_kappa_times_damping_s(),
222                                        0.);
223
224    rhs_cen->y_c_NonConst()->Set(0.);
225    rhs_cen->y_d_NonConst()->Set(0.);
226    rhs_cen->z_L_NonConst()->Set(avrg_compl);
227    rhs_cen->z_U_NonConst()->Set(avrg_compl);
228    rhs_cen->v_L_NonConst()->Set(avrg_compl);
229    rhs_cen->v_U_NonConst()->Set(avrg_compl);
230
231    // Get space for the centering step
232    SmartPtr<IteratesVector> step_cen = IpData().curr()->MakeNewIteratesVector(true);
233
234    // Now solve the primal-dual system to get the step
235    allow_inexact = true;
236    pd_solver_->Solve(1.0, 0.0,
237                      *rhs_cen,
238                      *step_cen,
239                      allow_inexact
240                     );
241
242    DBG_PRINT_VECTOR(2, "step_cen", *step_cen);
243
244    // Start the timing for the quality function search here
245    IpData().TimingStats().QualityFunctionSearch().Start();
246
247    // Some initializations
248    if (!initialized_) {
249      n_dual_ = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
250      n_pri_ = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
251      n_comp_ = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
252                IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
253      initialized_ = true;
254    }
255
256    count_qf_evals_ = 0;
257
258    // Compute some quantities used for the quality function evaluations
259    // (This way we try to avoid retrieving numbers from cache...
260
261    curr_slack_x_L_ = IpCq().curr_slack_x_L();
262    curr_slack_x_U_ = IpCq().curr_slack_x_U();
263    curr_slack_s_L_ = IpCq().curr_slack_s_L();
264    curr_slack_s_U_ = IpCq().curr_slack_s_U();
265
266    curr_z_L_ = IpData().curr()->z_L();
267    curr_z_U_ = IpData().curr()->z_U();
268    curr_v_L_ = IpData().curr()->v_L();
269    curr_v_U_ = IpData().curr()->v_U();
270
271    IpData().TimingStats().Task5().Start();
272    switch (quality_function_norm_) {
273      case NM_NORM_1:
274      curr_grad_lag_x_asum_ = IpCq().curr_grad_lag_x()->Asum();
275      curr_grad_lag_s_asum_ = IpCq().curr_grad_lag_s()->Asum();
276      curr_c_asum_ = IpCq().curr_c()->Asum();
277      curr_d_minus_s_asum_ = IpCq().curr_d_minus_s()->Asum();
278      break;
279      case NM_NORM_2_SQUARED:
280      case NM_NORM_2:
281      curr_grad_lag_x_nrm2_ = IpCq().curr_grad_lag_x()->Nrm2();
282      curr_grad_lag_s_nrm2_ = IpCq().curr_grad_lag_s()->Nrm2();
283      curr_c_nrm2_ = IpCq().curr_c()->Nrm2();
284      curr_d_minus_s_nrm2_ = IpCq().curr_d_minus_s()->Nrm2();
285      break;
286      case NM_NORM_MAX:
287      curr_grad_lag_x_amax_ = IpCq().curr_grad_lag_x()->Amax();
288      curr_grad_lag_s_amax_ = IpCq().curr_grad_lag_s()->Amax();
289      curr_c_amax_ = IpCq().curr_c()->Amax();
290      curr_d_minus_s_amax_ = IpCq().curr_d_minus_s()->Amax();
291      break;
292      default:
293      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
294    }
295    IpData().TimingStats().Task5().End();
296
297    // Some initializations
298    if (!initialized_) {
299      n_dual_ = IpData().curr()->x()->Dim() + IpData().curr()->s()->Dim();
300      n_pri_ = IpData().curr()->y_c()->Dim() + IpData().curr()->y_d()->Dim();
301      n_comp_ = IpData().curr()->z_L()->Dim() + IpData().curr()->z_U()->Dim() +
302                IpData().curr()->v_L()->Dim() + IpData().curr()->v_U()->Dim();
303      initialized_ = true;
304    }
305
306    count_qf_evals_ = 0;
307
308    // Compute some quantities used for the quality function evaluations
309    // (This way we try to avoid retrieving numbers from cache...
310
311    curr_slack_x_L_ = IpCq().curr_slack_x_L();
312    curr_slack_x_U_ = IpCq().curr_slack_x_U();
313    curr_slack_s_L_ = IpCq().curr_slack_s_L();
314    curr_slack_s_U_ = IpCq().curr_slack_s_U();
315
316    curr_z_L_ = IpData().curr()->z_L();
317    curr_z_U_ = IpData().curr()->z_U();
318    curr_v_L_ = IpData().curr()->v_L();
319    curr_v_U_ = IpData().curr()->v_U();
320
321    IpData().TimingStats().Task5().Start();
322    switch (quality_function_norm_) {
323      case NM_NORM_1:
324      curr_grad_lag_x_asum_ = IpCq().curr_grad_lag_x()->Asum();
325      curr_grad_lag_s_asum_ = IpCq().curr_grad_lag_s()->Asum();
326      curr_c_asum_ = IpCq().curr_c()->Asum();
327      curr_d_minus_s_asum_ = IpCq().curr_d_minus_s()->Asum();
328      break;
329      case NM_NORM_2_SQUARED:
330      case NM_NORM_2:
331      curr_grad_lag_x_nrm2_ = IpCq().curr_grad_lag_x()->Nrm2();
332      curr_grad_lag_s_nrm2_ = IpCq().curr_grad_lag_s()->Nrm2();
333      curr_c_nrm2_ = IpCq().curr_c()->Nrm2();
334      curr_d_minus_s_nrm2_ = IpCq().curr_d_minus_s()->Nrm2();
335      break;
336      case NM_NORM_MAX:
337      curr_grad_lag_x_amax_ = IpCq().curr_grad_lag_x()->Amax();
338      curr_grad_lag_s_amax_ = IpCq().curr_grad_lag_s()->Amax();
339      curr_c_amax_ = IpCq().curr_c()->Amax();
340      curr_d_minus_s_amax_ = IpCq().curr_d_minus_s()->Amax();
341      break;
342      default:
343      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
344    }
345    IpData().TimingStats().Task5().End();
346
347    // We now compute the step for the slack variables.  This safes
348    // time, because we then don't have to do this any more for each
349    // evaluation of the quality function
350    SmartPtr<Vector> step_aff_x_L = step_aff->z_L()->MakeNew();
351    SmartPtr<Vector> step_aff_x_U = step_aff->z_U()->MakeNew();
352    SmartPtr<Vector> step_aff_s_L = step_aff->v_L()->MakeNew();
353    SmartPtr<Vector> step_aff_s_U = step_aff->v_U()->MakeNew();
354    IpNLP().Px_L()->TransMultVector(1., *step_aff->x(), 0., *step_aff_x_L);
355    IpNLP().Px_U()->TransMultVector(-1., *step_aff->x(), 0., *step_aff_x_U);
356    IpNLP().Pd_L()->TransMultVector(1., *step_aff->s(), 0., *step_aff_s_L);
357    IpNLP().Pd_U()->TransMultVector(-1., *step_aff->s(), 0., *step_aff_s_U);
358    SmartPtr<Vector> step_cen_x_L = step_cen->z_L()->MakeNew();
359    SmartPtr<Vector> step_cen_x_U = step_cen->z_U()->MakeNew();
360    SmartPtr<Vector> step_cen_s_L = step_cen->v_L()->MakeNew();
361    SmartPtr<Vector> step_cen_s_U = step_cen->v_U()->MakeNew();
362    IpNLP().Px_L()->TransMultVector(1., *step_cen->x(), 0., *step_cen_x_L);
363    IpNLP().Px_U()->TransMultVector(-1., *step_cen->x(), 0., *step_cen_x_U);
364    IpNLP().Pd_L()->TransMultVector(1., *step_cen->s(), 0., *step_cen_s_L);
365    IpNLP().Pd_U()->TransMultVector(-1., *step_cen->s(), 0., *step_cen_s_U);
366
367    Number sigma;
368
369    // First we determine whether we want to search for a value of
370    // sigma larger or smaller than 1.  For this, we estimate the
371    // slope of the quality function at sigma=1.
372    Number qf_1 = CalculateQualityFunction(1.,
373                                           *step_aff_x_L,
374                                           *step_aff_x_U,
375                                           *step_aff_s_L,
376                                           *step_aff_s_U,
377                                           *step_aff->y_c(),
378                                           *step_aff->y_d(),
379                                           *step_aff->z_L(),
380                                           *step_aff->z_U(),
381                                           *step_aff->v_L(),
382                                           *step_aff->v_U(),
383                                           *step_cen_x_L,
384                                           *step_cen_x_U,
385                                           *step_cen_s_L,
386                                           *step_cen_s_U,
387                                           *step_cen->y_c(),
388                                           *step_cen->y_d(),
389                                           *step_cen->z_L(),
390                                           *step_cen->z_U(),
391                                           *step_cen->v_L(),
392                                           *step_cen->v_U());
393
394    Number sigma_1minus = 1.-Max(1e-4, quality_function_section_sigma_tol_);
395    Number qf_1minus = CalculateQualityFunction(sigma_1minus,
396                       *step_aff_x_L,
397                       *step_aff_x_U,
398                       *step_aff_s_L,
399                       *step_aff_s_U,
400                       *step_aff->y_c(),
401                       *step_aff->y_d(),
402                       *step_aff->z_L(),
403                       *step_aff->z_U(),
404                       *step_aff->v_L(),
405                       *step_aff->v_U(),
406                       *step_cen_x_L,
407                       *step_cen_x_U,
408                       *step_cen_s_L,
409                       *step_cen_s_U,
410                       *step_cen->y_c(),
411                       *step_cen->y_d(),
412                       *step_cen->z_L(),
413                       *step_cen->z_U(),
414                       *step_cen->v_L(),
415                       *step_cen->v_U());
416
417    if (qf_1minus > qf_1) {
418      // It seems that the quality function decreases for values
419      // larger than sigma, so perform golden section search for sigma
420      // > 1.
421      Number sigma_up = Min(sigma_max_, mu_max/avrg_compl);
422      Number sigma_lo = 1.;
423      if (sigma_lo >= sigma_up) {
424        sigma = sigma_up;
425      }
426      else {
427        // ToDo maybe we should use different tolerances for sigma>1
428        sigma = PerformGoldenSection(sigma_up, -100., sigma_lo, qf_1,
429                                     quality_function_section_sigma_tol_,
430                                     quality_function_section_qf_tol_,
431                                     *step_aff_x_L,
432                                     *step_aff_x_U,
433                                     *step_aff_s_L,
434                                     *step_aff_s_U,
435                                     *step_aff->y_c(),
436                                     *step_aff->y_d(),
437                                     *step_aff->z_L(),
438                                     *step_aff->z_U(),
439                                     *step_aff->v_L(),
440                                     *step_aff->v_U(),
441                                     *step_cen_x_L,
442                                     *step_cen_x_U,
443                                     *step_cen_s_L,
444                                     *step_cen_s_U,
445                                     *step_cen->y_c(),
446                                     *step_cen->y_d(),
447                                     *step_cen->z_L(),
448                                     *step_cen->z_U(),
449                                     *step_cen->v_L(),
450                                     *step_cen->v_U());
451      }
452    }
453    else {
454      // Search for sigma less than 1
455
456      Number sigma_lo = Max(sigma_min_, mu_min/avrg_compl);
457      Number sigma_up = Min(Max(sigma_lo, sigma_1minus), mu_max/avrg_compl);
458      if (sigma_lo >= sigma_up) {
459        // Skip the search, we are already at the minimum
460        sigma = sigma_lo;
461      }
462      else {
463        sigma = PerformGoldenSection(sigma_up, qf_1minus, sigma_lo, -100.,
464                                     quality_function_section_sigma_tol_,
465                                     quality_function_section_qf_tol_,
466                                     *step_aff_x_L,
467                                     *step_aff_x_U,
468                                     *step_aff_s_L,
469                                     *step_aff_s_U,
470                                     *step_aff->y_c(),
471                                     *step_aff->y_d(),
472                                     *step_aff->z_L(),
473                                     *step_aff->z_U(),
474                                     *step_aff->v_L(),
475                                     *step_aff->v_U(),
476                                     *step_cen_x_L,
477                                     *step_cen_x_U,
478                                     *step_cen_s_L,
479                                     *step_cen_s_U,
480                                     *step_cen->y_c(),
481                                     *step_cen->y_d(),
482                                     *step_cen->z_L(),
483                                     *step_cen->z_U(),
484                                     *step_cen->v_L(),
485                                     *step_cen->v_U());
486      }
487    }
488
489    //#define tracequalityfunction
490#ifdef tracequalityfunction
491    char fname[100];
492    sprintf(fname, "qf_values_%d.dat", IpData().iter_count());
493    FILE* fid = fopen(fname, "w");
494
495    Number sigma_1 = sigma_max_;
496    Number sigma_2 = 1e-9/avrg_compl;
497    Number sigma_trace = sigma_1;
498    while(sigma_trace > sigma_2) {
499      Number qf = CalculateQualityFunction(sigma_trace,
500                                           *step_aff_x_L,
501                                           *step_aff_x_U,
502                                           *step_aff_s_L,
503                                           *step_aff_s_U,
504                                           *step_aff->y_c(),
505                                           *step_aff->y_d(),
506                                           *step_aff->z_L(),
507                                           *step_aff->z_U(),
508                                           *step_aff->v_L(),
509                                           *step_aff->v_U(),
510                                           *step_cen_x_L,
511                                           *step_cen_x_U,
512                                           *step_cen_s_L,
513                                           *step_cen_s_U,
514                                           *step_cen->y_c(),
515                                           *step_cen->y_d(),
516                                           *step_cen->z_L(),
517                                           *step_cen->z_U(),
518                                           *step_cen->v_L(),
519                                           *step_cen->v_U());
520      fprintf(fid, "%9.2e %25.16e\n", sigma_trace, qf);
521      sigma_trace /= 1.1;
522    }
523    fclose(fid);
524#endif
525
526    // End timing of quality function search
527    IpData().TimingStats().QualityFunctionSearch().End();
528
529    Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE,
530                   "Sigma = %e\n", sigma);
531    Number mu = sigma*avrg_compl;
532
533    // Store the affine search direction (in case it is needed in the
534    // line search for a corrector step)
535    IpData().set_delta_aff(step_aff);
536    IpData().SetHaveAffineDeltas(true);
537
538    // Now construct the overall search direction here
539    SmartPtr<IteratesVector> step = IpData().curr()->MakeNewIteratesVector(true);
540    step->AddTwoVectors(sigma, *step_cen, 1.0, *IpData().delta_aff(), 0.0);
541
542    DBG_PRINT_VECTOR(2, "step", *step);
543    IpData().set_delta(step);
544    IpData().SetHaveDeltas(true);
545
546    ///////////////////////////////////////////////////////////////////////////
547    // Release memory for temporary vectors used in CalculateQualityFunction //
548    ///////////////////////////////////////////////////////////////////////////
549
550    tmp_step_x_L_ = NULL;
551    tmp_step_x_U_ = NULL;
552    tmp_step_s_L_ = NULL;
553    tmp_step_s_U_ = NULL;
554    tmp_step_z_L_ = NULL;
555    tmp_step_z_U_ = NULL;
556    tmp_step_v_L_ = NULL;
557    tmp_step_v_U_ = NULL;
558
559    tmp_slack_x_L_ = NULL;
560    tmp_slack_x_U_ = NULL;
561    tmp_slack_s_L_ = NULL;
562    tmp_slack_s_U_ = NULL;
563    tmp_z_L_ = NULL;
564    tmp_z_U_ = NULL;
565    tmp_v_L_ = NULL;
566    tmp_v_U_ = NULL;
567
568    curr_slack_x_L_ = NULL;
569    curr_slack_x_U_ = NULL;
570    curr_slack_s_L_ = NULL;
571    curr_slack_s_U_ = NULL;
572
573    // DELETEME
574    char ssigma[40];
575    sprintf(ssigma, " sigma=%8.2e", sigma);
576    IpData().Append_info_string(ssigma);
577    sprintf(ssigma, " qf=%d", count_qf_evals_);
578    IpData().Append_info_string(ssigma);
579    /*
580    sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
581    IpData().Append_info_string(ssigma);
582    if (sigma>1.) {
583      IpData().Append_info_string("LARGESIGMA");
584    }
585    */
586
587    return mu;
588  }
589
590  Number QualityFunctionMuOracle::CalculateQualityFunction
591  (Number sigma,
592   const Vector& step_aff_x_L,
593   const Vector& step_aff_x_U,
594   const Vector& step_aff_s_L,
595   const Vector& step_aff_s_U,
596   const Vector& step_aff_y_c,
597   const Vector& step_aff_y_d,
598   const Vector& step_aff_z_L,
599   const Vector& step_aff_z_U,
600   const Vector& step_aff_v_L,
601   const Vector& step_aff_v_U,
602   const Vector& step_cen_x_L,
603   const Vector& step_cen_x_U,
604   const Vector& step_cen_s_L,
605   const Vector& step_cen_s_U,
606   const Vector& step_cen_y_c,
607   const Vector& step_cen_y_d,
608   const Vector& step_cen_z_L,
609   const Vector& step_cen_z_U,
610   const Vector& step_cen_v_L,
611   const Vector& step_cen_v_U
612  )
613  {
614    DBG_START_METH("QualityFunctionMuOracle::CalculateQualityFunction",
615                   dbg_verbosity);
616    count_qf_evals_++;
617
618    IpData().TimingStats().Task1().Start();
619    tmp_step_x_L_->AddTwoVectors(1., step_aff_x_L, sigma, step_cen_x_L, 0.);
620    tmp_step_x_U_->AddTwoVectors(1., step_aff_x_U, sigma, step_cen_x_U, 0.);
621    tmp_step_s_L_->AddTwoVectors(1., step_aff_s_L, sigma, step_cen_s_L, 0.);
622    tmp_step_s_U_->AddTwoVectors(1., step_aff_s_U, sigma, step_cen_s_U, 0.);
623    tmp_step_z_L_->AddTwoVectors(1., step_aff_z_L, sigma, step_cen_z_L, 0.);
624    tmp_step_z_U_->AddTwoVectors(1., step_aff_z_U, sigma, step_cen_z_U, 0.);
625    tmp_step_v_L_->AddTwoVectors(1., step_aff_v_L, sigma, step_cen_v_L, 0.);
626    tmp_step_v_U_->AddTwoVectors(1., step_aff_v_U, sigma, step_cen_v_U, 0.);
627    IpData().TimingStats().Task1().End();
628
629    // Compute the fraction-to-the-boundary step sizes
630    IpData().TimingStats().Task2().Start();
631    Number tau = IpData().curr_tau();
632    Number alpha_primal = IpCq().uncached_slack_frac_to_the_bound(tau,
633                          *tmp_step_x_L_,
634                          *tmp_step_x_U_,
635                          *tmp_step_s_L_,
636                          *tmp_step_s_U_);
637
638    Number alpha_dual = IpCq().uncached_dual_frac_to_the_bound(tau,
639                        *tmp_step_z_L_,
640                        *tmp_step_z_U_,
641                        *tmp_step_v_L_,
642                        *tmp_step_v_U_);
643    IpData().TimingStats().Task2().End();
644
645    Number xi = 0.; // centrality measure
646
647    IpData().TimingStats().Task1().Start();
648    tmp_slack_x_L_->AddTwoVectors(1., *curr_slack_x_L_,
649                                  alpha_primal, *tmp_step_x_L_, 0.);
650    tmp_slack_x_U_->AddTwoVectors(1., *curr_slack_x_U_,
651                                  alpha_primal, *tmp_step_x_U_, 0.);
652    tmp_slack_s_L_->AddTwoVectors(1., *curr_slack_s_L_,
653                                  alpha_primal, *tmp_step_s_L_, 0.);
654    tmp_slack_s_U_->AddTwoVectors(1., *curr_slack_s_U_,
655                                  alpha_primal, *tmp_step_s_U_, 0.);
656
657    tmp_z_L_->AddTwoVectors(1., *curr_z_L_,
658                            alpha_dual, *tmp_step_z_L_, 0.);
659    tmp_z_U_->AddTwoVectors(1., *curr_z_U_,
660                            alpha_dual, *tmp_step_z_U_, 0.);
661    tmp_v_L_->AddTwoVectors(1., *curr_v_L_,
662                            alpha_dual, *tmp_step_v_L_, 0.);
663    tmp_v_U_->AddTwoVectors(1., *curr_v_U_,
664                            alpha_dual, *tmp_step_v_U_, 0.);
665    IpData().TimingStats().Task1().End();
666
667    IpData().TimingStats().Task3().Start();
668    tmp_slack_x_L_->ElementWiseMultiply(*tmp_z_L_);
669    tmp_slack_x_U_->ElementWiseMultiply(*tmp_z_U_);
670    tmp_slack_s_L_->ElementWiseMultiply(*tmp_v_L_);
671    tmp_slack_s_U_->ElementWiseMultiply(*tmp_v_U_);
672    IpData().TimingStats().Task3().End();
673
674    DBG_PRINT_VECTOR(2, "compl_x_L", *tmp_slack_x_L_);
675    DBG_PRINT_VECTOR(2, "compl_x_U", *tmp_slack_x_U_);
676    DBG_PRINT_VECTOR(2, "compl_s_L", *tmp_slack_s_L_);
677    DBG_PRINT_VECTOR(2, "compl_s_U", *tmp_slack_s_U_);
678
679    Number dual_inf=-1.;
680    Number primal_inf=-1.;
681    Number compl_inf=-1.;
682
683    IpData().TimingStats().Task5().Start();
684    switch (quality_function_norm_) {
685      case NM_NORM_1:
686      dual_inf = (1.-alpha_dual)*(curr_grad_lag_x_asum_ +
687                                  curr_grad_lag_s_asum_);
688
689      primal_inf = (1.-alpha_primal)*(curr_c_asum_ +
690                                      curr_d_minus_s_asum_);
691
692      compl_inf = tmp_slack_x_L_->Asum() + tmp_slack_x_U_->Asum() +
693                  tmp_slack_s_L_->Asum() + tmp_slack_s_U_->Asum();
694
695      dual_inf /= n_dual_;
696      if (n_pri_>0) {
697        primal_inf /= n_pri_;
698      }
699      DBG_ASSERT(n_comp_>0);
700      compl_inf /= n_comp_;
701      break;
702      case NM_NORM_2_SQUARED:
703      dual_inf =
704        pow(1.-alpha_dual, 2)*(pow(curr_grad_lag_x_nrm2_, 2) +
705                               pow(curr_grad_lag_s_nrm2_, 2));
706      primal_inf =
707        pow(1.-alpha_primal, 2)*(pow(curr_c_nrm2_, 2) +
708                                 pow(curr_d_minus_s_nrm2_, 2));
709      compl_inf =
710        pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
711        pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2);
712
713      dual_inf /= n_dual_;
714      if (n_pri_>0) {
715        primal_inf /= n_pri_;
716      }
717      DBG_ASSERT(n_comp_>0);
718      compl_inf /= n_comp_;
719      break;
720      case NM_NORM_MAX:
721      dual_inf =
722        (1.-alpha_dual)*Max(curr_grad_lag_x_amax_,
723                            curr_grad_lag_s_amax_);
724      primal_inf =
725        (1.-alpha_primal)*Max(curr_c_amax_,
726                              curr_d_minus_s_amax_);
727      compl_inf =
728        Max(tmp_slack_x_L_->Amax(), tmp_slack_x_U_->Amax(),
729            tmp_slack_s_L_->Amax(), tmp_slack_s_U_->Amax());
730      break;
731      case NM_NORM_2:
732      dual_inf =
733        (1.-alpha_dual)*sqrt(pow(curr_grad_lag_x_nrm2_, 2) +
734                             pow(curr_grad_lag_s_nrm2_, 2));
735      primal_inf =
736        (1.-alpha_primal)*sqrt(pow(curr_c_nrm2_, 2) +
737                               pow(curr_d_minus_s_nrm2_, 2));
738      compl_inf =
739        sqrt(pow(tmp_slack_x_L_->Nrm2(), 2) + pow(tmp_slack_x_U_->Nrm2(), 2) +
740             pow(tmp_slack_s_L_->Nrm2(), 2) + pow(tmp_slack_s_U_->Nrm2(), 2));
741
742      dual_inf /= sqrt((Number)n_dual_);
743      if (n_pri_>0) {
744        primal_inf /= sqrt((Number)n_pri_);
745      }
746      DBG_ASSERT(n_comp_>0);
747      compl_inf /= sqrt((Number)n_comp_);
748      break;
749      default:
750      DBG_ASSERT(false && "Unknown value for quality_function_norm_");
751    }
752    IpData().TimingStats().Task5().End();
753
754    Number quality_function = dual_inf + primal_inf + compl_inf;
755
756    if (quality_function_centrality_!=CEN_NONE) {
757      IpData().TimingStats().Task4().Start();
758      xi = IpCq().CalcCentralityMeasure(*tmp_slack_x_L_, *tmp_slack_x_U_,
759                                        *tmp_slack_s_L_, *tmp_slack_s_U_);
760      IpData().TimingStats().Task4().End();
761    }
762    switch (quality_function_centrality_) {
763      case CEN_NONE:
764      //Nothing
765      break;
766      case CEN_LOG:
767      quality_function -= compl_inf*log(xi);
768      break;
769      case CEN_RECIPROCAL:
770      quality_function += compl_inf/xi;
771      case CEN_CUBED_RECIPROCAL:
772      quality_function += compl_inf/pow(xi,3);
773      break;
774      default:
775      DBG_ASSERT(false && "Unknown value for quality_function_centrality_");
776    }
777
778    switch (quality_function_balancing_term_) {
779      case BT_NONE:
780      //Nothing
781      break;
782      case BT_CUBIC:
783      quality_function += pow(Max(0., Max(dual_inf,primal_inf)-compl_inf),3);
784      break;
785      default:
786      DBG_ASSERT(false && "Unknown value for quality_function_balancing term_");
787    }
788
789    Jnlst().Printf(J_MOREDETAILED, J_BARRIER_UPDATE,
790                   "sigma = %8.2e d_inf = %18.12e p_inf = %18.12e cmpl = %18.12e q = %18.12e a_pri = %8.2e a_dual = %8.2e xi = %8.2e\n", sigma, dual_inf, primal_inf, compl_inf, quality_function, alpha_primal, alpha_dual, xi);
791
792
793    return quality_function;
794    //return compl_inf;
795  }
796
797  Number
798  QualityFunctionMuOracle::PerformGoldenSection
799  (Number sigma_up_in,
800   Number q_up,
801   Number sigma_lo_in,
802   Number q_lo,
803   Number sigma_tol,
804   Number qf_tol,
805   const Vector& step_aff_x_L,
806   const Vector& step_aff_x_U,
807   const Vector& step_aff_s_L,
808   const Vector& step_aff_s_U,
809   const Vector& step_aff_y_c,
810   const Vector& step_aff_y_d,
811   const Vector& step_aff_z_L,
812   const Vector& step_aff_z_U,
813   const Vector& step_aff_v_L,
814   const Vector& step_aff_v_U,
815   const Vector& step_cen_x_L,
816   const Vector& step_cen_x_U,
817   const Vector& step_cen_s_L,
818   const Vector& step_cen_s_U,
819   const Vector& step_cen_y_c,
820   const Vector& step_cen_y_d,
821   const Vector& step_cen_z_L,
822   const Vector& step_cen_z_U,
823   const Vector& step_cen_v_L,
824   const Vector& step_cen_v_U
825  )
826  {
827    Number sigma_up = ScaleSigma(sigma_up_in);
828    Number sigma_lo = ScaleSigma(sigma_lo_in);
829
830    Number sigma;
831    Number gfac = (3.-sqrt(5.))/2.;
832    Number sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
833    Number sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
834
835    Number qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
836                                            step_aff_x_L,
837                                            step_aff_x_U,
838                                            step_aff_s_L,
839                                            step_aff_s_U,
840                                            step_aff_y_c,
841                                            step_aff_y_d,
842                                            step_aff_z_L,
843                                            step_aff_z_U,
844                                            step_aff_v_L,
845                                            step_aff_v_U,
846                                            step_cen_x_L,
847                                            step_cen_x_U,
848                                            step_cen_s_L,
849                                            step_cen_s_U,
850                                            step_cen_y_c,
851                                            step_cen_y_d,
852                                            step_cen_z_L,
853                                            step_cen_z_U,
854                                            step_cen_v_L,
855                                            step_cen_v_U);
856    Number qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
857                                            step_aff_x_L,
858                                            step_aff_x_U,
859                                            step_aff_s_L,
860                                            step_aff_s_U,
861                                            step_aff_y_c,
862                                            step_aff_y_d,
863                                            step_aff_z_L,
864                                            step_aff_z_U,
865                                            step_aff_v_L,
866                                            step_aff_v_U,
867                                            step_cen_x_L,
868                                            step_cen_x_U,
869                                            step_cen_s_L,
870                                            step_cen_s_U,
871                                            step_cen_y_c,
872                                            step_cen_y_d,
873                                            step_cen_z_L,
874                                            step_cen_z_U,
875                                            step_cen_v_L,
876                                            step_cen_v_U);
877
878    Index nsections = 0;
879    while ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
880           //while ((sigma_up-sigma_lo)>=sigma_tol &&  // Note we are using the non-relative criterion here for sigma
881           (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))>=qf_tol &&
882           nsections<quality_function_max_section_steps_) {
883      nsections++;
884      //printf("sigma_lo=%e sigma_mid1=%e sigma_mid2=%e sigma_up=%e\n",sigma_lo, sigma_mid1, sigma_mid2, sigma_up);
885      if (qmid1 > qmid2) {
886        sigma_lo = sigma_mid1;
887        q_lo = qmid1;
888        sigma_mid1 = sigma_mid2;
889        qmid1 = qmid2;
890        sigma_mid2 = sigma_lo + (1.-gfac)*(sigma_up-sigma_lo);
891        qmid2 = CalculateQualityFunction(UnscaleSigma(sigma_mid2),
892                                         step_aff_x_L,
893                                         step_aff_x_U,
894                                         step_aff_s_L,
895                                         step_aff_s_U,
896                                         step_aff_y_c,
897                                         step_aff_y_d,
898                                         step_aff_z_L,
899                                         step_aff_z_U,
900                                         step_aff_v_L,
901                                         step_aff_v_U,
902                                         step_cen_x_L,
903                                         step_cen_x_U,
904                                         step_cen_s_L,
905                                         step_cen_s_U,
906                                         step_cen_y_c,
907                                         step_cen_y_d,
908                                         step_cen_z_L,
909                                         step_cen_z_U,
910                                         step_cen_v_L,
911                                         step_cen_v_U);
912      }
913      else {
914        sigma_up = sigma_mid2;
915        q_up = qmid2;
916        sigma_mid2 = sigma_mid1;
917        qmid2 = qmid1;
918        sigma_mid1 = sigma_lo + gfac*(sigma_up-sigma_lo);
919        qmid1 = CalculateQualityFunction(UnscaleSigma(sigma_mid1),
920                                         step_aff_x_L,
921                                         step_aff_x_U,
922                                         step_aff_s_L,
923                                         step_aff_s_U,
924                                         step_aff_y_c,
925                                         step_aff_y_d,
926                                         step_aff_z_L,
927                                         step_aff_z_U,
928                                         step_aff_v_L,
929                                         step_aff_v_U,
930                                         step_cen_x_L,
931                                         step_cen_x_U,
932                                         step_cen_s_L,
933                                         step_cen_s_U,
934                                         step_cen_y_c,
935                                         step_cen_y_d,
936                                         step_cen_z_L,
937                                         step_cen_z_U,
938                                         step_cen_v_L,
939                                         step_cen_v_U);
940      }
941    }
942
943    if ((sigma_up-sigma_lo)>=sigma_tol*sigma_up &&
944        (1.-Min(q_lo, q_up, qmid1, qmid2)/Max(q_lo, q_up, qmid1, qmid2))<qf_tol) {
945      // The qf tolerance make it stop
946      IpData().Append_info_string("qf_tol ");
947      Number qf_min = Min(q_lo, q_up, qmid1, qmid2);
948      DBG_ASSERT(qf_min>-100.);
949      if (qf_min == q_lo) {
950        sigma = sigma_lo;
951      }
952      else if (qf_min == qmid1) {
953        sigma = sigma_mid1;
954      }
955      else if (qf_min == qmid2) {
956        sigma = sigma_mid2;
957      }
958      else {
959        sigma = sigma_up;
960      }
961    }
962    else {
963      Number q;
964      if (qmid1 < qmid2) {
965        sigma = sigma_mid1;
966        q = qmid1;
967      }
968      else {
969        sigma = sigma_mid2;
970        q = qmid2;
971      }
972      if (sigma_up == ScaleSigma(sigma_up_in)) {
973        Number qtmp;
974        if (q_up<0.) {
975          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_up),
976                                          step_aff_x_L,
977                                          step_aff_x_U,
978                                          step_aff_s_L,
979                                          step_aff_s_U,
980                                          step_aff_y_c,
981                                          step_aff_y_d,
982                                          step_aff_z_L,
983                                          step_aff_z_U,
984                                          step_aff_v_L,
985                                          step_aff_v_U,
986                                          step_cen_x_L,
987                                          step_cen_x_U,
988                                          step_cen_s_L,
989                                          step_cen_s_U,
990                                          step_cen_y_c,
991                                          step_cen_y_d,
992                                          step_cen_z_L,
993                                          step_cen_z_U,
994                                          step_cen_v_L,
995                                          step_cen_v_U);
996        }
997        else {
998          qtmp = q_up;
999        }
1000        if (qtmp < q) {
1001          sigma = sigma_up;
1002          q = qtmp;
1003        }
1004      }
1005      else if (sigma_lo == ScaleSigma(sigma_lo_in)) {
1006        Number qtmp;
1007        if (q_lo<0.) {
1008          qtmp = CalculateQualityFunction(UnscaleSigma(sigma_lo),
1009                                          step_aff_x_L,
1010                                          step_aff_x_U,
1011                                          step_aff_s_L,
1012                                          step_aff_s_U,
1013                                          step_aff_y_c,
1014                                          step_aff_y_d,
1015                                          step_aff_z_L,
1016                                          step_aff_z_U,
1017                                          step_aff_v_L,
1018                                          step_aff_v_U,
1019                                          step_cen_x_L,
1020                                          step_cen_x_U,
1021                                          step_cen_s_L,
1022                                          step_cen_s_U,
1023                                          step_cen_y_c,
1024                                          step_cen_y_d,
1025                                          step_cen_z_L,
1026                                          step_cen_z_U,
1027                                          step_cen_v_L,
1028                                          step_cen_v_U);
1029        }
1030        else {
1031          qtmp = q_lo;
1032        }
1033        if (qtmp < q) {
1034          sigma = sigma_lo;
1035          q = qtmp;
1036        }
1037      }
1038    }
1039
1040    return UnscaleSigma(sigma);
1041  }
1042
1043  /*
1044  Number QualityFunctionMuOracle::ScaleSigma(Number sigma) {return log(sigma);}
1045  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma) {return exp(scaled_sigma);}
1046  */
1047  Number QualityFunctionMuOracle::ScaleSigma(Number sigma)
1048  {
1049    return sigma;
1050  }
1051  Number QualityFunctionMuOracle::UnscaleSigma(Number scaled_sigma)
1052  {
1053    return scaled_sigma;
1054  }
1055
1056  /* AW: Tried search in the log space, but that was even worse than
1057     search in unscaled space */
1058  /*
1059  Number
1060  QualityFunctionMuOracle::PerformGoldenSectionLog
1061  (Number sigma_up,
1062   Number sigma_lo,
1063   Number tol,
1064   const Vector& step_aff_x_L,
1065   const Vector& step_aff_x_U,
1066   const Vector& step_aff_s_L,
1067   const Vector& step_aff_s_U,
1068   const Vector& step_aff_y_c,
1069   const Vector& step_aff_y_d,
1070   const Vector& step_aff_z_L,
1071   const Vector& step_aff_z_U,
1072   const Vector& step_aff_v_L,
1073   const Vector& step_aff_v_U,
1074   const Vector& step_cen_x_L,
1075   const Vector& step_cen_x_U,
1076   const Vector& step_cen_s_L,
1077   const Vector& step_cen_s_U,
1078   const Vector& step_cen_y_c,
1079   const Vector& step_cen_y_d,
1080   const Vector& step_cen_z_L,
1081   const Vector& step_cen_z_U,
1082   const Vector& step_cen_v_L,
1083   const Vector& step_cen_v_U
1084  )
1085  {
1086    Number log_sigma;
1087    Number log_sigma_up = log(sigma_up);
1088    Number log_sigma_lo = log(sigma_lo);
1089
1090    Number log_sigma_up_in = log_sigma_up;
1091    Number log_sigma_lo_in = log_sigma_lo;
1092    Number gfac = (3.-sqrt(5.))/2.;
1093    Number log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo);
1094    Number log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo);
1095
1096    Number qmid1 = CalculateQualityFunction(exp(log_sigma_mid1),
1097                                            step_aff_x_L,
1098                                            step_aff_x_U,
1099                                            step_aff_s_L,
1100                                            step_aff_s_U,
1101                                            step_aff_y_c,
1102                                            step_aff_y_d,
1103                                            step_aff_z_L,
1104                                            step_aff_z_U,
1105                                            step_aff_v_L,
1106                                            step_aff_v_U,
1107                                            step_cen_x_L,
1108                                            step_cen_x_U,
1109                                            step_cen_s_L,
1110                                            step_cen_s_U,
1111                                            step_cen_y_c,
1112                                            step_cen_y_d,
1113                                            step_cen_z_L,
1114                                            step_cen_z_U,
1115                                            step_cen_v_L,
1116                                            step_cen_v_U);
1117    Number qmid2 = CalculateQualityFunction(exp(log_sigma_mid2),
1118                                            step_aff_x_L,
1119                                            step_aff_x_U,
1120                                            step_aff_s_L,
1121                                            step_aff_s_U,
1122                                            step_aff_y_c,
1123                                            step_aff_y_d,
1124                                            step_aff_z_L,
1125                                            step_aff_z_U,
1126                                            step_aff_v_L,
1127                                            step_aff_v_U,
1128                                            step_cen_x_L,
1129                                            step_cen_x_U,
1130                                            step_cen_s_L,
1131                                            step_cen_s_U,
1132                                            step_cen_y_c,
1133                                            step_cen_y_d,
1134                                            step_cen_z_L,
1135                                            step_cen_z_U,
1136                                            step_cen_v_L,
1137                                            step_cen_v_U);
1138
1139    Index nsections = 0;
1140    while ((log_sigma_up-log_sigma_lo)>=tol*log_sigma_up && nsections<quality_function_max_section_steps_) {
1141      nsections++;
1142      if (qmid1 > qmid2) {
1143        log_sigma_lo = log_sigma_mid1;
1144        log_sigma_mid1 = log_sigma_mid2;
1145        qmid1 = qmid2;
1146        log_sigma_mid2 = log_sigma_lo + (1.-gfac)*(log_sigma_up-log_sigma_lo);
1147        qmid2 = CalculateQualityFunction(exp(log_sigma_mid2),
1148                                         step_aff_x_L,
1149                                         step_aff_x_U,
1150                                         step_aff_s_L,
1151                                         step_aff_s_U,
1152                                         step_aff_y_c,
1153                                         step_aff_y_d,
1154                                         step_aff_z_L,
1155                                         step_aff_z_U,
1156                                         step_aff_v_L,
1157                                         step_aff_v_U,
1158                                         step_cen_x_L,
1159                                         step_cen_x_U,
1160                                         step_cen_s_L,
1161                                         step_cen_s_U,
1162                                         step_cen_y_c,
1163                                         step_cen_y_d,
1164                                         step_cen_z_L,
1165                                         step_cen_z_U,
1166                                         step_cen_v_L,
1167                                         step_cen_v_U);
1168      }
1169      else {
1170        log_sigma_up = log_sigma_mid2;
1171        log_sigma_mid2 = log_sigma_mid1;
1172        qmid2 = qmid1;
1173        log_sigma_mid1 = log_sigma_lo + gfac*(log_sigma_up-log_sigma_lo);
1174        qmid1 = CalculateQualityFunction(exp(log_sigma_mid1),
1175                                         step_aff_x_L,
1176                                         step_aff_x_U,
1177                                         step_aff_s_L,
1178                                         step_aff_s_U,
1179                                         step_aff_y_c,
1180                                         step_aff_y_d,
1181                                         step_aff_z_L,
1182                                         step_aff_z_U,
1183                                         step_aff_v_L,
1184                                         step_aff_v_U,
1185                                         step_cen_x_L,
1186                                         step_cen_x_U,
1187                                         step_cen_s_L,
1188                                         step_cen_s_U,
1189                                         step_cen_y_c,
1190                                         step_cen_y_d,
1191                                         step_cen_z_L,
1192                                         step_cen_z_U,
1193                                         step_cen_v_L,
1194                                         step_cen_v_U);
1195      }
1196    }
1197
1198    Number q;
1199    if (qmid1 < qmid2) {
1200      log_sigma = log_sigma_mid1;
1201      q = qmid1;
1202    }
1203    else {
1204      log_sigma = log_sigma_mid2;
1205      q = qmid2;
1206    }
1207    if (log_sigma_up == log_sigma_up_in) {
1208      Number qtmp = CalculateQualityFunction(exp(log_sigma_up),
1209                                             step_aff_x_L,
1210                                             step_aff_x_U,
1211                                             step_aff_s_L,
1212                                             step_aff_s_U,
1213                                             step_aff_y_c,
1214                                             step_aff_y_d,
1215                                             step_aff_z_L,
1216                                             step_aff_z_U,
1217                                             step_aff_v_L,
1218                                             step_aff_v_U,
1219                                             step_cen_x_L,
1220                                             step_cen_x_U,
1221                                             step_cen_s_L,
1222                                             step_cen_s_U,
1223                                             step_cen_y_c,
1224                                             step_cen_y_d,
1225                                             step_cen_z_L,
1226                                             step_cen_z_U,
1227                                             step_cen_v_L,
1228                                             step_cen_v_U);
1229      if (qtmp < q) {
1230        log_sigma = log_sigma_up;
1231        q = qtmp;
1232      }
1233    }
1234    else if (log_sigma_lo == log_sigma_lo_in) {
1235      Number qtmp = CalculateQualityFunction(exp(log_sigma_lo),
1236                                             step_aff_x_L,
1237                                             step_aff_x_U,
1238                                             step_aff_s_L,
1239                                             step_aff_s_U,
1240                                             step_aff_y_c,
1241                                             step_aff_y_d,
1242                                             step_aff_z_L,
1243                                             step_aff_z_U,
1244                                             step_aff_v_L,
1245                                             step_aff_v_U,
1246                                             step_cen_x_L,
1247                                             step_cen_x_U,
1248                                             step_cen_s_L,
1249                                             step_cen_s_U,
1250                                             step_cen_y_c,
1251                                             step_cen_y_d,
1252                                             step_cen_z_L,
1253                                             step_cen_z_U,
1254                                             step_cen_v_L,
1255                                             step_cen_v_U);
1256      if (qtmp < q) {
1257        log_sigma = log_sigma_lo;
1258        q = qtmp;
1259      }
1260    }
1261
1262    return exp(log_sigma);
1263  }
1264  */
1265
1266
1267} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.