source: branches/dev/Algorithm/IpIpoptCalculatedQuantities.cpp @ 544

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

fixed calculation of 2-norm in CalculatedQuantities?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 108.5 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: IpIpoptCalculatedQuantities.cpp 544 2005-10-14 18:14:20Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9#include "IpIpoptCalculatedQuantities.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
21#include <limits>
22
23namespace Ipopt
24{
25#ifdef IP_DEBUG
26  static const Index dbg_verbosity = 0;
27#endif
28
29  IpoptCalculatedQuantities::IpoptCalculatedQuantities
30  (const SmartPtr<IpoptNLP>& ip_nlp,
31   const SmartPtr<IpoptData>& ip_data)
32      :
33      ip_nlp_(ip_nlp),
34      ip_data_(ip_data),
35
36      curr_slack_x_L_cache_(1),
37      curr_slack_x_U_cache_(1),
38      curr_slack_s_L_cache_(1),
39      curr_slack_s_U_cache_(1),
40      trial_slack_x_L_cache_(1),
41      trial_slack_x_U_cache_(1),
42      trial_slack_s_L_cache_(1),
43      trial_slack_s_U_cache_(1),
44      num_adjusted_slack_x_L_(0),
45      num_adjusted_slack_x_U_(0),
46      num_adjusted_slack_s_L_(0),
47      num_adjusted_slack_s_U_(0),
48
49      curr_f_cache_(2),
50      trial_f_cache_(5),
51      curr_grad_f_cache_(2),
52      trial_grad_f_cache_(1),
53
54      curr_barrier_obj_cache_(2),
55      trial_barrier_obj_cache_(5),
56      curr_grad_barrier_obj_x_cache_(1),
57      curr_grad_barrier_obj_s_cache_(1),
58      grad_kappa_times_damping_x_cache_(1),
59      grad_kappa_times_damping_s_cache_(1),
60
61      curr_c_cache_(1),
62      trial_c_cache_(2),
63      curr_d_cache_(1),
64      trial_d_cache_(2),
65      curr_d_minus_s_cache_(1),
66      trial_d_minus_s_cache_(1),
67      curr_jac_c_cache_(1),
68      trial_jac_c_cache_(1),
69      curr_jac_d_cache_(1),
70      trial_jac_d_cache_(1),
71      curr_jac_cT_times_vec_cache_(2),
72      trial_jac_cT_times_vec_cache_(1),
73      curr_jac_dT_times_vec_cache_(2),
74      trial_jac_dT_times_vec_cache_(1),
75      curr_jac_c_times_vec_cache_(1),
76      curr_jac_d_times_vec_cache_(1),
77      curr_constraint_violation_cache_(2),
78      trial_constraint_violation_cache_(5),
79      curr_nlp_constraint_violation_cache_(3),
80      unscaled_curr_nlp_constraint_violation_cache_(3),
81
82      curr_exact_hessian_cache_(1),
83
84      curr_grad_lag_x_cache_(1),
85      trial_grad_lag_x_cache_(1),
86      curr_grad_lag_s_cache_(1),
87      trial_grad_lag_s_cache_(1),
88      curr_grad_lag_with_damping_x_cache_(0),
89      curr_grad_lag_with_damping_s_cache_(0),
90      curr_compl_x_L_cache_(1),
91      curr_compl_x_U_cache_(1),
92      curr_compl_s_L_cache_(1),
93      curr_compl_s_U_cache_(1),
94      trial_compl_x_L_cache_(1),
95      trial_compl_x_U_cache_(1),
96      trial_compl_s_L_cache_(1),
97      trial_compl_s_U_cache_(1),
98      curr_relaxed_compl_x_L_cache_(1),
99      curr_relaxed_compl_x_U_cache_(1),
100      curr_relaxed_compl_s_L_cache_(1),
101      curr_relaxed_compl_s_U_cache_(1),
102      curr_primal_infeasibility_cache_(3),
103      trial_primal_infeasibility_cache_(3),
104      curr_dual_infeasibility_cache_(3),
105      trial_dual_infeasibility_cache_(3),
106      unscaled_curr_dual_infeasibility_cache_(3),
107      curr_complementarity_cache_(6),
108      trial_complementarity_cache_(6),
109      curr_centrality_measure_cache_(1),
110      curr_nlp_error_cache_(1),
111      unscaled_curr_nlp_error_cache_(1),
112      curr_barrier_error_cache_(1),
113      curr_primal_dual_system_error_cache_(1),
114      trial_primal_dual_system_error_cache_(3),
115
116      primal_frac_to_the_bound_cache_(5),
117      dual_frac_to_the_bound_cache_(5),
118
119      curr_sigma_x_cache_(1),
120      curr_sigma_s_cache_(1),
121
122      curr_avrg_compl_cache_(1),
123      trial_avrg_compl_cache_(1),
124      curr_gradBarrTDelta_cache_(1),
125
126      dampind_x_L_(NULL),
127      dampind_x_U_(NULL),
128      dampind_s_L_(NULL),
129      dampind_s_U_(NULL),
130
131      initialize_called_(false)
132  {
133    DBG_START_METH("IpoptCalculatedQuantities::IpoptCalculatedQuantities",
134                   dbg_verbosity);
135    DBG_ASSERT(IsValid(ip_nlp_) && IsValid(ip_data_));
136  }
137
138  IpoptCalculatedQuantities::~IpoptCalculatedQuantities()
139  {}
140
141  void IpoptCalculatedQuantities::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
142  {
143    roptions->SetRegisteringCategory("Convergence");
144    roptions->AddLowerBoundedNumberOption(
145      "s_max",
146      "Scaling threshold for the NLP error.",
147      0.0, true, 100.0,
148      "(See paragraph after Eqn. (6) in the implementation paper.)");
149
150    roptions->SetRegisteringCategory("NLP");
151    roptions->AddLowerBoundedNumberOption(
152      "kappa_d",
153      "Weight for linear damping term (to handle one-sided bounds).",
154      0.0, false, 1e-5,
155      "(see Section 3.7 in implementation paper.)");
156
157    roptions->SetRegisteringCategory("Line Search");
158    roptions->AddLowerBoundedNumberOption(
159      "slack_move",
160      "Correction size for very small slacks.",
161      0.0, false,
162      pow(std::numeric_limits<double>::epsilon(), 0.75),
163      "Due to numerical issues or the lack of an interior, the slack variables might "
164      "become very small.  If a slack becomes very small compared to machine "
165      "precision, the corresponding bound is moved slightly.  This parameter "
166      "determines how large the move should be.  Its default value is "
167      "mach_eps^{3/4}.  (See also end of Section 3.5 in implementation paper "
168      "- but actual implementation might be somewhat different.)");
169    roptions->SetRegisteringCategory("Convergence");
170    roptions->AddStringOption3(
171      "constraint_violation_norm_type",
172      "Norm to be used for the constraint violation.",
173      "1-norm",
174      "1-norm", "use the 1-norm",
175      "2-norm", "use the 2-norm",
176      "max-norm", "use the infinity norm",
177      "Determines which norm should be used when the algorithm computes the "
178      "constraint violation (for example, in the line search).");
179  }
180
181  bool IpoptCalculatedQuantities::Initialize(const Journalist& jnlst,
182      const OptionsList& options,
183      const std::string& prefix)
184  {
185    std::string svalue;
186    Index enum_int;
187
188    options.GetNumericValue("s_max", s_max_, prefix);
189    options.GetNumericValue("kappa_d", kappa_d_, prefix);
190    options.GetNumericValue("slack_move", slack_move_, prefix);
191    options.GetEnumValue("constraint_violation_norm_type", enum_int, prefix);
192    constr_viol_normtype_ = ENormType(enum_int);
193    // The following option is registered by OrigIpoptNLP
194    options.GetBoolValue("warm_start_same_structure",
195                         warm_start_same_structure_, prefix);
196
197    if (!warm_start_same_structure_) {
198      dampind_x_L_ = NULL;
199      dampind_x_U_ = NULL;
200      dampind_s_L_ = NULL;
201      dampind_s_U_ = NULL;
202
203      tmp_x_ = NULL;
204      tmp_s_ = NULL;
205      tmp_c_ = NULL;
206      tmp_d_ = NULL;
207      tmp_x_L_ = NULL;
208      tmp_x_U_ = NULL;
209      tmp_s_L_ = NULL;
210      tmp_s_U_ = NULL;
211    }
212
213    num_adjusted_slack_x_L_ = 0;
214    num_adjusted_slack_x_U_ = 0;
215    num_adjusted_slack_s_L_ = 0;
216    num_adjusted_slack_s_U_ = 0;
217
218    initialize_called_ = true;
219    return true;
220  }
221
222  ///////////////////////////////////////////////////////////////////////////
223  //                         Slack Calculations                            //
224  ///////////////////////////////////////////////////////////////////////////
225
226  SmartPtr<Vector>
227  IpoptCalculatedQuantities::CalcSlack_L(const Matrix& P,
228                                         const Vector& x,
229                                         const Vector& x_bound)
230  {
231    DBG_START_METH("IpoptCalculatedQuantities::CalcSlack_L",
232                   dbg_verbosity);
233    SmartPtr<Vector> result;
234    result = x_bound.MakeNew();
235    result->Copy(x_bound);
236    P.TransMultVector(1.0, x, -1.0, *result);
237    return result;
238  }
239
240  SmartPtr<Vector>
241  IpoptCalculatedQuantities::CalcSlack_U(const Matrix& P,
242                                         const Vector& x,
243                                         const Vector& x_bound)
244  {
245    DBG_START_METH("IpoptCalculatedQuantities::CalcSlack_U",
246                   dbg_verbosity);
247    SmartPtr<Vector> result;
248    result = x_bound.MakeNew();
249    result->Copy(x_bound);
250    P.TransMultVector(-1.0, x, 1.0, *result);
251    return result;
252  }
253
254  SmartPtr<const Vector> IpoptCalculatedQuantities::curr_slack_x_L()
255  {
256    DBG_START_METH("IpoptCalculatedQuantities::curr_slack_x_L()",
257                   dbg_verbosity);
258    SmartPtr<Vector> result;
259    SmartPtr<const Vector> x = ip_data_->curr()->x();
260    SmartPtr<const Vector> x_bound = ip_nlp_->x_L();
261    if (!curr_slack_x_L_cache_.GetCachedResult1Dep(result, *x)) {
262      if (!trial_slack_x_L_cache_.GetCachedResult1Dep(result, *x)) {
263        SmartPtr<const Matrix> P = ip_nlp_->Px_L();
264        DBG_PRINT_VECTOR(2,"x_L", *x_bound);
265        result = CalcSlack_L(*P, *x, *x_bound);
266        DBG_ASSERT(num_adjusted_slack_x_L_==0);
267        num_adjusted_slack_x_L_ =
268          CalculateSafeSlack(result, x_bound, x, ip_data_->curr()->z_L());
269      }
270      curr_slack_x_L_cache_.AddCachedResult1Dep(result, *x);
271    }
272    return ConstPtr(result);
273  }
274
275  SmartPtr<const Vector> IpoptCalculatedQuantities::curr_slack_x_U()
276  {
277    DBG_START_METH("IpoptCalculatedQuantities::curr_slack_x_U()",
278                   dbg_verbosity);
279
280    SmartPtr<Vector> result;
281    SmartPtr<const Vector> x = ip_data_->curr()->x();
282    SmartPtr<const Vector> x_bound = ip_nlp_->x_U();
283    if (!curr_slack_x_U_cache_.GetCachedResult1Dep(result, *x)) {
284      if (!trial_slack_x_U_cache_.GetCachedResult1Dep(result, *x)) {
285        SmartPtr<const Matrix> P = ip_nlp_->Px_U();
286        result = CalcSlack_U(*P, *x, *x_bound);
287        DBG_ASSERT(num_adjusted_slack_x_U_==0);
288        num_adjusted_slack_x_U_ =
289          CalculateSafeSlack(result, x_bound, x, ip_data_->curr()->z_U());
290      }
291      curr_slack_x_U_cache_.AddCachedResult1Dep(result, *x);
292    }
293    return ConstPtr(result);
294  }
295
296  SmartPtr<const Vector> IpoptCalculatedQuantities::curr_slack_s_L()
297  {
298    DBG_START_METH("IpoptCalculatedQuantities::curr_slack_s_L()",
299                   dbg_verbosity);
300
301    SmartPtr<Vector> result;
302    SmartPtr<const Vector> s = ip_data_->curr()->s();
303    SmartPtr<const Vector> s_bound = ip_nlp_->d_L();
304    if (!curr_slack_s_L_cache_.GetCachedResult1Dep(result, *s)) {
305      if (!trial_slack_s_L_cache_.GetCachedResult1Dep(result, *s)) {
306        SmartPtr<const Matrix> P = ip_nlp_->Pd_L();
307        result = CalcSlack_L(*P, *s, *s_bound);
308        DBG_ASSERT(num_adjusted_slack_s_L_==0);
309        num_adjusted_slack_s_L_ =
310          CalculateSafeSlack(result, s_bound, s, ip_data_->curr()->v_L());
311      }
312      curr_slack_s_L_cache_.AddCachedResult1Dep(result, *s);
313    }
314    return ConstPtr(result);
315  }
316
317  SmartPtr<const Vector> IpoptCalculatedQuantities::curr_slack_s_U()
318  {
319    DBG_START_METH("IpoptCalculatedQuantities::curr_slack_s_U()",
320                   dbg_verbosity);
321
322    SmartPtr<Vector> result;
323    SmartPtr<const Vector> s = ip_data_->curr()->s();
324    SmartPtr<const Vector> s_bound = ip_nlp_->d_U();
325    if (!curr_slack_s_U_cache_.GetCachedResult1Dep(result, *s)) {
326      if (!trial_slack_s_U_cache_.GetCachedResult1Dep(result, *s)) {
327        SmartPtr<const Matrix> P = ip_nlp_->Pd_U();
328        result = CalcSlack_U(*P, *s, *s_bound);
329        DBG_ASSERT(num_adjusted_slack_s_U_==0);
330        num_adjusted_slack_s_U_ =
331          CalculateSafeSlack(result, s_bound, s, ip_data_->curr()->v_U());
332        DBG_PRINT_VECTOR(2, "result", *result);
333        DBG_PRINT((1, "num_adjusted_slack_s_U = %d\n", num_adjusted_slack_s_U_));
334      }
335      curr_slack_s_U_cache_.AddCachedResult1Dep(result, *s);
336    }
337    return ConstPtr(result);
338  }
339
340  SmartPtr<const Vector> IpoptCalculatedQuantities::trial_slack_x_L()
341  {
342    DBG_START_METH("IpoptCalculatedQuantities::trial_slack_x_L()",
343                   dbg_verbosity);
344
345    num_adjusted_slack_x_L_ = 0;
346    SmartPtr<Vector> result;
347    SmartPtr<const Vector> x = ip_data_->trial()->x();
348    SmartPtr<const Vector> x_bound = ip_nlp_->x_L();
349    if (!trial_slack_x_L_cache_.GetCachedResult1Dep(result, *x)) {
350      if (!curr_slack_x_L_cache_.GetCachedResult1Dep(result, *x)) {
351        SmartPtr<const Matrix> P = ip_nlp_->Px_L();
352        result = CalcSlack_L(*P, *x, *x_bound);
353        DBG_ASSERT(num_adjusted_slack_x_L_==0);
354        num_adjusted_slack_x_L_ =
355          CalculateSafeSlack(result, x_bound, x, ip_data_->curr()->z_L());
356      }
357      trial_slack_x_L_cache_.AddCachedResult1Dep(result, *x);
358    }
359    return ConstPtr(result);
360  }
361
362  SmartPtr<const Vector> IpoptCalculatedQuantities::trial_slack_x_U()
363  {
364    DBG_START_METH("IpoptCalculatedQuantities::trial_slack_x_U()",
365                   dbg_verbosity);
366
367    num_adjusted_slack_x_U_ = 0;
368    SmartPtr<Vector> result;
369    SmartPtr<const Vector> x = ip_data_->trial()->x();
370    SmartPtr<const Vector> x_bound = ip_nlp_->x_U();
371    if (!trial_slack_x_U_cache_.GetCachedResult1Dep(result, *x)) {
372      if (!curr_slack_x_U_cache_.GetCachedResult1Dep(result, *x)) {
373        SmartPtr<const Matrix> P = ip_nlp_->Px_U();
374        result = CalcSlack_U(*P, *x, *x_bound);
375        DBG_ASSERT(num_adjusted_slack_x_U_==0);
376        num_adjusted_slack_x_U_ =
377          CalculateSafeSlack(result, x_bound, x, ip_data_->curr()->z_U());
378      }
379      trial_slack_x_U_cache_.AddCachedResult1Dep(result, *x);
380    }
381    return ConstPtr(result);
382  }
383
384  SmartPtr<const Vector> IpoptCalculatedQuantities::trial_slack_s_L()
385  {
386    DBG_START_METH("IpoptCalculatedQuantities::trial_slack_s_L()",
387                   dbg_verbosity);
388
389    num_adjusted_slack_s_L_ = 0;
390    SmartPtr<Vector> result;
391    SmartPtr<const Vector> s = ip_data_->trial()->s();
392    SmartPtr<const Vector> s_bound = ip_nlp_->d_L();
393    if (!trial_slack_s_L_cache_.GetCachedResult1Dep(result, *s)) {
394      if (!curr_slack_s_L_cache_.GetCachedResult1Dep(result, *s)) {
395        SmartPtr<const Matrix> P = ip_nlp_->Pd_L();
396        result = CalcSlack_L(*P, *s, *s_bound);
397        DBG_ASSERT(num_adjusted_slack_s_L_==0);
398        num_adjusted_slack_s_L_ =
399          CalculateSafeSlack(result, s_bound, s, ip_data_->curr()->v_L());
400      }
401      trial_slack_s_L_cache_.AddCachedResult1Dep(result, *s);
402    }
403    return ConstPtr(result);
404  }
405
406  SmartPtr<const Vector> IpoptCalculatedQuantities::trial_slack_s_U()
407  {
408    DBG_START_METH("IpoptCalculatedQuantities::trial_slack_s_U()",
409                   dbg_verbosity);
410
411    num_adjusted_slack_s_U_ = 0;
412    SmartPtr<Vector> result;
413    SmartPtr<const Vector> s = ip_data_->trial()->s();
414    SmartPtr<const Vector> s_bound = ip_nlp_->d_U();
415    if (!trial_slack_s_U_cache_.GetCachedResult1Dep(result, *s)) {
416      if (!curr_slack_s_U_cache_.GetCachedResult1Dep(result, *s)) {
417        SmartPtr<const Matrix> P = ip_nlp_->Pd_U();
418        result = CalcSlack_U(*P, *s, *s_bound);
419        DBG_ASSERT(num_adjusted_slack_s_U_==0);
420        num_adjusted_slack_s_U_ =
421          CalculateSafeSlack(result, s_bound, s, ip_data_->curr()->v_U());
422        DBG_PRINT((1, "num_adjusted_slack_s_U = %d\n", num_adjusted_slack_s_U_));
423        DBG_PRINT_VECTOR(2, "trial_slack_s_U", *result);
424      }
425      trial_slack_s_U_cache_.AddCachedResult1Dep(result, *s);
426    }
427    return ConstPtr(result);
428  }
429
430  Index IpoptCalculatedQuantities::
431  CalculateSafeSlack(SmartPtr<Vector>& slack,
432                     const SmartPtr<const Vector>& bound,
433                     const SmartPtr<const Vector>& curr_point,
434                     const SmartPtr<const Vector>& multiplier)
435  {
436    DBG_START_METH("IpoptCalculatedQuantities::CalculateSafeSlack", dbg_verbosity);
437    DBG_ASSERT(initialize_called_);
438    Index retval = 0;
439    if (slack->Dim() > 0) {
440      Number min_slack = slack->Min();
441      // TODO we need to make sure that this also works for non-monotone MUs
442      Number s_min = std::numeric_limits<Number>::epsilon()
443                     * Min(1., ip_data_->curr_mu());
444      DBG_PRINT((1,"s_min = %g, min_slack=%g\n", s_min, min_slack));
445      if (min_slack < s_min) {
446        // Need to correct the slacks and calculate new bounds...
447        SmartPtr<Vector> t = slack->MakeNew();
448        t->Copy(*slack);
449        t->AddScalar(-s_min);
450        t->ElementWiseSgn();
451
452        SmartPtr<Vector> zero_vec = t->MakeNew();
453        zero_vec->Set(0.0);
454        t->ElementWiseMin(*zero_vec);
455        t->Scal(-1.0);
456        retval = (Index)t->Asum();
457        DBG_PRINT((1,"Number of slack corrections = %d\n", retval));
458        DBG_PRINT_VECTOR(2, "t(sgn)", *t);
459
460        // ToDo AW: I added the follwing line b/c I found a case where
461        // slack was negative and this correction produced 0
462        slack->ElementWiseMax(*zero_vec);
463
464        SmartPtr<Vector> t2 = t->MakeNew();
465        t2->Set(ip_data_->curr_mu());
466        t2->ElementWiseDivide(*multiplier);
467
468        SmartPtr<Vector> s_min_vec = t2->MakeNew();
469        s_min_vec->Set(s_min);
470
471        t2->ElementWiseMax(*s_min_vec);
472        t2->Axpy(-1.0, *slack);
473        DBG_PRINT_VECTOR(2, "tw(smin,mu/mult)", *t2);
474
475        t->ElementWiseMultiply(*t2);
476        t->Axpy(1.0, *slack);
477
478        SmartPtr<Vector> t_max = t2;
479        t_max->Set(1.0);
480        SmartPtr<Vector> abs_bound = bound->MakeNew();
481        abs_bound->Copy(*bound);
482        abs_bound->ElementWiseAbs();
483        t_max->ElementWiseMax(*abs_bound);
484        DBG_PRINT_VECTOR(2, "t_max1", *t_max);
485        DBG_PRINT_VECTOR(2, "slack", *slack);
486        t_max->AddOneVector(1.0, *slack, slack_move_);
487        DBG_PRINT_VECTOR(2, "t_max2", *t_max);
488
489        t->ElementWiseMin(*t_max);
490        DBG_PRINT_VECTOR(2, "new_slack", *t);
491
492        slack = t;
493        return retval;
494      }
495    }
496
497    return retval;
498  }
499
500  Index
501  IpoptCalculatedQuantities::AdjustedTrialSlacks()
502  {
503    DBG_START_METH("IpoptCalculatedQuantities::AdjustedTrialSlacks()",
504                   dbg_verbosity);
505    Index result =  (num_adjusted_slack_x_L_ +
506                     num_adjusted_slack_x_U_ +
507                     num_adjusted_slack_s_L_ +
508                     num_adjusted_slack_s_U_);
509    DBG_PRINT((1,"result = %d\n", result));
510    return result;
511  }
512
513  void
514  IpoptCalculatedQuantities::ResetAdjustedTrialSlacks()
515  {
516    DBG_START_METH("IpoptCalculatedQuantities::ResetAdjustedTrialSlacks()",
517                   dbg_verbosity);
518    num_adjusted_slack_x_L_
519    = num_adjusted_slack_x_U_
520      = num_adjusted_slack_s_L_
521        = num_adjusted_slack_s_U_ = 0;
522  }
523
524  ///////////////////////////////////////////////////////////////////////////
525  //                          Objective Function                           //
526  ///////////////////////////////////////////////////////////////////////////
527
528  Number
529  IpoptCalculatedQuantities::curr_f()
530  {
531    DBG_START_METH("IpoptCalculatedQuantities::curr_f()",
532                   dbg_verbosity);
533    Number result;
534    SmartPtr<const Vector> x = ip_data_->curr()->x();
535    DBG_PRINT_VECTOR(2,"curr_x",*x);
536    DBG_PRINT((1, "curr_x tag = %d\n", x->GetTag()));
537
538    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
539    std::vector<const TaggedObject*> tdeps(1);
540    tdeps[0] = GetRawPtr(x);
541    std::vector<Number> sdeps(1);
542    if (objective_depends_on_mu) {
543      sdeps[0] = ip_data_->curr_mu();
544    }
545    else {
546      sdeps[0] = -1.;
547    }
548
549    if (!curr_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
550      if (!trial_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
551        DBG_PRINT((2,"evaluate curr f\n"));
552        if (objective_depends_on_mu) {
553          result = ip_nlp_->f(*x, ip_data_->curr_mu());
554        }
555        else {
556          result = ip_nlp_->f(*x);
557        }
558      }
559      curr_f_cache_.AddCachedResult(result, tdeps, sdeps);
560    }
561    DBG_PRINT((1,"result (curr_f) = %e\n", result));
562    return result;
563  }
564
565  Number
566  IpoptCalculatedQuantities::unscaled_curr_f()
567  {
568    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(curr_f());
569  }
570
571  Number
572  IpoptCalculatedQuantities::trial_f()
573  {
574    DBG_START_METH("IpoptCalculatedQuantities::trial_f()",
575                   dbg_verbosity);
576    Number result;
577    SmartPtr<const Vector> x = ip_data_->trial()->x();
578    DBG_PRINT_VECTOR(2,"trial_x",*x);
579    DBG_PRINT((1, "trial_x tag = %d\n", x->GetTag()));
580
581    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
582    std::vector<const TaggedObject*> tdeps(1);
583    tdeps[0] = GetRawPtr(x);
584    std::vector<Number> sdeps(1);
585    if (objective_depends_on_mu) {
586      sdeps[0] = ip_data_->curr_mu();
587    }
588    else {
589      sdeps[0] = -1.;
590    }
591
592    if (!trial_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
593      if (!curr_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
594        DBG_PRINT((2,"evaluate trial f\n"));
595        if (objective_depends_on_mu) {
596          result = ip_nlp_->f(*x, ip_data_->curr_mu());
597        }
598        else {
599          result = ip_nlp_->f(*x);
600        }
601      }
602      trial_f_cache_.AddCachedResult(result, tdeps, sdeps);
603    }
604    DBG_PRINT((1,"result (trial_f) = %e\n", result));
605    return result;
606  }
607
608  Number
609  IpoptCalculatedQuantities::unscaled_trial_f()
610  {
611    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(trial_f());
612  }
613
614  SmartPtr<const Vector>
615  IpoptCalculatedQuantities::curr_grad_f()
616  {
617    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_f()",
618                   dbg_verbosity);
619    SmartPtr<const Vector> result;
620    SmartPtr<const Vector> x = ip_data_->curr()->x();
621
622    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
623    std::vector<const TaggedObject*> tdeps(1);
624    tdeps[0] = GetRawPtr(x);
625    std::vector<Number> sdeps(1);
626    if (objective_depends_on_mu) {
627      sdeps[0] = ip_data_->curr_mu();
628    }
629    else {
630      sdeps[0] = -1.;
631    }
632
633    if (!curr_grad_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
634      if (!trial_grad_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
635        if (objective_depends_on_mu) {
636          result = ip_nlp_->grad_f(*x, ip_data_->curr_mu());
637        }
638        else {
639          result = ip_nlp_->grad_f(*x);
640        }
641      }
642      curr_grad_f_cache_.AddCachedResult(result, tdeps, sdeps);
643    }
644    return result;
645  }
646
647  SmartPtr<const Vector>
648  IpoptCalculatedQuantities::trial_grad_f()
649  {
650    DBG_START_METH("IpoptCalculatedQuantities::trial_grad_f()",
651                   dbg_verbosity);
652    SmartPtr<const Vector> result;
653    SmartPtr<const Vector> x = ip_data_->trial()->x();
654
655    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
656    std::vector<const TaggedObject*> tdeps(1);
657    tdeps[0] = GetRawPtr(x);
658    std::vector<Number> sdeps(1);
659    if (objective_depends_on_mu) {
660      sdeps[0] = ip_data_->curr_mu();
661    }
662    else {
663      sdeps[0] = -1.;
664    }
665
666    if (!trial_grad_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
667      if (!curr_grad_f_cache_.GetCachedResult(result, tdeps, sdeps)) {
668        if (objective_depends_on_mu) {
669          result = ip_nlp_->grad_f(*x, ip_data_->curr_mu());
670        }
671        else {
672          result = ip_nlp_->grad_f(*x);
673        }
674      }
675      trial_grad_f_cache_.AddCachedResult(result, tdeps, sdeps);
676    }
677    return result;
678  }
679
680  ///////////////////////////////////////////////////////////////////////////
681  //                    Barrier Objective Function                         //
682  ///////////////////////////////////////////////////////////////////////////
683  Number
684  IpoptCalculatedQuantities::CalcBarrierTerm(Number mu,
685      const Vector& slack_x_L,
686      const Vector& slack_x_U,
687      const Vector& slack_s_L,
688      const Vector& slack_s_U)
689  {
690    DBG_START_METH("IpoptCalculatedQuantities::CalcBarrierTerm",
691                   dbg_verbosity);
692    DBG_ASSERT(initialize_called_);
693
694    DBG_PRINT_VECTOR(2, "slack_x_L", slack_x_L);
695    DBG_PRINT_VECTOR(2, "slack_x_U", slack_x_U);
696    DBG_PRINT_VECTOR(2, "slack_s_L", slack_s_L);
697    DBG_PRINT_VECTOR(2, "slack_s_U", slack_s_U);
698
699    Number retval=0.;
700    retval += slack_x_L.SumLogs();
701    DBG_PRINT((1, "BarrierTerm after x_L = %25.16e\n", retval));
702    retval += slack_x_U.SumLogs();
703    DBG_PRINT((1, "BarrierTerm after x_U = %25.16e\n", retval));
704    retval += slack_s_L.SumLogs();
705    DBG_PRINT((1, "BarrierTerm after s_L = %25.16e\n", retval));
706    retval += slack_s_U.SumLogs();
707    DBG_PRINT((1, "BarrierTerm after s_U = %25.16e\n", retval));
708    retval *= -mu;
709
710    DBG_PRINT((1, "BarrierTerm without damping = %25.16e\n", retval));
711
712    // Include the linear damping term if kappa_d is nonzero.
713    if (kappa_d_>0) {
714      SmartPtr<const Vector> dampind_x_L;
715      SmartPtr<const Vector> dampind_x_U;
716      SmartPtr<const Vector> dampind_s_L;
717      SmartPtr<const Vector> dampind_s_U;
718      ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
719
720      Tmp_x_L().Copy(slack_x_L);
721      Tmp_x_L().ElementWiseMultiply(*dampind_x_L);
722      retval += kappa_d_ * mu * Tmp_x_L().Asum();
723      Tmp_x_U().Copy(slack_x_U);
724      Tmp_x_U().ElementWiseMultiply(*dampind_x_U);
725      retval += kappa_d_ * mu * Tmp_x_U().Asum();
726      Tmp_s_L().Copy(slack_s_L);
727      Tmp_s_L().ElementWiseMultiply(*dampind_s_L);
728      retval += kappa_d_ * mu * Tmp_s_L().Asum();
729      Tmp_s_U().Copy(slack_s_U);
730      Tmp_s_U().ElementWiseMultiply(*dampind_s_U);
731      retval += kappa_d_ * mu * Tmp_s_U().Asum();
732    }
733
734    DBG_PRINT((1, "BarrierTerm with damping = %25.16e\n", retval));
735
736    DBG_ASSERT(IsFiniteNumber(retval));
737    return retval;
738  }
739
740  Number
741  IpoptCalculatedQuantities::curr_barrier_obj()
742  {
743    DBG_START_METH("IpoptCalculatedQuantities::curr_barrier_obj()",
744                   dbg_verbosity);
745    Number result;
746
747    SmartPtr<const Vector> x = ip_data_->curr()->x();
748    SmartPtr<const Vector> s = ip_data_->curr()->s();
749    DBG_PRINT_VECTOR(2,"curr_x",*x);
750    DBG_PRINT_VECTOR(2,"curr_s",*s);
751    std::vector<const TaggedObject*> tdeps(2);
752    tdeps[0] = GetRawPtr(x);
753    tdeps[1] = GetRawPtr(s);
754
755    Number mu = ip_data_->curr_mu();
756    DBG_PRINT((1,"curr_mu=%e\n",mu));
757    std::vector<Number> sdeps(1);
758    sdeps[0] = mu;
759
760    if (!curr_barrier_obj_cache_.GetCachedResult(result, tdeps, sdeps)) {
761      if (!trial_barrier_obj_cache_.GetCachedResult(result, tdeps, sdeps)) {
762        result = curr_f();
763        DBG_PRINT((1,"curr_F=%e\n",result));
764        result += CalcBarrierTerm(mu,
765                                  *curr_slack_x_L(),
766                                  *curr_slack_x_U(),
767                                  *curr_slack_s_L(),
768                                  *curr_slack_s_U());
769      }
770      curr_barrier_obj_cache_.AddCachedResult(result, tdeps, sdeps);
771    }
772    DBG_ASSERT(IsFiniteNumber(result));
773    return result;
774  }
775
776  Number
777  IpoptCalculatedQuantities::trial_barrier_obj()
778  {
779    DBG_START_METH("IpoptCalculatedQuantities::trial_barrier_obj()",
780                   dbg_verbosity);
781    Number result;
782
783    SmartPtr<const Vector> x = ip_data_->trial()->x();
784    SmartPtr<const Vector> s = ip_data_->trial()->s();
785    DBG_PRINT_VECTOR(2,"trial_x",*x);
786    DBG_PRINT_VECTOR(2,"trial_s",*s);
787    std::vector<const TaggedObject*> tdeps(2);
788    tdeps[0] = GetRawPtr(x);
789    tdeps[1] = GetRawPtr(s);
790
791    Number mu = ip_data_->curr_mu();
792    DBG_PRINT((1,"trial_mu=%e\n",mu));
793    std::vector<Number> sdeps(1);
794    sdeps[0] = mu;
795
796    if (!trial_barrier_obj_cache_.GetCachedResult(result, tdeps, sdeps)) {
797      if (!curr_barrier_obj_cache_.GetCachedResult(result, tdeps, sdeps)) {
798        result = trial_f();
799        DBG_PRINT((1,"trial_F=%e\n",result));
800        DBG_PRINT_VECTOR(2, "trial_slack_s_U", *trial_slack_s_U());
801        result += CalcBarrierTerm(ip_data_->curr_mu(),
802                                  *trial_slack_x_L(),
803                                  *trial_slack_x_U(),
804                                  *trial_slack_s_L(),
805                                  *trial_slack_s_U());
806      }
807      trial_barrier_obj_cache_.AddCachedResult(result, tdeps, sdeps);
808    }
809    DBG_ASSERT(IsFiniteNumber(result));
810    return result;
811  }
812
813  SmartPtr<const Vector>
814  IpoptCalculatedQuantities::curr_grad_barrier_obj_x()
815  {
816    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_barrier_obj_x()",
817                   dbg_verbosity);
818    DBG_ASSERT(initialize_called_);
819    SmartPtr<const Vector> result;
820
821    SmartPtr<const Vector> x = ip_data_->curr()->x();
822    std::vector<const TaggedObject*> tdeps(1);
823    tdeps[0] = GetRawPtr(x);
824    Number mu = ip_data_->curr_mu();
825    std::vector<Number> sdeps(1);
826    sdeps[0] = mu;
827    DBG_PRINT((1,"curr_mu=%e\n",mu));
828
829    if (!curr_grad_barrier_obj_x_cache_.GetCachedResult(result, tdeps, sdeps)) {
830      SmartPtr<Vector> tmp1 = x->MakeNew();
831      tmp1->Copy(*curr_grad_f());
832
833      Tmp_x_L().Set(1.);
834      ip_nlp_->Px_L()->AddMSinvZ(-mu, *curr_slack_x_L(), Tmp_x_L(), *tmp1);
835
836      Tmp_x_U().Set(1.);
837      ip_nlp_->Px_U()->AddMSinvZ(mu, *curr_slack_x_U(), Tmp_x_U(), *tmp1);
838
839      DBG_PRINT_VECTOR(2, "Barrier_Grad_x without damping", *tmp1);
840
841      // Take care of linear damping terms
842      if (kappa_d_>0.) {
843        SmartPtr<const Vector> dampind_x_L;
844        SmartPtr<const Vector> dampind_x_U;
845        SmartPtr<const Vector> dampind_s_L;
846        SmartPtr<const Vector> dampind_s_U;
847        ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
848
849        DBG_PRINT((1, "kappa_d*mu = %e\n", kappa_d_*mu));
850        DBG_PRINT_VECTOR(2, "dampind_x_L", *dampind_x_L);
851        ip_nlp_->Px_L()->MultVector(kappa_d_*mu, *dampind_x_L, 1., *tmp1);
852        ip_nlp_->Px_U()->MultVector(-kappa_d_*mu, *dampind_x_U, 1., *tmp1);
853      }
854
855      DBG_PRINT_VECTOR(2, "Barrier_Grad_x with damping", *tmp1);
856
857      result = ConstPtr(tmp1);
858
859      curr_grad_barrier_obj_x_cache_.AddCachedResult(result, tdeps, sdeps);
860    }
861
862    return result;
863  }
864
865  SmartPtr<const Vector>
866  IpoptCalculatedQuantities::grad_kappa_times_damping_x()
867  {
868    DBG_START_METH("IpoptCalculatedQuantities::grad_kappa_times_damping_x()",
869                   dbg_verbosity);
870    DBG_ASSERT(initialize_called_);
871    SmartPtr<const Vector> result;
872    SmartPtr<const Vector> x = ip_data_->curr()->x();
873
874    std::vector<const TaggedObject*> tdeps(2);
875    tdeps[0] = GetRawPtr(ip_nlp_->Px_L());
876    tdeps[1] = GetRawPtr(ip_nlp_->Px_U());
877    std::vector<Number> sdeps(1);
878    sdeps[0] = kappa_d_;
879    if (!grad_kappa_times_damping_x_cache_.GetCachedResult(result, tdeps, sdeps)) {
880      SmartPtr<Vector> tmp1 = x->MakeNew();
881      if (kappa_d_>0.) {
882        SmartPtr<const Vector> dampind_x_L;
883        SmartPtr<const Vector> dampind_x_U;
884        SmartPtr<const Vector> dampind_s_L;
885        SmartPtr<const Vector> dampind_s_U;
886        ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
887
888        ip_nlp_->Px_L()->MultVector(kappa_d_, *dampind_x_L, 0., *tmp1);
889        ip_nlp_->Px_U()->MultVector(-kappa_d_, *dampind_x_U, 1., *tmp1);
890      }
891      else {
892        tmp1->Set(0.);
893      }
894      result = ConstPtr(tmp1);
895
896      grad_kappa_times_damping_x_cache_.AddCachedResult(result, tdeps, sdeps);
897    }
898
899    return result;
900  }
901
902  SmartPtr<const Vector>
903  IpoptCalculatedQuantities::curr_grad_barrier_obj_s()
904  {
905    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_barrier_obj_s()",
906                   dbg_verbosity);
907    DBG_ASSERT(initialize_called_);
908    SmartPtr<const Vector> result;
909
910    SmartPtr<const Vector> s = ip_data_->curr()->s();
911    std::vector<const TaggedObject*> tdeps(1);
912    tdeps[0] = GetRawPtr(s);
913    Number mu = ip_data_->curr_mu();
914    std::vector<Number> sdeps(1);
915    sdeps[0] = mu;
916    DBG_PRINT((1,"curr_mu=%e\n",mu));
917
918    if (!curr_grad_barrier_obj_s_cache_.GetCachedResult(result, tdeps, sdeps)) {
919      SmartPtr<Vector> tmp1 = s->MakeNew();
920
921      Tmp_s_L().Set(-mu);
922      Tmp_s_L().ElementWiseDivide(*curr_slack_s_L());
923      ip_nlp_->Pd_L()->MultVector(1., Tmp_s_L(), 0., *tmp1);
924
925      Tmp_s_U().Set(1.);
926      ip_nlp_->Pd_U()->AddMSinvZ(mu, *curr_slack_s_U(), Tmp_s_U(), *tmp1);
927
928      DBG_PRINT_VECTOR(2, "Barrier_Grad_s without damping", *tmp1);
929
930      // Take care of linear damping terms
931      if (kappa_d_>0.) {
932        SmartPtr<const Vector> dampind_x_L;
933        SmartPtr<const Vector> dampind_x_U;
934        SmartPtr<const Vector> dampind_s_L;
935        SmartPtr<const Vector> dampind_s_U;
936        ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
937
938        DBG_PRINT((1, "kappa_d*mu = %e\n", kappa_d_*mu));
939        DBG_PRINT_VECTOR(2, "dampind_s_L", *dampind_s_L);
940        DBG_PRINT_VECTOR(2, "dampind_s_U", *dampind_s_U);
941        ip_nlp_->Pd_L()->MultVector(kappa_d_*mu, *dampind_s_L, 1., *tmp1);
942        ip_nlp_->Pd_U()->MultVector(-kappa_d_*mu, *dampind_s_U, 1., *tmp1);
943      }
944
945      DBG_PRINT_VECTOR(2, "Barrier_Grad_s with damping", *tmp1);
946
947      result = ConstPtr(tmp1);
948
949      curr_grad_barrier_obj_s_cache_.AddCachedResult(result, tdeps, sdeps);
950    }
951
952    return result;
953  }
954
955  SmartPtr<const Vector>
956  IpoptCalculatedQuantities::grad_kappa_times_damping_s()
957  {
958    DBG_START_METH("IpoptCalculatedQuantities::grad_kappa_times_damping_s()",
959                   dbg_verbosity);
960    DBG_ASSERT(initialize_called_);
961    SmartPtr<const Vector> result;
962    SmartPtr<const Vector> s = ip_data_->curr()->s();
963
964    std::vector<const TaggedObject*> tdeps(2);
965    tdeps[0] = GetRawPtr(ip_nlp_->Pd_L());
966    tdeps[1] = GetRawPtr(ip_nlp_->Pd_U());
967    std::vector<Number> sdeps(1);
968    sdeps[0] = kappa_d_;
969    if (!grad_kappa_times_damping_s_cache_.GetCachedResult(result, tdeps, sdeps)) {
970      SmartPtr<Vector> tmp1 = s->MakeNew();
971      if (kappa_d_>0.) {
972        SmartPtr<const Vector> dampind_x_L;
973        SmartPtr<const Vector> dampind_x_U;
974        SmartPtr<const Vector> dampind_s_L;
975        SmartPtr<const Vector> dampind_s_U;
976        ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
977
978        ip_nlp_->Pd_L()->MultVector(kappa_d_, *dampind_s_L, 0., *tmp1);
979        ip_nlp_->Pd_U()->MultVector(-kappa_d_, *dampind_s_U, 1., *tmp1);
980      }
981      else {
982        tmp1->Set(0.);
983      }
984      result = ConstPtr(tmp1);
985
986      grad_kappa_times_damping_s_cache_.AddCachedResult(result, tdeps, sdeps);
987    }
988
989    return result;
990  }
991
992  void
993  IpoptCalculatedQuantities::ComputeDampingIndicators(
994    SmartPtr<const Vector>& dampind_x_L,
995    SmartPtr<const Vector>& dampind_x_U,
996    SmartPtr<const Vector>& dampind_s_L,
997    SmartPtr<const Vector>& dampind_s_U)
998  {
999    DBG_START_METH("IpoptCalculatedQuantities::ComputeDampingFilters()",
1000                   dbg_verbosity);
1001
1002    // Assume that all indicators have to be computed if one of the
1003    // SmartPtrs is still zero.
1004    if (IsNull(dampind_x_L_)) {
1005      // First for x
1006      Tmp_x_L().Set(1.0);
1007      ip_nlp_->Px_L()->MultVector(1.0, Tmp_x_L(), 0.0, Tmp_x());
1008      Tmp_x_U().Set(1.0);
1009      ip_nlp_->Px_U()->MultVector(-1.0, Tmp_x_U(), 1.0, Tmp_x());
1010
1011      dampind_x_L_ = ip_nlp_->x_L()->MakeNew();
1012      ip_nlp_->Px_L()->TransMultVector(1.0, Tmp_x(), 0.0, *dampind_x_L_);
1013
1014      dampind_x_U_ = ip_nlp_->x_U()->MakeNew();
1015      ip_nlp_->Px_U()->TransMultVector(-1.0, Tmp_x(), 0.0, *dampind_x_U_);
1016
1017      // Now for s
1018      Tmp_s_L().Set(1.0);
1019      ip_nlp_->Pd_L()->MultVector(1.0, Tmp_s_L(), 0.0, Tmp_s());
1020      Tmp_s_U().Set(1.0);
1021      ip_nlp_->Pd_U()->MultVector(-1.0, Tmp_s_U(), 1.0, Tmp_s());
1022
1023      dampind_s_L_ = ip_nlp_->d_L()->MakeNew();
1024      ip_nlp_->Pd_L()->TransMultVector(1.0, Tmp_s(), 0.0, *dampind_s_L_);
1025
1026      dampind_s_U_ = ip_nlp_->d_U()->MakeNew();
1027      ip_nlp_->Pd_U()->TransMultVector(-1.0, Tmp_s(), 0.0, *dampind_s_U_);
1028
1029      DBG_PRINT_VECTOR(2, "dampind_x_L_", *dampind_x_L_);
1030      DBG_PRINT_VECTOR(2, "dampind_x_U_", *dampind_x_U_);
1031      DBG_PRINT_VECTOR(2, "dampind_s_L_", *dampind_s_L_);
1032      DBG_PRINT_VECTOR(2, "dampind_s_U_", *dampind_s_U_);
1033    }
1034
1035    dampind_x_L = ConstPtr(dampind_x_L_);
1036    dampind_x_U = ConstPtr(dampind_x_U_);
1037    dampind_s_L = ConstPtr(dampind_s_L_);
1038    dampind_s_U = ConstPtr(dampind_s_U_);
1039  }
1040
1041  ///////////////////////////////////////////////////////////////////////////
1042  //                                Constraints                            //
1043  ///////////////////////////////////////////////////////////////////////////
1044
1045  SmartPtr<const Vector>
1046  IpoptCalculatedQuantities::curr_c()
1047  {
1048    DBG_START_METH("IpoptCalculatedQuantities::curr_c()",
1049                   dbg_verbosity);
1050    SmartPtr<const Vector> result;
1051    SmartPtr<const Vector> x = ip_data_->curr()->x();
1052
1053    if (!curr_c_cache_.GetCachedResult1Dep(result, *x)) {
1054      if (!trial_c_cache_.GetCachedResult1Dep(result, *x)) {
1055        result = ip_nlp_->c(*x);
1056      }
1057      curr_c_cache_.AddCachedResult1Dep(result, *x);
1058    }
1059    return result;
1060  }
1061
1062  SmartPtr<const Vector>
1063  IpoptCalculatedQuantities::unscaled_curr_c()
1064  {
1065    return ip_nlp_->NLP_scaling()->unapply_vector_scaling_c(curr_c());
1066  }
1067
1068  SmartPtr<const Vector>
1069  IpoptCalculatedQuantities::trial_c()
1070  {
1071    DBG_START_METH("IpoptCalculatedQuantities::trial_c()",
1072                   dbg_verbosity);
1073    SmartPtr<const Vector> result;
1074    SmartPtr<const Vector> x = ip_data_->trial()->x();
1075
1076    if (!trial_c_cache_.GetCachedResult1Dep(result, *x)) {
1077      if (!curr_c_cache_.GetCachedResult1Dep(result, *x)) {
1078        result = ip_nlp_->c(*x);
1079      }
1080      trial_c_cache_.AddCachedResult1Dep(result, *x);
1081    }
1082    return result;
1083  }
1084
1085  SmartPtr<const Vector>
1086  IpoptCalculatedQuantities::curr_d()
1087  {
1088    DBG_START_METH("IpoptCalculatedQuantities::curr_d()",
1089                   dbg_verbosity);
1090    SmartPtr<const Vector> result;
1091    SmartPtr<const Vector> x = ip_data_->curr()->x();
1092
1093    if (!curr_d_cache_.GetCachedResult1Dep(result, *x)) {
1094      if (!trial_d_cache_.GetCachedResult1Dep(result, *x)) {
1095        result = ip_nlp_->d(*x);
1096      }
1097      curr_d_cache_.AddCachedResult1Dep(result, *x);
1098    }
1099    return result;
1100  }
1101
1102  SmartPtr<const Vector>
1103  IpoptCalculatedQuantities::unscaled_curr_d()
1104  {
1105    return ip_nlp_->NLP_scaling()->unapply_vector_scaling_d(curr_d());
1106  }
1107
1108  SmartPtr<const Vector>
1109  IpoptCalculatedQuantities::trial_d()
1110  {
1111    DBG_START_METH("IpoptCalculatedQuantities::trial_d()",
1112                   dbg_verbosity);
1113    SmartPtr<const Vector> result;
1114    SmartPtr<const Vector> x = ip_data_->trial()->x();
1115
1116    if (!trial_d_cache_.GetCachedResult1Dep(result, *x)) {
1117      if (!curr_d_cache_.GetCachedResult1Dep(result, *x)) {
1118        result = ip_nlp_->d(*x);
1119      }
1120      trial_d_cache_.AddCachedResult1Dep(result, *x);
1121    }
1122    return result;
1123  }
1124
1125  SmartPtr<const Vector>
1126  IpoptCalculatedQuantities::curr_d_minus_s()
1127  {
1128    DBG_START_METH("IpoptCalculatedQuantities::curr_d_minus_s()",
1129                   dbg_verbosity);
1130    SmartPtr<const Vector> result;
1131
1132    SmartPtr<const Vector> x = ip_data_->curr()->x();
1133    SmartPtr<const Vector> s = ip_data_->curr()->s();
1134
1135    if (!curr_d_minus_s_cache_.GetCachedResult2Dep(result, *x, *s)) {
1136      if (!trial_d_minus_s_cache_.GetCachedResult2Dep(result, *x, *s)) {
1137        SmartPtr<Vector> tmp = s->MakeNew();
1138        tmp->AddTwoVectors(1., *curr_d(), -1., *s, 0.);
1139        result = ConstPtr(tmp);
1140      }
1141      curr_d_minus_s_cache_.AddCachedResult2Dep(result, *x, *s);
1142    }
1143
1144    return result;
1145  }
1146
1147  SmartPtr<const Vector>
1148  IpoptCalculatedQuantities::trial_d_minus_s()
1149  {
1150    DBG_START_METH("IpoptCalculatedQuantities::trial_d_minus_s()",
1151                   dbg_verbosity);
1152    SmartPtr<const Vector> result;
1153
1154    SmartPtr<const Vector> x = ip_data_->trial()->x();
1155    SmartPtr<const Vector> s = ip_data_->trial()->s();
1156
1157    if (!trial_d_minus_s_cache_.GetCachedResult2Dep(result, *x, *s)) {
1158      if (!curr_d_minus_s_cache_.GetCachedResult2Dep(result, *x, *s)) {
1159        SmartPtr<Vector> tmp = s->MakeNew();
1160        tmp->AddTwoVectors(1., *trial_d(), -1., *s, 0.);
1161        result = ConstPtr(tmp);
1162      }
1163      trial_d_minus_s_cache_.AddCachedResult2Dep(result, *x, *s);
1164    }
1165
1166    return result;
1167  }
1168
1169  SmartPtr<const Matrix>
1170  IpoptCalculatedQuantities::curr_jac_c()
1171  {
1172    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_c()",
1173                   dbg_verbosity);
1174    SmartPtr<const Matrix> result;
1175    SmartPtr<const Vector> x = ip_data_->curr()->x();
1176
1177    if (!curr_jac_c_cache_.GetCachedResult1Dep(result, *x)) {
1178      if (!trial_jac_c_cache_.GetCachedResult1Dep(result, *x)) {
1179        result = ip_nlp_->jac_c(*x);
1180      }
1181      curr_jac_c_cache_.AddCachedResult1Dep(result, *x);
1182    }
1183    return result;
1184  }
1185
1186  SmartPtr<const Matrix>
1187  IpoptCalculatedQuantities::trial_jac_c()
1188  {
1189    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_c()",
1190                   dbg_verbosity);
1191    SmartPtr<const Matrix> result;
1192    SmartPtr<const Vector> x = ip_data_->trial()->x();
1193
1194    if (!trial_jac_c_cache_.GetCachedResult1Dep(result, *x)) {
1195      if (!curr_jac_c_cache_.GetCachedResult1Dep(result, *x)) {
1196        result = ip_nlp_->jac_c(*x);
1197      }
1198      trial_jac_c_cache_.AddCachedResult1Dep(result, *x);
1199    }
1200    return result;
1201  }
1202
1203  SmartPtr<const Matrix>
1204  IpoptCalculatedQuantities::curr_jac_d()
1205  {
1206    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_d()",
1207                   dbg_verbosity);
1208    SmartPtr<const Matrix> result;
1209    SmartPtr<const Vector> x = ip_data_->curr()->x();
1210
1211    if (!curr_jac_d_cache_.GetCachedResult1Dep(result, *x)) {
1212      if (!trial_jac_d_cache_.GetCachedResult1Dep(result, *x)) {
1213        result = ip_nlp_->jac_d(*x);
1214      }
1215      curr_jac_d_cache_.AddCachedResult1Dep(result, *x);
1216    }
1217    return result;
1218  }
1219
1220  SmartPtr<const Matrix>
1221  IpoptCalculatedQuantities::trial_jac_d()
1222  {
1223    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_d()",
1224                   dbg_verbosity);
1225    SmartPtr<const Matrix> result;
1226    SmartPtr<const Vector> x = ip_data_->trial()->x();
1227
1228    if (!trial_jac_d_cache_.GetCachedResult1Dep(result, *x)) {
1229      if (!curr_jac_d_cache_.GetCachedResult1Dep(result, *x)) {
1230        result = ip_nlp_->jac_d(*x);
1231      }
1232      trial_jac_d_cache_.AddCachedResult1Dep(result, *x);
1233    }
1234    return result;
1235  }
1236
1237  SmartPtr<const Vector>
1238  IpoptCalculatedQuantities::curr_jac_c_times_vec(const Vector& vec)
1239  {
1240    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_c_times_vec",
1241                   dbg_verbosity);
1242    SmartPtr<const Vector> result;
1243    SmartPtr<const Vector> x = ip_data_->curr()->x();
1244
1245    if (!curr_jac_c_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1246      SmartPtr<Vector> tmp = ip_data_->curr()->y_c()->MakeNew();
1247      curr_jac_c()->MultVector(1.0, vec, 0., *tmp);
1248      result = ConstPtr(tmp);
1249      curr_jac_c_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1250    }
1251
1252    return result;
1253  }
1254
1255  SmartPtr<const Vector>
1256  IpoptCalculatedQuantities::curr_jac_d_times_vec(const Vector& vec)
1257  {
1258    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_d_times_vec()",
1259                   dbg_verbosity);
1260    SmartPtr<const Vector> result;
1261    SmartPtr<const Vector> x = ip_data_->curr()->x();
1262
1263    if (!curr_jac_d_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1264      SmartPtr<Vector> tmp = ip_data_->curr()->s()->MakeNew();
1265      curr_jac_d()->MultVector(1.0, vec, 0., *tmp);
1266      result = ConstPtr(tmp);
1267      curr_jac_d_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1268    }
1269
1270    return result;
1271  }
1272
1273  SmartPtr<const Vector>
1274  IpoptCalculatedQuantities::curr_jac_cT_times_curr_y_c()
1275  {
1276    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_cT_times_curr_y_c()",
1277                   dbg_verbosity);
1278    return curr_jac_cT_times_vec(*ip_data_->curr()->y_c());
1279  }
1280
1281  SmartPtr<const Vector>
1282  IpoptCalculatedQuantities::trial_jac_cT_times_trial_y_c()
1283  {
1284    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_cT_times_trial_y_c()",
1285                   dbg_verbosity);
1286    return trial_jac_cT_times_vec(*ip_data_->trial()->y_c());
1287  }
1288
1289  SmartPtr<const Vector>
1290  IpoptCalculatedQuantities::curr_jac_dT_times_curr_y_d()
1291  {
1292    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_dT_times_curr_y_d()",
1293                   dbg_verbosity);
1294    return curr_jac_dT_times_vec(*ip_data_->curr()->y_d());
1295  }
1296
1297  SmartPtr<const Vector>
1298  IpoptCalculatedQuantities::trial_jac_dT_times_trial_y_d()
1299  {
1300    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_dT_times_trial_y_d()",
1301                   dbg_verbosity);
1302    return trial_jac_dT_times_vec(*ip_data_->trial()->y_d());
1303  }
1304
1305  SmartPtr<const Vector>
1306  IpoptCalculatedQuantities::curr_jac_cT_times_vec(const Vector& vec)
1307  {
1308    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_cT_times_vec",
1309                   dbg_verbosity);
1310    SmartPtr<const Vector> result;
1311    SmartPtr<const Vector> x = ip_data_->curr()->x();
1312
1313    if (!curr_jac_cT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1314      if (!trial_jac_cT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1315        SmartPtr<Vector> tmp = x->MakeNew();
1316        curr_jac_c()->TransMultVector(1.0, vec, 0., *tmp);
1317        result = ConstPtr(tmp);
1318      }
1319      curr_jac_cT_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1320    }
1321
1322    return result;
1323  }
1324
1325  SmartPtr<const Vector>
1326  IpoptCalculatedQuantities::trial_jac_cT_times_vec(const Vector& vec)
1327  {
1328    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_cT_times_vec",
1329                   dbg_verbosity);
1330    SmartPtr<const Vector> result;
1331    SmartPtr<const Vector> x = ip_data_->trial()->x();
1332
1333    if (!trial_jac_cT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1334      if (!curr_jac_cT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1335        SmartPtr<Vector> tmp = x->MakeNew();
1336        trial_jac_c()->TransMultVector(1.0, vec, 0., *tmp);
1337        result = ConstPtr(tmp);
1338      }
1339      trial_jac_cT_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1340    }
1341
1342    return result;
1343  }
1344
1345  SmartPtr<const Vector>
1346  IpoptCalculatedQuantities::curr_jac_dT_times_vec(const Vector& vec)
1347  {
1348    DBG_START_METH("IpoptCalculatedQuantities::curr_jac_dT_times_vec()",
1349                   dbg_verbosity);
1350    SmartPtr<const Vector> result;
1351    SmartPtr<const Vector> x = ip_data_->curr()->x();
1352
1353    if (!curr_jac_dT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1354      if (!trial_jac_dT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1355        SmartPtr<Vector> tmp = x->MakeNew();
1356        curr_jac_d()->TransMultVector(1.0, vec, 0., *tmp);
1357        result = ConstPtr(tmp);
1358      }
1359      curr_jac_dT_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1360    }
1361
1362    return result;
1363  }
1364
1365  SmartPtr<const Vector>
1366  IpoptCalculatedQuantities::trial_jac_dT_times_vec(const Vector& vec)
1367  {
1368    DBG_START_METH("IpoptCalculatedQuantities::trial_jac_dT_times_vec()",
1369                   dbg_verbosity);
1370    SmartPtr<const Vector> result;
1371    SmartPtr<const Vector> x = ip_data_->trial()->x();
1372
1373    if (!trial_jac_dT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1374      if (!curr_jac_dT_times_vec_cache_.GetCachedResult2Dep(result, *x, vec)) {
1375        SmartPtr<Vector> tmp = x->MakeNew();
1376        trial_jac_d()->TransMultVector(1.0, vec, 0., *tmp);
1377        result = ConstPtr(tmp);
1378      }
1379      trial_jac_dT_times_vec_cache_.AddCachedResult2Dep(result, *x, vec);
1380    }
1381
1382    return result;
1383  }
1384
1385  Number
1386  IpoptCalculatedQuantities::curr_constraint_violation()
1387  {
1388    DBG_START_METH("IpoptCalculatedQuantities::curr_constraint_violation()",
1389                   dbg_verbosity);
1390    return curr_primal_infeasibility(constr_viol_normtype_);
1391  }
1392
1393  Number
1394  IpoptCalculatedQuantities::trial_constraint_violation()
1395  {
1396    DBG_START_METH("IpoptCalculatedQuantities::trial_constraint_violation()",
1397                   dbg_verbosity);
1398    return trial_primal_infeasibility(constr_viol_normtype_);
1399  }
1400
1401  Number
1402  IpoptCalculatedQuantities::curr_nlp_constraint_violation
1403  (ENormType NormType)
1404  {
1405    DBG_START_METH("IpoptCalculatedQuantities::curr_nlp_constraint_violation()",
1406                   dbg_verbosity);
1407    Number result;
1408
1409    SmartPtr<const Vector> x = ip_data_->curr()->x();
1410
1411    std::vector<const TaggedObject*> deps(1);
1412    deps[0] = GetRawPtr(x);
1413    std::vector<Number> sdeps(1);
1414    sdeps[0] = (Number)NormType;
1415
1416    if (!curr_nlp_constraint_violation_cache_.GetCachedResult(result, deps, sdeps)) {
1417      SmartPtr<const Vector> c = curr_c();
1418      SmartPtr<const Vector> d = curr_d();
1419
1420      SmartPtr<Vector> d_viol_L = ip_nlp_->d_L()->MakeNewCopy();
1421      ip_nlp_->Pd_L()->TransMultVector(-1., *d, 1., *d_viol_L);
1422      SmartPtr<Vector> tmp = d_viol_L->MakeNew();
1423      tmp->Set(0.);
1424      d_viol_L->ElementWiseMax(*tmp);
1425      DBG_PRINT_VECTOR(2, "d_viol_L", *d_viol_L);
1426
1427      SmartPtr<Vector> d_viol_U = ip_nlp_->d_U()->MakeNewCopy();
1428      ip_nlp_->Pd_U()->TransMultVector(-1., *d, 1., *d_viol_U);
1429      tmp = d_viol_U->MakeNew();
1430      tmp->Set(0.);
1431      d_viol_U->ElementWiseMin(*tmp);
1432      DBG_PRINT_VECTOR(2, "d_viol_U", *d_viol_U);
1433
1434      std::vector<SmartPtr<const Vector> > vecs(3);
1435      vecs[0] = c;
1436      vecs[1] = GetRawPtr(d_viol_L);
1437      vecs[2] = GetRawPtr(d_viol_U);
1438      result = CalcNormOfType(NormType, vecs);
1439      curr_nlp_constraint_violation_cache_.AddCachedResult(result, deps, sdeps);
1440    }
1441
1442    return result;
1443  }
1444
1445  Number
1446  IpoptCalculatedQuantities::unscaled_curr_nlp_constraint_violation
1447  (ENormType NormType)
1448  {
1449    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_nlp_constraint_violation()",
1450                   dbg_verbosity);
1451    Number result;
1452
1453    SmartPtr<const Vector> x = ip_data_->curr()->x();
1454
1455    std::vector<const TaggedObject*> deps(1);
1456    deps[0] = GetRawPtr(x);
1457    std::vector<Number> sdeps(1);
1458    sdeps[0] = (Number)NormType;
1459
1460    if (!unscaled_curr_nlp_constraint_violation_cache_.GetCachedResult(result, deps, sdeps)) {
1461      SmartPtr<const Vector> c = unscaled_curr_c();
1462
1463      SmartPtr<const Vector> d = curr_d();
1464
1465      SmartPtr<Vector> d_viol_L = ip_nlp_->d_L()->MakeNewCopy();
1466      ip_nlp_->Pd_L()->TransMultVector(-1., *d, 1., *d_viol_L);
1467      SmartPtr<Vector> tmp = d_viol_L->MakeNew();
1468      tmp->Set(0.);
1469      d_viol_L->ElementWiseMax(*tmp);
1470      DBG_PRINT_VECTOR(2, "d_viol_L", *d_viol_L);
1471
1472      SmartPtr<Vector> d_viol_U = ip_nlp_->d_U()->MakeNewCopy();
1473      ip_nlp_->Pd_U()->TransMultVector(-1., *d, 1., *d_viol_U);
1474      tmp = d_viol_U->MakeNew();
1475      tmp->Set(0.);
1476      d_viol_U->ElementWiseMin(*tmp);
1477      DBG_PRINT_VECTOR(2, "d_viol_U", *d_viol_U);
1478
1479      std::vector<SmartPtr<const Vector> > vecs(3);
1480      vecs[0] = c;
1481      vecs[1] = GetRawPtr(d_viol_L);
1482      vecs[2] = GetRawPtr(d_viol_U);
1483      result = CalcNormOfType(NormType, vecs);
1484      unscaled_curr_nlp_constraint_violation_cache_.AddCachedResult(result, deps, sdeps);
1485    }
1486
1487    return result;
1488  }
1489
1490  ///////////////////////////////////////////////////////////////////////////
1491  //                Exact Hessian using second derivatives                 //
1492  ///////////////////////////////////////////////////////////////////////////
1493
1494  SmartPtr<const SymMatrix>
1495  IpoptCalculatedQuantities::curr_exact_hessian()
1496  {
1497    DBG_START_METH("IpoptCalculatedQuantities::curr_exact_hessian()",
1498                   dbg_verbosity);
1499
1500    SmartPtr<const SymMatrix> result;
1501    SmartPtr<const Vector> x = ip_data_->curr()->x();
1502    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
1503    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
1504
1505    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
1506    std::vector<const TaggedObject*> tdeps(3);
1507    tdeps[0] = GetRawPtr(x);
1508    tdeps[1] = GetRawPtr(y_c);
1509    tdeps[2] = GetRawPtr(y_d);
1510    std::vector<Number> sdeps(1);
1511    if (objective_depends_on_mu) {
1512      sdeps[0] = ip_data_->curr_mu();
1513    }
1514    else {
1515      sdeps[0] = -1.;
1516    }
1517
1518    if (!curr_exact_hessian_cache_.GetCachedResult(result, tdeps, sdeps)) {
1519      if (objective_depends_on_mu) {
1520        result = ip_nlp_->h(*x, 1.0, *y_c, *y_d, ip_data_->curr_mu());
1521      }
1522      else {
1523        result = ip_nlp_->h(*x, 1.0, *y_c, *y_d);
1524      }
1525      curr_exact_hessian_cache_.AddCachedResult(result, tdeps, sdeps);
1526    }
1527
1528    return result;
1529  }
1530
1531  ///////////////////////////////////////////////////////////////////////////
1532  //                             Zero Hessian                              //
1533  ///////////////////////////////////////////////////////////////////////////
1534
1535  SmartPtr<const SymMatrix>
1536  IpoptCalculatedQuantities::zero_hessian()
1537  {
1538    DBG_START_METH("IpoptCalculatedQuantities::zero_hessian()",
1539                   dbg_verbosity);
1540
1541    SmartPtr<const Vector> x = ip_data_->curr()->x();
1542    Tmp_c().Set(0.);
1543    Tmp_d().Set(0.);
1544
1545    SmartPtr<const SymMatrix> h ;
1546    bool objective_depends_on_mu = ip_nlp_->objective_depends_on_mu();
1547    if (objective_depends_on_mu) {
1548      h = ip_nlp_->h(*x, 0.0, Tmp_c(), Tmp_d(), 0.);
1549    }
1550    else {
1551      h = ip_nlp_->h(*x, 0.0, Tmp_c(), Tmp_d());
1552    }
1553
1554    DBG_PRINT_MATRIX(2, "zero_hessian", *h);
1555
1556    return h;
1557  }
1558
1559  ///////////////////////////////////////////////////////////////////////////
1560  //                  Optimality Error and its components                  //
1561  ///////////////////////////////////////////////////////////////////////////
1562
1563  SmartPtr<const Vector>
1564  IpoptCalculatedQuantities::curr_grad_lag_x()
1565  {
1566    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_lag_x()",
1567                   dbg_verbosity);
1568    SmartPtr<const Vector> result;
1569
1570    SmartPtr<const Vector> x = ip_data_->curr()->x();
1571    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
1572    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
1573    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
1574    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
1575
1576    std::vector<const TaggedObject*> deps(5);
1577    deps[0] = GetRawPtr(x);
1578    deps[1] = GetRawPtr(y_c);
1579    deps[2] = GetRawPtr(y_d);
1580    deps[3] = GetRawPtr(z_L);
1581    deps[4] = GetRawPtr(z_U);
1582
1583    if (!curr_grad_lag_x_cache_.GetCachedResult(result, deps)) {
1584      if (!trial_grad_lag_x_cache_.GetCachedResult(result, deps)) {
1585        SmartPtr<Vector> tmp = x->MakeNew();
1586        DBG_PRINT_VECTOR(2,"curr_grad_f",*curr_grad_f());
1587        tmp->Copy(*curr_grad_f());
1588        tmp->AddTwoVectors(1., *curr_jac_cT_times_curr_y_c(),
1589                           1., *curr_jac_dT_times_curr_y_d(), 1.);
1590        DBG_PRINT_VECTOR(2,"jac_cT*y_c",*curr_jac_cT_times_curr_y_c());
1591        DBG_PRINT_VECTOR(2,"jac_dT*y_d",*curr_jac_dT_times_curr_y_d());
1592        ip_nlp_->Px_L()->MultVector(-1., *z_L, 1., *tmp);
1593        ip_nlp_->Px_U()->MultVector(1., *z_U, 1., *tmp);
1594        result = ConstPtr(tmp);
1595      }
1596      curr_grad_lag_x_cache_.AddCachedResult(result, deps);
1597    }
1598
1599    return result;
1600  }
1601
1602  SmartPtr<const Vector>
1603  IpoptCalculatedQuantities::trial_grad_lag_x()
1604  {
1605    DBG_START_METH("IpoptCalculatedQuantities::trial_grad_lag_x()",
1606                   dbg_verbosity);
1607    SmartPtr<const Vector> result;
1608
1609    SmartPtr<const Vector> x = ip_data_->trial()->x();
1610    SmartPtr<const Vector> y_c = ip_data_->trial()->y_c();
1611    SmartPtr<const Vector> y_d = ip_data_->trial()->y_d();
1612    SmartPtr<const Vector> z_L = ip_data_->trial()->z_L();
1613    SmartPtr<const Vector> z_U = ip_data_->trial()->z_U();
1614
1615    std::vector<const TaggedObject*> deps(5);
1616    deps[0] = GetRawPtr(x);
1617    deps[1] = GetRawPtr(y_c);
1618    deps[2] = GetRawPtr(y_d);
1619    deps[3] = GetRawPtr(z_L);
1620    deps[4] = GetRawPtr(z_U);
1621
1622    if (!trial_grad_lag_x_cache_.GetCachedResult(result, deps)) {
1623      if (!curr_grad_lag_x_cache_.GetCachedResult(result, deps)) {
1624        SmartPtr<Vector> tmp = x->MakeNew();
1625        DBG_PRINT_VECTOR(2,"trial_grad_f",*trial_grad_f());
1626        tmp->Copy(*trial_grad_f());
1627        tmp->AddTwoVectors(1., *trial_jac_cT_times_trial_y_c(),
1628                           1., *trial_jac_dT_times_trial_y_d(), 1.);
1629        ip_nlp_->Px_L()->MultVector(-1., *z_L, 1., *tmp);
1630        ip_nlp_->Px_U()->MultVector(1., *z_U, 1., *tmp);
1631        result = ConstPtr(tmp);
1632      }
1633      trial_grad_lag_x_cache_.AddCachedResult(result, deps);
1634    }
1635
1636    return result;
1637  }
1638
1639  SmartPtr<const Vector>
1640  IpoptCalculatedQuantities::curr_grad_lag_s()
1641  {
1642    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_lag_s()",
1643                   dbg_verbosity);
1644    SmartPtr<const Vector> result;
1645
1646    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
1647    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
1648    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
1649
1650    std::vector<const TaggedObject*> deps(3);
1651    deps[0] = GetRawPtr(y_d);
1652    deps[1] = GetRawPtr(v_L);
1653    deps[2] = GetRawPtr(v_U);
1654
1655    if (!curr_grad_lag_s_cache_.GetCachedResult(result, deps)) {
1656      if (!trial_grad_lag_s_cache_.GetCachedResult(result, deps)) {
1657        SmartPtr<Vector> tmp = y_d->MakeNew();
1658        ip_nlp_->Pd_U()->MultVector(1., *v_U, 0., *tmp);
1659        ip_nlp_->Pd_L()->MultVector(-1., *v_L, 1., *tmp);
1660        tmp->Axpy(-1., *y_d);
1661        result = ConstPtr(tmp);
1662      }
1663      curr_grad_lag_s_cache_.AddCachedResult(result, deps);
1664    }
1665
1666    return result;
1667  }
1668
1669  SmartPtr<const Vector>
1670  IpoptCalculatedQuantities::trial_grad_lag_s()
1671  {
1672    DBG_START_METH("IpoptCalculatedQuantities::trial_grad_lag_s()",
1673                   dbg_verbosity);
1674    SmartPtr<const Vector> result;
1675
1676    SmartPtr<const Vector> y_d = ip_data_->trial()->y_d();
1677    SmartPtr<const Vector> v_L = ip_data_->trial()->v_L();
1678    SmartPtr<const Vector> v_U = ip_data_->trial()->v_U();
1679
1680    std::vector<const TaggedObject*> deps(3);
1681    deps[0] = GetRawPtr(y_d);
1682    deps[1] = GetRawPtr(v_L);
1683    deps[2] = GetRawPtr(v_U);
1684
1685    if (!trial_grad_lag_s_cache_.GetCachedResult(result, deps)) {
1686      if (!curr_grad_lag_s_cache_.GetCachedResult(result, deps)) {
1687        SmartPtr<Vector> tmp = y_d->MakeNew();
1688        ip_nlp_->Pd_U()->MultVector(1., *v_U, 0., *tmp);
1689        ip_nlp_->Pd_L()->MultVector(-1., *v_L, 1., *tmp);
1690        tmp->Axpy(-1., *y_d);
1691        result = ConstPtr(tmp);
1692      }
1693      trial_grad_lag_s_cache_.AddCachedResult(result, deps);
1694    }
1695
1696    return result;
1697  }
1698
1699  SmartPtr<const Vector>
1700  IpoptCalculatedQuantities::curr_grad_lag_with_damping_x()
1701  {
1702    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_lag_with_damping_x()",
1703                   dbg_verbosity);
1704
1705    /* If no damping is used, just return the gradient of the regular
1706       Lagrangian function */
1707    if (kappa_d_==0.) {
1708      return curr_grad_lag_x();
1709    }
1710
1711    SmartPtr<const Vector> result;
1712
1713    SmartPtr<const Vector> x = ip_data_->curr()->x();
1714    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
1715    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
1716    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
1717    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
1718    Number mu = ip_data_->curr_mu();
1719
1720    std::vector<const TaggedObject*> deps(5);
1721    deps[0] = GetRawPtr(x);
1722    deps[1] = GetRawPtr(y_c);
1723    deps[2] = GetRawPtr(y_d);
1724    deps[3] = GetRawPtr(z_L);
1725    deps[4] = GetRawPtr(z_U);
1726    std::vector<Number> sdeps(1);
1727    sdeps[0] = mu;
1728
1729    if (!curr_grad_lag_with_damping_x_cache_.GetCachedResult(result, deps, sdeps)) {
1730      SmartPtr<Vector> tmp = x->MakeNew();
1731      tmp->Copy(*curr_grad_lag_x());
1732
1733      SmartPtr<const Vector> dampind_x_L;
1734      SmartPtr<const Vector> dampind_x_U;
1735      SmartPtr<const Vector> dampind_s_L;
1736      SmartPtr<const Vector> dampind_s_U;
1737      ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
1738
1739      ip_nlp_->Px_L()->MultVector(kappa_d_*mu, *dampind_x_L, 1., *tmp);
1740      ip_nlp_->Px_U()->MultVector(-kappa_d_*mu, *dampind_x_U, 1., *tmp);
1741
1742      result = ConstPtr(tmp);
1743      curr_grad_lag_with_damping_x_cache_.AddCachedResult(result, deps, sdeps);
1744    }
1745
1746    return result;
1747  }
1748
1749  SmartPtr<const Vector>
1750  IpoptCalculatedQuantities::curr_grad_lag_with_damping_s()
1751  {
1752    DBG_START_METH("IpoptCalculatedQuantities::curr_grad_lag_with_damping_s()",
1753                   dbg_verbosity);
1754
1755    /* If no damping is used, just return the gradient of the regular
1756       Lagrangian function */
1757    if (kappa_d_==0.) {
1758      return curr_grad_lag_s();
1759    }
1760
1761    SmartPtr<const Vector> result;
1762
1763    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
1764    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
1765    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
1766    Number mu = ip_data_->curr_mu();
1767
1768    std::vector<const TaggedObject*> deps(3);
1769    deps[0] = GetRawPtr(y_d);
1770    deps[1] = GetRawPtr(v_L);
1771    deps[2] = GetRawPtr(v_U);
1772    std::vector<Number> sdeps(1);
1773    sdeps[0] = mu;
1774
1775    if (!curr_grad_lag_with_damping_s_cache_.GetCachedResult(result, deps, sdeps)) {
1776      SmartPtr<Vector> tmp = y_d->MakeNew();
1777      tmp->Copy(*curr_grad_lag_s());
1778
1779      SmartPtr<const Vector> dampind_x_L;
1780      SmartPtr<const Vector> dampind_x_U;
1781      SmartPtr<const Vector> dampind_s_L;
1782      SmartPtr<const Vector> dampind_s_U;
1783      ComputeDampingIndicators(dampind_x_L, dampind_x_U, dampind_s_L, dampind_s_U);
1784
1785      ip_nlp_->Pd_L()->MultVector(kappa_d_*mu, *dampind_s_L, 1., *tmp);
1786      ip_nlp_->Pd_U()->MultVector(-kappa_d_*mu, *dampind_s_U, 1., *tmp);
1787
1788      result = ConstPtr(tmp);
1789      curr_grad_lag_with_damping_s_cache_.AddCachedResult(result, deps, sdeps);
1790    }
1791
1792    return result;
1793  }
1794
1795  SmartPtr<const Vector>
1796  IpoptCalculatedQuantities::CalcCompl(const Vector& slack,
1797                                       const Vector& mult)
1798  {
1799    DBG_START_METH("IpoptCalculatedQuantities::CalcCompl()",
1800                   dbg_verbosity);
1801    SmartPtr<Vector> result = slack.MakeNew();
1802    result->Copy(slack);
1803    result->ElementWiseMultiply(mult);
1804    return ConstPtr(result);
1805  }
1806
1807  SmartPtr<const Vector>
1808  IpoptCalculatedQuantities::curr_compl_x_L()
1809  {
1810    DBG_START_METH("IpoptCalculatedQuantities::curr_compl_x_L()",
1811                   dbg_verbosity);
1812    SmartPtr<const Vector> result;
1813
1814    SmartPtr<const Vector> slack = curr_slack_x_L();
1815    SmartPtr<const Vector> mult = ip_data_->curr()->z_L();
1816    DBG_PRINT_VECTOR(2, "slack_x_L", *slack);
1817    DBG_PRINT_VECTOR(2, "z_L", *mult);
1818
1819    if (!curr_compl_x_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1820      if (!trial_compl_x_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1821        result = CalcCompl(*slack, *mult);
1822      }
1823      curr_compl_x_L_cache_.AddCachedResult2Dep(result, *slack, *mult);
1824    }
1825    return result;
1826  }
1827
1828  SmartPtr<const Vector>
1829  IpoptCalculatedQuantities::trial_compl_x_L()
1830  {
1831    DBG_START_METH("IpoptCalculatedQuantities::trial_compl_x_L()",
1832                   dbg_verbosity);
1833    SmartPtr<const Vector> result;
1834
1835    SmartPtr<const Vector> slack = trial_slack_x_L();
1836    SmartPtr<const Vector> mult = ip_data_->trial()->z_L();
1837    DBG_PRINT_VECTOR(2, "slack_x_L", *slack);
1838    DBG_PRINT_VECTOR(2, "z_L", *mult);
1839
1840    if (!trial_compl_x_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1841      if (!curr_compl_x_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1842        result = CalcCompl(*slack, *mult);
1843      }
1844      trial_compl_x_L_cache_.AddCachedResult2Dep(result, *slack, *mult);
1845    }
1846    return result;
1847  }
1848
1849  SmartPtr<const Vector>
1850  IpoptCalculatedQuantities::curr_compl_x_U()
1851  {
1852    DBG_START_METH("IpoptCalculatedQuantities::curr_compl_x_U()",
1853                   dbg_verbosity);
1854    SmartPtr<const Vector> result;
1855
1856    SmartPtr<const Vector> slack = curr_slack_x_U();
1857    SmartPtr<const Vector> mult = ip_data_->curr()->z_U();
1858
1859    if (!curr_compl_x_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1860      if (!trial_compl_x_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1861        result = CalcCompl(*slack, *mult);
1862      }
1863      curr_compl_x_U_cache_.AddCachedResult2Dep(result, *slack, *mult);
1864    }
1865    return result;
1866  }
1867
1868  SmartPtr<const Vector>
1869  IpoptCalculatedQuantities::trial_compl_x_U()
1870  {
1871    DBG_START_METH("IpoptCalculatedQuantities::trial_compl_x_U()",
1872                   dbg_verbosity);
1873    SmartPtr<const Vector> result;
1874
1875    SmartPtr<const Vector> slack = trial_slack_x_U();
1876    SmartPtr<const Vector> mult = ip_data_->trial()->z_U();
1877
1878    if (!trial_compl_x_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1879      if (!curr_compl_x_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1880        result = CalcCompl(*slack, *mult);
1881      }
1882      trial_compl_x_U_cache_.AddCachedResult2Dep(result, *slack, *mult);
1883    }
1884    return result;
1885  }
1886
1887  SmartPtr<const Vector>
1888  IpoptCalculatedQuantities::curr_compl_s_L()
1889  {
1890    DBG_START_METH("IpoptCalculatedQuantities::curr_compl_s_L()",
1891                   dbg_verbosity);
1892    SmartPtr<const Vector> result;
1893
1894    SmartPtr<const Vector> slack = curr_slack_s_L();
1895    SmartPtr<const Vector> mult = ip_data_->curr()->v_L();
1896
1897    if (!curr_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1898      if (!trial_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1899        result = CalcCompl(*slack, *mult);
1900      }
1901      curr_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult);
1902    }
1903    return result;
1904  }
1905
1906  SmartPtr<const Vector>
1907  IpoptCalculatedQuantities::trial_compl_s_L()
1908  {
1909    DBG_START_METH("IpoptCalculatedQuantities::trial_compl_s_L()",
1910                   dbg_verbosity);
1911    SmartPtr<const Vector> result;
1912
1913    SmartPtr<const Vector> slack = trial_slack_s_L();
1914    SmartPtr<const Vector> mult = ip_data_->trial()->v_L();
1915
1916    if (!trial_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1917      if (!curr_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1918        result = CalcCompl(*slack, *mult);
1919      }
1920      trial_compl_s_L_cache_.GetCachedResult2Dep(result, *slack, *mult);
1921    }
1922    return result;
1923  }
1924
1925  SmartPtr<const Vector>
1926  IpoptCalculatedQuantities::curr_compl_s_U()
1927  {
1928    DBG_START_METH("IpoptCalculatedQuantities::curr_compl_s_U()",
1929                   dbg_verbosity);
1930    SmartPtr<const Vector> result;
1931
1932    SmartPtr<const Vector> slack = curr_slack_s_U();
1933    SmartPtr<const Vector> mult = ip_data_->curr()->v_U();
1934
1935    if (!curr_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1936      if (!trial_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1937        result = CalcCompl(*slack, *mult);
1938      }
1939      curr_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult);
1940    }
1941    return result;
1942  }
1943
1944  SmartPtr<const Vector>
1945  IpoptCalculatedQuantities::trial_compl_s_U()
1946  {
1947    DBG_START_METH("IpoptCalculatedQuantities::trial_compl_s_U()",
1948                   dbg_verbosity);
1949    SmartPtr<const Vector> result;
1950
1951    SmartPtr<const Vector> slack = trial_slack_s_U();
1952    SmartPtr<const Vector> mult = ip_data_->trial()->v_U();
1953
1954    if (!trial_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1955      if (!curr_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult)) {
1956        result = CalcCompl(*slack, *mult);
1957      }
1958      trial_compl_s_U_cache_.GetCachedResult2Dep(result, *slack, *mult);
1959    }
1960    return result;
1961  }
1962
1963  SmartPtr<const Vector>
1964  IpoptCalculatedQuantities::curr_relaxed_compl_x_L()
1965  {
1966    DBG_START_METH("IpoptCalculatedQuantities::curr_relaxed_compl_x_L()",
1967                   dbg_verbosity);
1968    SmartPtr<const Vector> result;
1969
1970    SmartPtr<const Vector> slack = curr_slack_x_L();
1971    SmartPtr<const Vector> mult = ip_data_->curr()->z_L();
1972    std::vector<const TaggedObject*> tdeps(2);
1973    tdeps[0] = GetRawPtr(slack);
1974    tdeps[1] = GetRawPtr(mult);
1975
1976    Number mu = ip_data_->curr_mu();
1977    std::vector<Number> sdeps(1);
1978    sdeps[0] = mu;
1979
1980    if (!curr_relaxed_compl_x_L_cache_.GetCachedResult(result, tdeps, sdeps)) {
1981      SmartPtr<Vector> tmp = slack->MakeNew();
1982      tmp->Copy(*curr_compl_x_L());
1983      tmp->AddScalar(-mu);
1984      result = ConstPtr(tmp);
1985      curr_relaxed_compl_x_L_cache_.AddCachedResult(result, tdeps, sdeps);
1986    }
1987    return result;
1988  }
1989
1990  SmartPtr<const Vector>
1991  IpoptCalculatedQuantities::curr_relaxed_compl_x_U()
1992  {
1993    DBG_START_METH("IpoptCalculatedQuantities::curr_relaxed_compl_x_U()",
1994                   dbg_verbosity);
1995    SmartPtr<const Vector> result;
1996
1997    SmartPtr<const Vector> slack = curr_slack_x_U();
1998    SmartPtr<const Vector> mult = ip_data_->curr()->z_U();
1999    std::vector<const TaggedObject*> tdeps(2);
2000    tdeps[0] = GetRawPtr(slack);
2001    tdeps[1] = GetRawPtr(mult);
2002
2003    Number mu = ip_data_->curr_mu();
2004    std::vector<Number> sdeps(1);
2005    sdeps[0] = mu;
2006
2007    if (!curr_relaxed_compl_x_U_cache_.GetCachedResult(result, tdeps, sdeps)) {
2008      SmartPtr<Vector> tmp = slack->MakeNew();
2009      tmp->Copy(*curr_compl_x_U());
2010      tmp->AddScalar(-mu);
2011      result = ConstPtr(tmp);
2012      curr_relaxed_compl_x_U_cache_.AddCachedResult(result, tdeps, sdeps);
2013    }
2014    return result;
2015  }
2016
2017  SmartPtr<const Vector>
2018  IpoptCalculatedQuantities::curr_relaxed_compl_s_L()
2019  {
2020    DBG_START_METH("IpoptCalculatedQuantities::curr_relaxed_compl_s_L()",
2021                   dbg_verbosity);
2022    SmartPtr<const Vector> result;
2023
2024    SmartPtr<const Vector> slack = curr_slack_s_L();
2025    SmartPtr<const Vector> mult = ip_data_->curr()->v_L();
2026    std::vector<const TaggedObject*> tdeps(2);
2027    tdeps[0] = GetRawPtr(slack);
2028    tdeps[1] = GetRawPtr(mult);
2029
2030    Number mu = ip_data_->curr_mu();
2031    std::vector<Number> sdeps(1);
2032    sdeps[0] = mu;
2033
2034    if (!curr_relaxed_compl_s_L_cache_.GetCachedResult(result, tdeps, sdeps)) {
2035      SmartPtr<Vector> tmp = slack->MakeNew();
2036      tmp->Copy(*curr_compl_s_L());
2037      tmp->AddScalar(-mu);
2038      result = ConstPtr(tmp);
2039      curr_relaxed_compl_s_L_cache_.AddCachedResult(result, tdeps, sdeps);
2040    }
2041    return result;
2042  }
2043
2044  SmartPtr<const Vector>
2045  IpoptCalculatedQuantities::curr_relaxed_compl_s_U()
2046  {
2047    DBG_START_METH("IpoptCalculatedQuantities::curr_relaxed_compl_s_U()",
2048                   dbg_verbosity);
2049    SmartPtr<const Vector> result;
2050
2051    SmartPtr<const Vector> slack = curr_slack_s_U();
2052    SmartPtr<const Vector> mult = ip_data_->curr()->v_U();
2053    std::vector<const TaggedObject*> tdeps(2);
2054    tdeps[0] = GetRawPtr(slack);
2055    tdeps[1] = GetRawPtr(mult);
2056
2057    Number mu = ip_data_->curr_mu();
2058    std::vector<Number> sdeps(1);
2059    sdeps[0] = mu;
2060
2061    if (!curr_relaxed_compl_s_U_cache_.GetCachedResult(result, tdeps, sdeps)) {
2062      SmartPtr<Vector> tmp = slack->MakeNew();
2063      tmp->Copy(*curr_compl_s_U());
2064      tmp->AddScalar(-mu);
2065      result = ConstPtr(tmp);
2066      curr_relaxed_compl_s_U_cache_.AddCachedResult(result, tdeps, sdeps);
2067    }
2068    return result;
2069  }
2070
2071  Number
2072  IpoptCalculatedQuantities::CalcNormOfType
2073  (ENormType NormType,
2074   const Vector& vec1, const Vector& vec2)
2075  {
2076    switch (NormType) {
2077      case NORM_1 :
2078      return vec1.Asum() + vec2.Asum();
2079      case NORM_2 :
2080      return sqrt(pow(vec1.Nrm2(),2) + pow(vec2.Nrm2(),2));
2081      case NORM_MAX :
2082      return Max(vec1.Amax(), vec2.Amax());
2083      default:
2084      DBG_ASSERT("Unknown NormType.");
2085      return 0.0;
2086    }
2087  }
2088
2089  Number
2090  IpoptCalculatedQuantities::CalcNormOfType
2091  (ENormType NormType,
2092   std::vector<SmartPtr<const Vector> > vecs)
2093  {
2094    Number result=0.;
2095
2096    switch (NormType) {
2097      case NORM_1 :
2098      for (Index i=0; i<(Index)vecs.size(); i++) {
2099        result += vecs[i]->Asum();
2100      }
2101      break;
2102      case NORM_2 :
2103      for (Index i=0; i<(Index)vecs.size(); i++) {
2104        Number nrm = vecs[i]->Nrm2();
2105        result += nrm*nrm;
2106      }
2107      result = sqrt(result);
2108      break;
2109      case NORM_MAX :
2110      for (Index i=0; i<(Index)vecs.size(); i++) {
2111        result = Max(result, vecs[i]->Amax());
2112      }
2113      break;
2114      default:
2115      DBG_ASSERT("Unknown NormType.");
2116    }
2117
2118    return result;
2119  }
2120
2121  Number
2122  IpoptCalculatedQuantities::curr_primal_infeasibility
2123  (ENormType NormType)
2124  {
2125    DBG_START_METH("IpoptCalculatedQuantities::curr_primal_infeasibility()",
2126                   dbg_verbosity);
2127    Number result;
2128
2129    SmartPtr<const Vector> x = ip_data_->curr()->x();
2130    SmartPtr<const Vector> s = ip_data_->curr()->s();
2131
2132    DBG_PRINT_VECTOR(2, "x to eval", *x);
2133    DBG_PRINT_VECTOR(2, "s to eval", *s);
2134    DBG_PRINT((1,"NormType = %d\n", NormType))
2135
2136    std::vector<const TaggedObject*> deps(2);
2137    deps[0] = GetRawPtr(x);
2138    deps[1] = GetRawPtr(s);
2139    std::vector<Number> sdeps(1);
2140    sdeps[0] = (Number)NormType;
2141
2142    if (!curr_primal_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2143      if (!trial_primal_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2144        DBG_PRINT((1,"Recomputing recomputing infeasibility.\n"));
2145        SmartPtr<const Vector> c = curr_c();
2146        SmartPtr<const Vector> d_minus_s = curr_d_minus_s();
2147
2148        DBG_PRINT_VECTOR(2,"c", *c);
2149        DBG_PRINT_VECTOR(2,"d_minus_s", *d_minus_s);
2150
2151        result = CalcNormOfType(NormType, *c, *d_minus_s);
2152
2153      }
2154      curr_primal_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
2155    }
2156
2157    DBG_PRINT((1,"result = %e\n",result));
2158    return result;
2159  }
2160
2161  Number
2162  IpoptCalculatedQuantities::trial_primal_infeasibility
2163  (ENormType NormType)
2164  {
2165    DBG_START_METH("IpoptCalculatedQuantities::trial_primal_infeasibility()",
2166                   dbg_verbosity);
2167    Number result;
2168
2169    SmartPtr<const Vector> x = ip_data_->trial()->x();
2170    SmartPtr<const Vector> s = ip_data_->trial()->s();
2171
2172    DBG_PRINT_VECTOR(2, "x to eval", *x);
2173    DBG_PRINT_VECTOR(2, "s to eval", *s);
2174    DBG_PRINT((1,"NormType = %d\n", NormType))
2175
2176    std::vector<const TaggedObject*> deps(2);
2177    deps[0] = GetRawPtr(x);
2178    deps[1] = GetRawPtr(s);
2179    std::vector<Number> sdeps(1);
2180    sdeps[0] = (Number)NormType;
2181
2182    if (!trial_primal_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2183      if (!curr_primal_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2184        DBG_PRINT((1,"Recomputing recomputing infeasibility.\n"));
2185        SmartPtr<const Vector> c = trial_c();
2186        SmartPtr<const Vector> d_minus_s = trial_d_minus_s();
2187
2188        DBG_PRINT_VECTOR(2,"c", *c);
2189        DBG_PRINT_VECTOR(2,"d_minus_s", *d_minus_s);
2190
2191        result = CalcNormOfType(NormType, *c, *d_minus_s);
2192      }
2193      trial_primal_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
2194    }
2195
2196    DBG_PRINT((1,"result = %e\n",result));
2197    return result;
2198  }
2199
2200  Number
2201  IpoptCalculatedQuantities::curr_dual_infeasibility
2202  (ENormType NormType)
2203  {
2204    DBG_START_METH("IpoptCalculatedQuantities::curr_dual_infeasibility()",
2205                   dbg_verbosity);
2206    Number result;
2207
2208    SmartPtr<const Vector> x = ip_data_->curr()->x();
2209    SmartPtr<const Vector> s = ip_data_->curr()->s();
2210    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2211    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2212    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2213    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2214    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2215    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2216
2217    std::vector<const TaggedObject*> deps(8);
2218    deps[0] = GetRawPtr(x);
2219    deps[1] = GetRawPtr(s);
2220    deps[2] = GetRawPtr(y_c);
2221    deps[3] = GetRawPtr(y_d);
2222    deps[4] = GetRawPtr(z_L);
2223    deps[5] = GetRawPtr(z_U);
2224    deps[6] = GetRawPtr(v_L);
2225    deps[7] = GetRawPtr(v_U);
2226    std::vector<Number> sdeps(1);
2227    sdeps[0] = (Number)NormType;
2228
2229    if (!curr_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2230      if (!trial_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2231        SmartPtr<const Vector> grad_lag_x = curr_grad_lag_x();
2232        SmartPtr<const Vector> grad_lag_s = curr_grad_lag_s();
2233        DBG_PRINT_VECTOR(2,"grad_lag_x", *grad_lag_x);
2234        DBG_PRINT_VECTOR(2,"grad_lag_s", *grad_lag_s);
2235
2236        result = CalcNormOfType(NormType, *grad_lag_x, *grad_lag_s);
2237      }
2238      curr_dual_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
2239    }
2240
2241    return result;
2242  }
2243
2244  Number
2245  IpoptCalculatedQuantities::trial_dual_infeasibility
2246  (ENormType NormType)
2247  {
2248    DBG_START_METH("IpoptCalculatedQuantities::trial_dual_infeasibility()",
2249                   dbg_verbosity);
2250    Number result;
2251
2252    SmartPtr<const Vector> x = ip_data_->trial()->x();
2253    SmartPtr<const Vector> s = ip_data_->trial()->s();
2254    SmartPtr<const Vector> y_c = ip_data_->trial()->y_c();
2255    SmartPtr<const Vector> y_d = ip_data_->trial()->y_d();
2256    SmartPtr<const Vector> z_L = ip_data_->trial()->z_L();
2257    SmartPtr<const Vector> z_U = ip_data_->trial()->z_U();
2258    SmartPtr<const Vector> v_L = ip_data_->trial()->v_L();
2259    SmartPtr<const Vector> v_U = ip_data_->trial()->v_U();
2260
2261    std::vector<const TaggedObject*> deps(8);
2262    deps[0] = GetRawPtr(x);
2263    deps[1] = GetRawPtr(s);
2264    deps[2] = GetRawPtr(y_c);
2265    deps[3] = GetRawPtr(y_d);
2266    deps[4] = GetRawPtr(z_L);
2267    deps[5] = GetRawPtr(z_U);
2268    deps[6] = GetRawPtr(v_L);
2269    deps[7] = GetRawPtr(v_U);
2270    std::vector<Number> sdeps(1);
2271    sdeps[0] = (Number)NormType;
2272
2273    if (!trial_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2274      if (!curr_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2275        SmartPtr<const Vector> grad_lag_x = trial_grad_lag_x();
2276        SmartPtr<const Vector> grad_lag_s = trial_grad_lag_s();
2277        DBG_PRINT_VECTOR(2,"grad_lag_x", *grad_lag_x);
2278        DBG_PRINT_VECTOR(2,"grad_lag_s", *grad_lag_s);
2279
2280        result = CalcNormOfType(NormType, *grad_lag_x, *grad_lag_s);
2281
2282      }
2283      trial_dual_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
2284    }
2285
2286    return result;
2287  }
2288
2289  Number
2290  IpoptCalculatedQuantities::unscaled_curr_dual_infeasibility
2291  (ENormType NormType)
2292  {
2293    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_dual_infeasibility()",
2294                   dbg_verbosity);
2295    Number result;
2296
2297    SmartPtr<const Vector> x = ip_data_->curr()->x();
2298    SmartPtr<const Vector> s = ip_data_->curr()->s();
2299    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2300    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2301    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2302    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2303    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2304    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2305
2306    std::vector<const TaggedObject*> deps(8);
2307    deps[0] = GetRawPtr(x);
2308    deps[1] = GetRawPtr(s);
2309    deps[2] = GetRawPtr(y_c);
2310    deps[3] = GetRawPtr(y_d);
2311    deps[4] = GetRawPtr(z_L);
2312    deps[5] = GetRawPtr(z_U);
2313    deps[6] = GetRawPtr(v_L);
2314    deps[7] = GetRawPtr(v_U);
2315    std::vector<Number> sdeps(1);
2316    sdeps[0] = (Number)NormType;
2317
2318    if (!unscaled_curr_dual_infeasibility_cache_.GetCachedResult(result, deps, sdeps)) {
2319      SmartPtr<const Vector> grad_lag_x =
2320        ip_nlp_->NLP_scaling()->unapply_grad_obj_scaling(curr_grad_lag_x());
2321
2322      Number obj_unscal = ip_nlp_->NLP_scaling()->unapply_obj_scaling(1.);
2323      SmartPtr<const Vector> grad_lag_s;
2324      if (obj_unscal != 1.) {
2325        SmartPtr<Vector> tmp =
2326          ip_nlp_->NLP_scaling()->apply_vector_scaling_d_NonConst(ConstPtr(curr_grad_lag_s()));
2327        tmp->Scal(obj_unscal);
2328        grad_lag_s = ConstPtr(tmp);
2329      }
2330      else {
2331        grad_lag_s = ip_nlp_->NLP_scaling()->apply_vector_scaling_d(curr_grad_lag_s());
2332      }
2333
2334      result = CalcNormOfType(NormType, *grad_lag_x, *grad_lag_s);
2335      unscaled_curr_dual_infeasibility_cache_.AddCachedResult(result, deps, sdeps);
2336    }
2337
2338    return result;
2339  }
2340
2341  Number
2342  IpoptCalculatedQuantities::curr_complementarity
2343  (Number mu, ENormType NormType)
2344  {
2345    DBG_START_METH("IpoptCalculatedQuantities::curr_complementarity()",
2346                   dbg_verbosity);
2347    Number result;
2348
2349    SmartPtr<const Vector> x = ip_data_->curr()->x();
2350    SmartPtr<const Vector> s = ip_data_->curr()->s();
2351    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2352    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2353    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2354    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2355
2356    std::vector<const TaggedObject*> deps(6);
2357    deps[0] = GetRawPtr(x);
2358    deps[1] = GetRawPtr(s);
2359    deps[2] = GetRawPtr(z_L);
2360    deps[3] = GetRawPtr(z_U);
2361    deps[4] = GetRawPtr(v_L);
2362    deps[5] = GetRawPtr(v_U);
2363    std::vector<Number> sdeps(2);
2364    sdeps[0] = (Number)NormType;
2365    sdeps[1] = mu;
2366
2367    if (!curr_complementarity_cache_.GetCachedResult(result, deps, sdeps)) {
2368      if (!trial_complementarity_cache_.GetCachedResult(result, deps, sdeps)) {
2369
2370        std::vector<SmartPtr<const Vector> > vecs(4);
2371        SmartPtr<const Vector> compl_x_L = curr_compl_x_L();
2372        SmartPtr<const Vector> compl_x_U = curr_compl_x_U();
2373        SmartPtr<const Vector> compl_s_L = curr_compl_s_L();
2374        SmartPtr<const Vector> compl_s_U = curr_compl_s_U();
2375
2376        if (mu==.0) {
2377          vecs[0] = GetRawPtr(compl_x_L);
2378          vecs[1] = GetRawPtr(compl_x_U);
2379          vecs[2] = GetRawPtr(compl_s_L);
2380          vecs[3] = GetRawPtr(compl_s_U);
2381        }
2382        else {
2383          SmartPtr<Vector> tmp = compl_x_L->MakeNew();
2384          tmp->Copy(*compl_x_L);
2385          tmp->AddScalar(-mu);
2386          vecs[0] = GetRawPtr(tmp);
2387          tmp = compl_x_U->MakeNew();
2388          tmp->Copy(*compl_x_U);
2389          tmp->AddScalar(-mu);
2390          vecs[1] = GetRawPtr(tmp);
2391          tmp = compl_s_L->MakeNew();
2392          tmp->Copy(*compl_s_L);
2393          tmp->AddScalar(-mu);
2394          vecs[2] = GetRawPtr(tmp);
2395          tmp = compl_s_U->MakeNew();
2396          tmp->Copy(*compl_s_U);
2397          tmp->AddScalar(-mu);
2398          vecs[3] = GetRawPtr(tmp);
2399        }
2400
2401        result = CalcNormOfType(NormType, vecs);
2402      }
2403
2404      curr_complementarity_cache_.AddCachedResult(result, deps, sdeps);
2405    }
2406
2407    return result;
2408  }
2409
2410  Number
2411  IpoptCalculatedQuantities::trial_complementarity
2412  (Number mu, ENormType NormType)
2413  {
2414    DBG_START_METH("IpoptCalculatedQuantities::trial_complementarity()",
2415                   dbg_verbosity);
2416    Number result;
2417
2418    SmartPtr<const Vector> x = ip_data_->trial()->x();
2419    SmartPtr<const Vector> s = ip_data_->trial()->s();
2420    SmartPtr<const Vector> z_L = ip_data_->trial()->z_L();
2421    SmartPtr<const Vector> z_U = ip_data_->trial()->z_U();
2422    SmartPtr<const Vector> v_L = ip_data_->trial()->v_L();
2423    SmartPtr<const Vector> v_U = ip_data_->trial()->v_U();
2424
2425    std::vector<const TaggedObject*> deps(6);
2426    deps[0] = GetRawPtr(x);
2427    deps[1] = GetRawPtr(s);
2428    deps[2] = GetRawPtr(z_L);
2429    deps[3] = GetRawPtr(z_U);
2430    deps[4] = GetRawPtr(v_L);
2431    deps[5] = GetRawPtr(v_U);
2432    std::vector<Number> sdeps(2);
2433    sdeps[0] = (Number)NormType;
2434    sdeps[1] = mu;
2435
2436    if (!trial_complementarity_cache_.GetCachedResult(result, deps, sdeps)) {
2437      if (!curr_complementarity_cache_.GetCachedResult(result, deps, sdeps)) {
2438
2439        std::vector<SmartPtr<const Vector> > vecs(4);
2440        SmartPtr<const Vector> compl_x_L = trial_compl_x_L();
2441        SmartPtr<const Vector> compl_x_U = trial_compl_x_U();
2442        SmartPtr<const Vector> compl_s_L = trial_compl_s_L();
2443        SmartPtr<const Vector> compl_s_U = trial_compl_s_U();
2444
2445        if (mu==.0) {
2446          vecs[0] = GetRawPtr(compl_x_L);
2447          vecs[1] = GetRawPtr(compl_x_U);
2448          vecs[2] = GetRawPtr(compl_s_L);
2449          vecs[3] = GetRawPtr(compl_s_U);
2450        }
2451        else {
2452          SmartPtr<Vector> tmp = compl_x_L->MakeNew();
2453          tmp->Copy(*compl_x_L);
2454          tmp->AddScalar(-mu);
2455          vecs[0] = GetRawPtr(tmp);
2456          tmp = compl_x_U->MakeNew();
2457          tmp->Copy(*compl_x_U);
2458          tmp->AddScalar(-mu);
2459          vecs[1] = GetRawPtr(tmp);
2460          tmp = compl_s_L->MakeNew();
2461          tmp->Copy(*compl_s_L);
2462          tmp->AddScalar(-mu);
2463          vecs[2] = GetRawPtr(tmp);
2464          tmp = compl_s_U->MakeNew();
2465          tmp->Copy(*compl_s_U);
2466          tmp->AddScalar(-mu);
2467          vecs[3] = GetRawPtr(tmp);
2468        }
2469
2470        result = CalcNormOfType(NormType, vecs);
2471      }
2472
2473      trial_complementarity_cache_.AddCachedResult(result, deps, sdeps);
2474    }
2475
2476    return result;
2477  }
2478
2479  Number
2480  IpoptCalculatedQuantities::unscaled_curr_complementarity
2481  (Number mu, ENormType NormType)
2482  {
2483    return ip_nlp_->NLP_scaling()->unapply_obj_scaling(curr_complementarity(mu, NormType));
2484  }
2485
2486  Number
2487  IpoptCalculatedQuantities::CalcCentralityMeasure(const Vector& compl_x_L,
2488      const Vector& compl_x_U,
2489      const Vector& compl_s_L,
2490      const Vector& compl_s_U)
2491  {
2492    DBG_START_METH("IpoptCalculatedQuantities::CalcCentralityMeasure()",
2493                   dbg_verbosity);
2494
2495    Number MinCompl = std::numeric_limits<Number>::max();
2496    bool have_bounds = false;
2497
2498    Index n_compl_x_L = compl_x_L.Dim();
2499    Index n_compl_x_U = compl_x_U.Dim();
2500    Index n_compl_s_L = compl_s_L.Dim();
2501    Index n_compl_s_U = compl_s_U.Dim();
2502
2503    // Compute the Minimum of all complementarities
2504    if( n_compl_x_L>0 ) {
2505      if( have_bounds ) {
2506        MinCompl = Min(MinCompl, compl_x_L.Min());
2507      }
2508      else {
2509        MinCompl = compl_x_L.Min();
2510      }
2511      have_bounds = true;
2512    }
2513    if( n_compl_x_U>0 ) {
2514      if( have_bounds ) {
2515        MinCompl = Min(MinCompl, compl_x_U.Min());
2516      }
2517      else {
2518        MinCompl = compl_x_U.Min();
2519      }
2520      have_bounds = true;
2521    }
2522    if( n_compl_s_L>0 ) {
2523      if( have_bounds ) {
2524        MinCompl = Min(MinCompl, compl_s_L.Min());
2525      }
2526      else {
2527        MinCompl = compl_s_L.Min();
2528      }
2529      have_bounds = true;
2530    }
2531    if( n_compl_s_U>0 ) {
2532      if( have_bounds ) {
2533        MinCompl = Min(MinCompl, compl_s_U.Min());
2534      }
2535      else {
2536        MinCompl = compl_s_U.Min();
2537      }
2538      have_bounds = true;
2539    }
2540
2541    // If there are no bounds, just return 0.;
2542    if (!have_bounds) {
2543      return 0.;
2544    }
2545
2546    DBG_PRINT_VECTOR(2, "compl_x_L", compl_x_L);
2547    DBG_PRINT_VECTOR(2, "compl_x_U", compl_x_U);
2548    DBG_PRINT_VECTOR(2, "compl_s_L", compl_s_L);
2549    DBG_PRINT_VECTOR(2, "compl_s_U", compl_s_U);
2550
2551    DBG_ASSERT(MinCompl>0. && "There is a zero complementarity entry");
2552
2553    Number avrg_compl = (compl_x_L.Asum() + compl_x_U.Asum() +
2554                         compl_s_L.Asum() + compl_s_U.Asum());
2555    DBG_PRINT((1,"sum_compl = %25.16e\n", avrg_compl));
2556    avrg_compl /= (n_compl_x_L + n_compl_x_U + n_compl_s_L + n_compl_s_U);
2557    DBG_PRINT((1,"avrg_compl = %25.16e\n", avrg_compl));
2558    DBG_PRINT((1,"MinCompl = %25.16e\n", MinCompl));
2559
2560    Number xi = MinCompl/avrg_compl;
2561    // The folloking line added for the case that avrg_compl is
2562    // slighly smaller than MinCompl, due to numerical roundoff
2563    xi = Min(1., xi);
2564
2565    return xi;
2566  }
2567
2568  Number
2569  IpoptCalculatedQuantities::curr_centrality_measure()
2570  {
2571    Number result;
2572
2573    SmartPtr<const Vector> x = ip_data_->curr()->x();
2574    SmartPtr<const Vector> s = ip_data_->curr()->s();
2575    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2576    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2577    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2578    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2579
2580    std::vector<const TaggedObject*> tdeps(6);
2581    tdeps[0] = GetRawPtr(x);
2582    tdeps[1] = GetRawPtr(s);
2583    tdeps[2] = GetRawPtr(z_L);
2584    tdeps[3] = GetRawPtr(z_U);
2585    tdeps[4] = GetRawPtr(v_L);
2586    tdeps[5] = GetRawPtr(v_U);
2587
2588    if (!curr_centrality_measure_cache_.GetCachedResult(result, tdeps)) {
2589      SmartPtr<const Vector> compl_x_L = curr_compl_x_L();
2590      SmartPtr<const Vector> compl_x_U = curr_compl_x_U();
2591      SmartPtr<const Vector> compl_s_L = curr_compl_s_L();
2592      SmartPtr<const Vector> compl_s_U = curr_compl_s_U();
2593
2594      result = CalcCentralityMeasure(*compl_x_L, *compl_x_U,
2595                                     *compl_s_L, *compl_s_U);
2596
2597      curr_centrality_measure_cache_.AddCachedResult(result, tdeps);
2598    }
2599    return result;
2600  }
2601
2602  Number
2603  IpoptCalculatedQuantities::curr_nlp_error()
2604  {
2605    DBG_START_METH("IpoptCalculatedQuantities::curr_nlp_error()",
2606                   dbg_verbosity);
2607    DBG_ASSERT(initialize_called_);
2608    Number result;
2609
2610    SmartPtr<const Vector> x = ip_data_->curr()->x();
2611    SmartPtr<const Vector> s = ip_data_->curr()->s();
2612    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2613    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2614    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2615    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2616    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2617    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2618
2619    std::vector<const TaggedObject*> tdeps(8);
2620    tdeps[0] = GetRawPtr(x);
2621    tdeps[1] = GetRawPtr(s);
2622    tdeps[2] = GetRawPtr(y_c);
2623    tdeps[3] = GetRawPtr(y_d);
2624    tdeps[4] = GetRawPtr(z_L);
2625    tdeps[5] = GetRawPtr(z_U);
2626    tdeps[6] = GetRawPtr(v_L);
2627    tdeps[7] = GetRawPtr(v_U);
2628
2629    if (!curr_nlp_error_cache_.GetCachedResult(result, tdeps)) {
2630      Number s_d = 0;
2631      Number s_c = 0;
2632      ComputeOptimalityErrorScaling(*ip_data_->curr()->y_c(), *ip_data_->curr()->y_d(),
2633                                    *ip_data_->curr()->z_L(), *ip_data_->curr()->z_U(),
2634                                    *ip_data_->curr()->v_L(), *ip_data_->curr()->v_U(),
2635                                    s_max_,
2636                                    s_d, s_c);
2637      DBG_PRINT((1, "s_d = %lf, s_c = %lf\n", s_d, s_c));
2638
2639      // Dual infeasibility
2640      DBG_PRINT((1, "curr_dual_infeasibility(NORM_MAX) = %8.2e\n",
2641                 curr_dual_infeasibility(NORM_MAX)));
2642      result = curr_dual_infeasibility(NORM_MAX)/s_d;
2643      /*
2644      // Primal infeasibility
2645      DBG_PRINT((1, "curr_primal_infeasibility(NORM_MAX) = %8.2e\n",
2646                 curr_primal_infeasibility(NORM_MAX)));
2647      result = Max(result, curr_primal_infeasibility(NORM_MAX));
2648      */
2649      result = Max(result, curr_nlp_constraint_violation(NORM_MAX));
2650      // Complementarity
2651      DBG_PRINT((1, "curr_complementarity(0., NORM_MAX) = %8.2e\n",
2652                 curr_complementarity(0., NORM_MAX)));
2653      result = Max(result, curr_complementarity(0., NORM_MAX)/s_c);
2654
2655      curr_nlp_error_cache_.AddCachedResult(result, tdeps);
2656    }
2657
2658    return result;
2659  }
2660
2661  Number
2662  IpoptCalculatedQuantities::unscaled_curr_nlp_error()
2663  {
2664    DBG_START_METH("IpoptCalculatedQuantities::unscaled_curr_nlp_error()",
2665                   dbg_verbosity);
2666    DBG_ASSERT(initialize_called_);
2667    Number result;
2668
2669    SmartPtr<const Vector> x = ip_data_->curr()->x();
2670    SmartPtr<const Vector> s = ip_data_->curr()->s();
2671    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2672    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2673    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2674    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2675    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2676    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2677
2678    std::vector<const TaggedObject*> tdeps(8);
2679    tdeps[0] = GetRawPtr(x);
2680    tdeps[1] = GetRawPtr(s);
2681    tdeps[2] = GetRawPtr(y_c);
2682    tdeps[3] = GetRawPtr(y_d);
2683    tdeps[4] = GetRawPtr(z_L);
2684    tdeps[5] = GetRawPtr(z_U);
2685    tdeps[6] = GetRawPtr(v_L);
2686    tdeps[7] = GetRawPtr(v_U);
2687
2688    if (!unscaled_curr_nlp_error_cache_.GetCachedResult(result, tdeps)) {
2689
2690      // Dual infeasibility
2691      result = unscaled_curr_dual_infeasibility(NORM_MAX);
2692      // Constraint violation
2693      result = Max(result, unscaled_curr_nlp_constraint_violation(NORM_MAX));
2694      // Complementarity (ToDo use unscaled?)
2695      DBG_PRINT((1, "curr_complementarity(0., NORM_MAX) = %8.2e\n",
2696                 curr_complementarity(0., NORM_MAX)));
2697      result = Max(result, unscaled_curr_complementarity(0., NORM_MAX));
2698
2699      unscaled_curr_nlp_error_cache_.AddCachedResult(result, tdeps);
2700    }
2701
2702    return result;
2703  }
2704
2705  Number
2706  IpoptCalculatedQuantities::curr_barrier_error()
2707  {
2708    DBG_START_METH("IpoptCalculatedQuantities::curr_barrier_error()",
2709                   dbg_verbosity);
2710    DBG_ASSERT(initialize_called_);
2711    Number result;
2712
2713    SmartPtr<const Vector> x = ip_data_->curr()->x();
2714    SmartPtr<const Vector> s = ip_data_->curr()->s();
2715    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2716    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2717    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2718    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2719    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2720    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2721    Number mu = ip_data_->curr_mu();
2722
2723    std::vector<const TaggedObject*> tdeps(8);
2724    tdeps[0] = GetRawPtr(x);
2725    tdeps[1] = GetRawPtr(s);
2726    tdeps[2] = GetRawPtr(y_c);
2727    tdeps[3] = GetRawPtr(y_d);
2728    tdeps[4] = GetRawPtr(z_L);
2729    tdeps[5] = GetRawPtr(z_U);
2730    tdeps[6] = GetRawPtr(v_L);
2731    tdeps[7] = GetRawPtr(v_U);
2732    std::vector<Number> sdeps(1);
2733    sdeps[0] = mu;
2734
2735    if (!curr_barrier_error_cache_.GetCachedResult(result, tdeps, sdeps)) {
2736      Number s_d = 0;
2737      Number s_c = 0;
2738      ComputeOptimalityErrorScaling(*ip_data_->curr()->y_c(), *ip_data_->curr()->y_d(),
2739                                    *ip_data_->curr()->z_L(), *ip_data_->curr()->z_U(),
2740                                    *ip_data_->curr()->v_L(), *ip_data_->curr()->v_U(),
2741                                    s_max_,
2742                                    s_d, s_c);
2743      DBG_PRINT((1, "s_d = %lf, s_c = %lf\n", s_d, s_c));
2744
2745      // Primal infeasibility
2746      result = curr_dual_infeasibility(NORM_MAX)/s_d;
2747      // Dual infeasibility
2748      result = Max(result, curr_primal_infeasibility(NORM_MAX));
2749      // Complementarity
2750      result = Max(result, curr_complementarity(mu, NORM_MAX)/s_c);
2751
2752      curr_barrier_error_cache_.AddCachedResult(result, tdeps);
2753    }
2754
2755    return result;
2756  }
2757
2758  Number
2759  IpoptCalculatedQuantities::curr_primal_dual_system_error(Number mu)
2760  {
2761    DBG_START_METH("IpoptCalculatedQuantities::curr_primal_dual_system_error()",
2762                   dbg_verbosity);
2763    DBG_ASSERT(initialize_called_);
2764    Number result;
2765
2766    SmartPtr<const Vector> x = ip_data_->curr()->x();
2767    SmartPtr<const Vector> s = ip_data_->curr()->s();
2768    SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
2769    SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
2770    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2771    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2772    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2773    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2774
2775    std::vector<const TaggedObject*> tdeps(8);
2776    tdeps[0] = GetRawPtr(x);
2777    tdeps[1] = GetRawPtr(s);
2778    tdeps[2] = GetRawPtr(y_c);
2779    tdeps[3] = GetRawPtr(y_d);
2780    tdeps[4] = GetRawPtr(z_L);
2781    tdeps[5] = GetRawPtr(z_U);
2782    tdeps[6] = GetRawPtr(v_L);
2783    tdeps[7] = GetRawPtr(v_U);
2784    std::vector<Number> sdeps(1);
2785    sdeps[0] = mu;
2786
2787    if (!curr_primal_dual_system_error_cache_.GetCachedResult(result, tdeps, sdeps)) {
2788      if (!trial_primal_dual_system_error_cache_.GetCachedResult(result, tdeps, sdeps)) {
2789        // For now we use the 1 norm, and scale each component by the number of entries...
2790        Index n_dual = x->Dim() + s->Dim();
2791        Number dual_inf = curr_dual_infeasibility(NORM_1)/((Number)n_dual);
2792
2793        Index n_primal = y_c->Dim() + y_d->Dim();
2794        Number primal_inf = 0.;
2795        if (n_primal>0) {
2796          primal_inf = curr_primal_infeasibility(NORM_1)/((Number)n_primal);
2797        }
2798
2799        Index n_cmpl = z_L->Dim() + z_U->Dim() + v_L->Dim() + v_U->Dim();
2800        Number cmpl = 0.;
2801        if (n_cmpl>0) {
2802          cmpl = curr_complementarity(mu, NORM_1)/((Number)n_cmpl);
2803        }
2804
2805        result = dual_inf + primal_inf + cmpl;
2806      }
2807      curr_primal_dual_system_error_cache_.AddCachedResult(result, tdeps, sdeps);
2808    }
2809
2810    return result;
2811  }
2812
2813  Number
2814  IpoptCalculatedQuantities::trial_primal_dual_system_error(Number mu)
2815  {
2816    DBG_START_METH("IpoptCalculatedQuantities::trial_primal_dual_system_error()",
2817                   dbg_verbosity);
2818    DBG_ASSERT(initialize_called_);
2819    Number result;
2820
2821    SmartPtr<const Vector> x = ip_data_->trial()->x();
2822    SmartPtr<const Vector> s = ip_data_->trial()->s();
2823    SmartPtr<const Vector> y_c = ip_data_->trial()->y_c();
2824    SmartPtr<const Vector> y_d = ip_data_->trial()->y_d();
2825    SmartPtr<const Vector> z_L = ip_data_->trial()->z_L();
2826    SmartPtr<const Vector> z_U = ip_data_->trial()->z_U();
2827    SmartPtr<const Vector> v_L = ip_data_->trial()->v_L();
2828    SmartPtr<const Vector> v_U = ip_data_->trial()->v_U();
2829
2830    std::vector<const TaggedObject*> tdeps(8);
2831    tdeps[0] = GetRawPtr(x);
2832    tdeps[1] = GetRawPtr(s);
2833    tdeps[2] = GetRawPtr(y_c);
2834    tdeps[3] = GetRawPtr(y_d);
2835    tdeps[4] = GetRawPtr(z_L);
2836    tdeps[5] = GetRawPtr(z_U);
2837    tdeps[6] = GetRawPtr(v_L);
2838    tdeps[7] = GetRawPtr(v_U);
2839    std::vector<Number> sdeps(1);
2840    sdeps[0] = mu;
2841
2842    if (!trial_primal_dual_system_error_cache_.GetCachedResult(result, tdeps, sdeps)) {
2843      if (!curr_primal_dual_system_error_cache_.GetCachedResult(result, tdeps, sdeps)) {
2844        // For now we use the 1 norm, and scale each component by the number of entries...
2845        Index n_dual = x->Dim() + s->Dim();
2846        Number dual_inf = trial_dual_infeasibility(NORM_1)/((Number)n_dual);
2847
2848        Index n_primal = y_c->Dim() + y_d->Dim();
2849        Number primal_inf = 0.;
2850        if (n_primal>0) {
2851          primal_inf = trial_primal_infeasibility(NORM_1)/((Number)n_primal);
2852        }
2853
2854        Index n_cmpl = z_L->Dim() + z_U->Dim() + v_L->Dim() + v_U->Dim();
2855        Number cmpl = 0.;
2856        if (n_cmpl>0) {
2857          cmpl = trial_complementarity(mu, NORM_1)/((Number)n_cmpl);
2858        }
2859
2860        result = dual_inf + primal_inf + cmpl;
2861      }
2862      trial_primal_dual_system_error_cache_.AddCachedResult(result, tdeps, sdeps);
2863    }
2864
2865    return result;
2866  }
2867
2868  ///////////////////////////////////////////////////////////////////////////
2869  //                Fraction-to-the-boundary step sizes                    //
2870  ///////////////////////////////////////////////////////////////////////////
2871
2872  Number
2873  IpoptCalculatedQuantities::CalcFracToBound(const Vector& slack_L,
2874      Vector& tmp_L,
2875      const Matrix& P_L,
2876      const Vector& slack_U,
2877      Vector& tmp_U,
2878      const Matrix& P_U,
2879      const Vector& delta,
2880      Number tau)
2881  {
2882    DBG_START_METH("IpoptCalculatedQuantities::CalcFracToBound",
2883                   dbg_verbosity);
2884
2885    Number alpha_L = 1.0;
2886    Number alpha_U = 1.0;
2887    if (slack_L.Dim() > 0) {
2888      P_L.TransMultVector(1.0, delta, 0.0, tmp_L);
2889      alpha_L = slack_L.FracToBound(tmp_L, tau);
2890    }
2891
2892    if (slack_U.Dim() > 0) {
2893      P_U.TransMultVector(-1.0, delta, 0.0, tmp_U);
2894      alpha_U = slack_U.FracToBound(tmp_U, tau);
2895    }
2896
2897    DBG_PRINT((1,"alpha_L = %lf, alpha_U = %lf\n", alpha_L, alpha_U));
2898    DBG_ASSERT(alpha_L >= 0.0 && alpha_L <= 1.0
2899               && alpha_U >=0.0 && alpha_U <= 1.0);
2900
2901    return Min(alpha_L, alpha_U);
2902  }
2903
2904  Number
2905  IpoptCalculatedQuantities::primal_frac_to_the_bound(Number tau,
2906      const Vector& delta_x,
2907      const Vector& delta_s)
2908  {
2909    DBG_START_METH("IpoptCalculatedQuantities::primal_frac_to_the_bound",
2910                   dbg_verbosity);
2911    Number result;
2912    SmartPtr<const Vector> x = ip_data_->curr()->x();
2913    SmartPtr<const Vector> s = ip_data_->curr()->s();
2914    std::vector<const TaggedObject*> tdeps(4);
2915    tdeps[0] = GetRawPtr(x);
2916    tdeps[1] = GetRawPtr(s);
2917    tdeps[2] = &delta_x;
2918    tdeps[3] = &delta_s;
2919
2920    std::vector<Number> sdeps(1);
2921    sdeps[0] = tau;
2922
2923    if (!primal_frac_to_the_bound_cache_.GetCachedResult(result, tdeps, sdeps)) {
2924      result = Min(CalcFracToBound(*curr_slack_x_L(), Tmp_x_L(), *ip_nlp_->Px_L(),
2925                                   *curr_slack_x_U(), Tmp_x_U(), *ip_nlp_->Px_U(),
2926                                   delta_x, tau),
2927                   CalcFracToBound(*curr_slack_s_L(), Tmp_s_L(), *ip_nlp_->Pd_L(),
2928                                   *curr_slack_s_U(), Tmp_s_U(), *ip_nlp_->Pd_U(),
2929                                   delta_s, tau));
2930
2931      primal_frac_to_the_bound_cache_.AddCachedResult(result, tdeps, sdeps);
2932    }
2933
2934    return result;
2935  }
2936
2937  Number
2938  IpoptCalculatedQuantities::curr_primal_frac_to_the_bound(Number tau)
2939  {
2940    DBG_START_METH("IpoptCalculatedQuantities::curr_primal_frac_to_the_bound()",
2941                   dbg_verbosity);
2942    return primal_frac_to_the_bound(tau, *ip_data_->delta()->x(),
2943                                    *ip_data_->delta()->s());
2944  }
2945
2946  Number
2947  IpoptCalculatedQuantities::uncached_dual_frac_to_the_bound(
2948    Number tau,
2949    const Vector& delta_z_L,
2950    const Vector& delta_z_U,
2951    const Vector& delta_v_L,
2952    const Vector& delta_v_U)
2953  {
2954    DBG_START_METH("IpoptCalculatedQuantities::uncached_dual_frac_to_the_bound",
2955                   dbg_verbosity);
2956    Number result;
2957
2958    result = ip_data_->curr()->z_L()->FracToBound(delta_z_L, tau);
2959    result = Min(result, ip_data_->curr()->z_U()->FracToBound(delta_z_U, tau));
2960    result = Min(result, ip_data_->curr()->v_L()->FracToBound(delta_v_L, tau));
2961    result = Min(result, ip_data_->curr()->v_U()->FracToBound(delta_v_U, tau));
2962
2963    return result;
2964  }
2965
2966  Number
2967  IpoptCalculatedQuantities::dual_frac_to_the_bound(
2968    Number tau,
2969    const Vector& delta_z_L,
2970    const Vector& delta_z_U,
2971    const Vector& delta_v_L,
2972    const Vector& delta_v_U)
2973  {
2974    DBG_START_METH("IpoptCalculatedQuantities::dual_frac_to_the_bound",
2975                   dbg_verbosity);
2976    Number result;
2977
2978    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
2979    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
2980    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
2981    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
2982    std::vector<const TaggedObject*> tdeps(8);
2983    tdeps[0] = GetRawPtr(z_L);
2984    tdeps[1] = GetRawPtr(z_U);
2985    tdeps[2] = GetRawPtr(v_L);
2986    tdeps[3] = GetRawPtr(v_U);
2987    tdeps[4] = &delta_z_L;
2988    tdeps[5] = &delta_z_U;
2989    tdeps[6] = &delta_v_L;
2990    tdeps[7] = &delta_v_U;
2991
2992    std::vector<Number> sdeps(1);
2993    sdeps[0] = tau;
2994
2995    if (!dual_frac_to_the_bound_cache_.GetCachedResult(result, tdeps, sdeps)) {
2996      result = z_L->FracToBound(delta_z_L, tau);
2997      result = Min(result, z_U->FracToBound(delta_z_U, tau));
2998      result = Min(result, v_L->FracToBound(delta_v_L, tau));
2999      result = Min(result, v_U->FracToBound(delta_v_U, tau));
3000
3001      dual_frac_to_the_bound_cache_.AddCachedResult(result, tdeps, sdeps);
3002    }
3003
3004    return result;
3005  }
3006
3007  Number
3008  IpoptCalculatedQuantities::curr_dual_frac_to_the_bound(Number tau)
3009  {
3010    DBG_START_METH("IpoptCalculatedQuantities::curr_dual_frac_to_the_bound()",
3011                   dbg_verbosity);
3012    return dual_frac_to_the_bound(tau, *ip_data_->delta()->z_L(),
3013                                  *ip_data_->delta()->z_U(),
3014                                  *ip_data_->delta()->v_L(),
3015                                  *ip_data_->delta()->v_U());
3016  }
3017
3018  Number
3019  IpoptCalculatedQuantities::uncached_slack_frac_to_the_bound(
3020    Number tau,
3021    const Vector& delta_x_L,
3022    const Vector& delta_x_U,
3023    const Vector& delta_s_L,
3024    const Vector& delta_s_U)
3025  {
3026    DBG_START_METH("IpoptCalculatedQuantities::slack_frac_to_the_bound",
3027                   dbg_verbosity);
3028    Number result;
3029
3030    SmartPtr<const Vector> x_L = curr_slack_x_L();
3031    SmartPtr<const Vector> x_U = curr_slack_x_U();
3032    SmartPtr<const Vector> s_L = curr_slack_s_L();
3033    SmartPtr<const Vector> s_U = curr_slack_s_U();
3034
3035    result = x_L->FracToBound(delta_x_L, tau);
3036    result = Min(result, x_U->FracToBound(delta_x_U, tau));
3037    result = Min(result, s_L->FracToBound(delta_s_L, tau));
3038    result = Min(result, s_U->FracToBound(delta_s_U, tau));
3039
3040    return result;
3041  }
3042
3043  ///////////////////////////////////////////////////////////////////////////
3044  //                             Sigma Matrices                            //
3045  ///////////////////////////////////////////////////////////////////////////
3046
3047  SmartPtr<const Vector>
3048  IpoptCalculatedQuantities::curr_sigma_x()
3049  {
3050    DBG_START_METH("IpoptCalculatedQuantities::curr_sigma_x()",
3051                   dbg_verbosity);
3052    SmartPtr<const Vector> result;
3053    SmartPtr<const Vector> x = ip_data_->curr()->x();
3054    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
3055    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
3056
3057    if (!curr_sigma_x_cache_.GetCachedResult3Dep(result, *x, *z_L, *z_U)) {
3058      SmartPtr<Vector> sigma = x->MakeNew();
3059
3060      sigma->Set(0.);
3061      ip_nlp_->Px_L()->AddMSinvZ(1., *curr_slack_x_L(), *z_L, *sigma);
3062      ip_nlp_->Px_U()->AddMSinvZ(1., *curr_slack_x_U(), *z_U, *sigma);
3063
3064      DBG_PRINT_VECTOR(2,"sigma_x", *sigma);
3065
3066      result = ConstPtr(sigma);
3067      curr_sigma_x_cache_.AddCachedResult3Dep(result, *x, *z_L, *z_U);
3068    }
3069
3070    return result;
3071  }
3072
3073  SmartPtr<const Vector>
3074  IpoptCalculatedQuantities::curr_sigma_s()
3075  {
3076    DBG_START_METH("IpoptCalculatedQuantities::curr_sigma_s()",
3077                   dbg_verbosity);
3078    SmartPtr<const Vector> result;
3079    SmartPtr<const Vector> s = ip_data_->curr()->s();
3080    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
3081    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
3082
3083    if (!curr_sigma_s_cache_.GetCachedResult3Dep(result, *s, *v_L, *v_U)) {
3084      SmartPtr<Vector> sigma = s->MakeNew();
3085
3086      sigma->Set(0.);
3087      ip_nlp_->Pd_L()->AddMSinvZ(1., *curr_slack_s_L(), *v_L, *sigma);
3088      ip_nlp_->Pd_U()->AddMSinvZ(1., *curr_slack_s_U(), *v_U, *sigma);
3089
3090      result = ConstPtr(sigma);
3091      curr_sigma_s_cache_.AddCachedResult3Dep(result, *s, *v_L, *v_U);
3092    }
3093
3094    return result;
3095  }
3096
3097  Number
3098  IpoptCalculatedQuantities::curr_avrg_compl()
3099  {
3100    DBG_START_METH("IpoptCalculatedQuantities::curr_avrg_compl()",
3101                   dbg_verbosity);
3102
3103    Number result;
3104
3105    SmartPtr<const Vector> x = ip_data_->curr()->x();
3106    SmartPtr<const Vector> s = ip_data_->curr()->s();
3107    SmartPtr<const Vector> z_L = ip_data_->curr()->z_L();
3108    SmartPtr<const Vector> z_U = ip_data_->curr()->z_U();
3109    SmartPtr<const Vector> v_L = ip_data_->curr()->v_L();
3110    SmartPtr<const Vector> v_U = ip_data_->curr()->v_U();
3111
3112    std::vector<const TaggedObject*> tdeps(6);
3113    tdeps[0] = GetRawPtr(x);
3114    tdeps[1] = GetRawPtr(s);
3115    tdeps[2] = GetRawPtr(z_L);
3116    tdeps[3] = GetRawPtr(z_U);
3117    tdeps[4] = GetRawPtr(v_L);
3118    tdeps[5] = GetRawPtr(v_U);
3119
3120    if (!curr_avrg_compl_cache_.GetCachedResult(result, tdeps)) {
3121      if (!trial_avrg_compl_cache_.GetCachedResult(result, tdeps)) {
3122
3123        SmartPtr<const Vector> slack_x_L = curr_slack_x_L();
3124        SmartPtr<const Vector> slack_x_U = curr_slack_x_U();
3125        SmartPtr<const Vector> slack_s_L = curr_slack_s_L();
3126        SmartPtr<const Vector> slack_s_U = curr_slack_s_U();
3127
3128        Index ncomps = z_L->Dim() + z_U->Dim() + v_L->Dim() + v_U->Dim();
3129
3130        if (ncomps>0) {
3131          result = z_L->Dot(*slack_x_L);
3132          result += z_U->Dot(*slack_x_U);
3133          result += v_L->Dot(*slack_s_L);
3134          result += v_U->Dot(*slack_s_U);
3135
3136          result /= (Number)ncomps;
3137        }
3138        else {
3139          result = 0.;
3140        }
3141      }
3142
3143      curr_avrg_compl_cache_.AddCachedResult(result, tdeps);
3144    }
3145
3146    return result;
3147  }
3148
3149  Number
3150  IpoptCalculatedQuantities::trial_avrg_compl()
3151  {
3152    DBG_START_METH("IpoptCalculatedQuantities::trial_avrg_compl()",
3153                   dbg_verbosity);
3154
3155    Number result;
3156
3157    SmartPtr<const Vector> x = ip_data_->trial()->x();
3158    SmartPtr<const Vector> s = ip_data_->trial()->s();
3159    SmartPtr<const Vector> z_L = ip_data_->trial()->z_L();
3160    SmartPtr<const Vector> z_U = ip_data_->trial()->z_U();
3161    SmartPtr<const Vector> v_L = ip_data_->trial()->v_L();
3162    SmartPtr<const Vector> v_U = ip_data_->trial()->v_U();
3163
3164    std::vector<const TaggedObject*> tdeps(6);
3165    tdeps[0] = GetRawPtr(x);
3166    tdeps[1] = GetRawPtr(s);
3167    tdeps[2] = GetRawPtr(z_L);
3168    tdeps[3] = GetRawPtr(z_U);
3169    tdeps[4] = GetRawPtr(v_L);
3170    tdeps[5] = GetRawPtr(v_U);
3171
3172    if (!trial_avrg_compl_cache_.GetCachedResult(result, tdeps)) {
3173      if (!curr_avrg_compl_cache_.GetCachedResult(result, tdeps)) {
3174
3175        SmartPtr<const Vector> slack_x_L = trial_slack_x_L();
3176        SmartPtr<const Vector> slack_x_U = trial_slack_x_U();
3177        SmartPtr<const Vector> slack_s_L = trial_slack_s_L();
3178        SmartPtr<const Vector> slack_s_U = trial_slack_s_U();
3179
3180        Index ncomps = z_L->Dim() + z_U->Dim() + v_L->Dim() + v_U->Dim();
3181
3182        if (ncomps>0) {
3183          result = z_L->Dot(*slack_x_L);
3184          result += z_U->Dot(*slack_x_U);
3185          result += v_L->Dot(*slack_s_L);
3186          result += v_U->Dot(*slack_s_U);
3187
3188          result /= (Number)ncomps;
3189        }
3190        else {
3191          result = 0.;
3192        }
3193      }
3194
3195      trial_avrg_compl_cache_.AddCachedResult(result, tdeps);
3196    }
3197
3198    return result;
3199  }
3200
3201  void IpoptCalculatedQuantities::ComputeOptimalityErrorScaling(const Vector& y_c, const Vector& y_d,
3202      const Vector& z_L, const Vector& z_U,
3203      const Vector& v_L, const Vector& v_U,
3204      Number s_max,
3205      Number& s_d, Number& s_c)
3206  {
3207    DBG_ASSERT(initialize_called_);
3208
3209    s_c = z_L.Asum() + z_U.Asum() + v_L.Asum() + v_U.Asum();
3210    Number n = (z_L.Dim() + z_U.Dim() + v_L.Dim() + v_U.Dim());
3211    if (n == 0) {
3212      s_c = 1.0;
3213    }
3214    else {
3215      s_c = s_c / n;
3216      s_c = Max(s_max, s_c)/s_max;
3217    }
3218
3219    s_d = y_c.Asum() + y_d.Asum() + z_L.Asum() + z_U.Asum() + v_L.Asum() + v_U.Asum();
3220    n = (y_c.Dim() + y_d.Dim() + z_L.Dim() + z_U.Dim() + v_L.Dim() + v_U.Dim());
3221    if ( n == 0 ) {
3222      s_d = 1.0;
3223    }
3224    else {
3225      s_d = s_d / n;
3226      s_d = Max(s_max, s_d)/s_max;
3227    }
3228  }
3229
3230  Number IpoptCalculatedQuantities::curr_gradBarrTDelta()
3231  {
3232    DBG_START_METH("IpoptCalculatedQuantities::curr_gradBarrTDelta()",
3233                   dbg_verbosity);
3234
3235    Number result;
3236
3237    SmartPtr<const Vector> x = ip_data_->curr()->x();
3238    SmartPtr<const Vector> s = ip_data_->curr()->s();
3239    SmartPtr<const Vector> delta_x = ip_data_->delta()->x();
3240    SmartPtr<const Vector> delta_s = ip_data_->delta()->s();
3241    std::vector<const TaggedObject*> tdeps(4);
3242    tdeps[0] = GetRawPtr(x);
3243    tdeps[1] = GetRawPtr(s);
3244    tdeps[2] = GetRawPtr(delta_x);
3245    tdeps[3] = GetRawPtr(delta_s);
3246    Number mu = ip_data_->curr_mu();
3247    std::vector<Number> sdeps(1);
3248    sdeps[0] = mu;
3249    DBG_PRINT((1,"curr_mu=%e\n",mu));
3250
3251    if (!curr_gradBarrTDelta_cache_.GetCachedResult(result, tdeps, sdeps)) {
3252      result = curr_grad_barrier_obj_x()->Dot(*delta_x) +
3253               curr_grad_barrier_obj_s()->Dot(*delta_s);
3254
3255      curr_gradBarrTDelta_cache_.AddCachedResult(result, tdeps, sdeps);
3256    }
3257    return result;
3258  }
3259
3260  bool IpoptCalculatedQuantities::IsSquareProblem() const
3261  {
3262    return (ip_data_->curr()->x()->Dim() == ip_data_->curr()->y_c()->Dim());
3263  }
3264
3265  Vector& IpoptCalculatedQuantities::Tmp_x()
3266  {
3267    if (!IsValid(tmp_x_)) {
3268      tmp_x_ = ip_data_->curr()->x()->MakeNew();
3269    }
3270    return *tmp_x_;
3271  }
3272
3273  Vector& IpoptCalculatedQuantities::Tmp_s()
3274  {
3275    if (!IsValid(tmp_s_)) {
3276      tmp_s_ = ip_data_->curr()->s()->MakeNew();
3277    }
3278    return *tmp_s_;
3279  }
3280
3281  Vector& IpoptCalculatedQuantities::Tmp_c()
3282  {
3283    if (!IsValid(tmp_c_)) {
3284      tmp_c_ = ip_data_->curr()->y_c()->MakeNew();
3285    }
3286    return *tmp_c_;
3287  }
3288
3289  Vector& IpoptCalculatedQuantities::Tmp_d()
3290  {
3291    if (!IsValid(tmp_d_)) {
3292      tmp_d_ = ip_data_->curr()->y_d()->MakeNew();
3293    }
3294    return *tmp_d_;
3295  }
3296
3297  Vector& IpoptCalculatedQuantities::Tmp_x_L()
3298  {
3299    if (!IsValid(tmp_x_L_)) {
3300      tmp_x_L_ = ip_nlp_->x_L()->MakeNew();
3301    }
3302    return *tmp_x_L_;
3303  }
3304
3305  Vector& IpoptCalculatedQuantities::Tmp_x_U()
3306  {
3307    if (!IsValid(tmp_x_U_)) {
3308      tmp_x_U_ = ip_nlp_->x_U()->MakeNew();
3309    }
3310    return *tmp_x_U_;
3311  }
3312
3313  Vector& IpoptCalculatedQuantities::Tmp_s_L()
3314  {
3315    if (!IsValid(tmp_s_L_)) {
3316      tmp_s_L_ = ip_nlp_->d_L()->MakeNew();
3317    }
3318    return *tmp_s_L_;
3319  }
3320
3321  Vector& IpoptCalculatedQuantities::Tmp_s_U()
3322  {
3323    if (!IsValid(tmp_s_U_)) {
3324      tmp_s_U_ = ip_nlp_->d_U()->MakeNew();
3325    }
3326    return *tmp_s_U_;
3327  }
3328
3329
3330} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.