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

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

adapted for Portland compilers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.6 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: IpRestoIterationOutput.cpp 567 2005-11-01 23:27:11Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter              IBM    2004-09-23
8
9#include "IpRestoIterationOutput.hpp"
10#include "IpRestoIpoptNLP.hpp"
11
12#ifdef HAVE_CMATH
13# include <cmath>
14#else
15# ifdef HAVE_MATH_H
16#  include <math.h>
17# else
18#  error "don't have header file for math"
19# endif
20#endif
21
22namespace Ipopt
23{
24  RestoIterationOutput::RestoIterationOutput(const SmartPtr<OrigIterationOutput>& resto_orig_iteration_output)
25      :
26      resto_orig_iteration_output_(resto_orig_iteration_output)
27  {}
28
29  RestoIterationOutput::~RestoIterationOutput()
30  {}
31
32  bool RestoIterationOutput::InitializeImpl(const OptionsList& options,
33      const std::string& prefix)
34  {
35    options.GetBoolValue("print_info_string", print_info_string_, prefix);
36
37    bool retval = true;
38    if (IsValid(resto_orig_iteration_output_)) {
39      retval = resto_orig_iteration_output_->Initialize(Jnlst(), IpNLP(),
40               IpData(), IpCq(),
41               options, prefix);
42    }
43    return retval;
44  }
45
46  void RestoIterationOutput::WriteOutput()
47  {
48    // Get pointers to the Original NLP objects
49    const RestoIpoptNLP* resto_ipopt_nlp =
50      dynamic_cast<const RestoIpoptNLP*>(&IpNLP());
51    DBG_ASSERT(resto_ipopt_nlp);
52
53    SmartPtr<IpoptData> orig_ip_data = &resto_ipopt_nlp->OrigIpData();
54    SmartPtr<IpoptNLP> orig_ip_nlp = &resto_ipopt_nlp->OrigIpNLP();
55    SmartPtr<IpoptCalculatedQuantities> orig_ip_cq =
56      &resto_ipopt_nlp->OrigIpCq();
57
58    // Set the iteration counter for the original NLP to the current value
59    Index iter = IpData().iter_count();
60    orig_ip_data->Set_iter_count(iter);
61
62    // If a resto_orig_iteration_output object was given, first do the
63    // WriteOutput method with that one
64    if (IsValid(resto_orig_iteration_output_)) {
65      resto_orig_iteration_output_->WriteOutput();
66    }
67
68    //////////////////////////////////////////////////////////////////////
69    //         First print the summary line for the iteration           //
70    //////////////////////////////////////////////////////////////////////
71
72    std::string header =
73      "iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls\n";
74    Jnlst().Printf(J_DETAILED, J_MAIN,
75                   "\n\n**************************************************\n");
76    Jnlst().Printf(J_DETAILED, J_MAIN,
77                   "*** Summary of Iteration %d for original NLP:", IpData().iter_count());
78    Jnlst().Printf(J_DETAILED, J_MAIN,
79                   "\n**************************************************\n\n");
80    if (iter%10 == 0 && !IsValid(resto_orig_iteration_output_)) {
81      // output the header
82      Jnlst().Printf(J_ITERSUMMARY, J_MAIN, header.c_str());
83    }
84    else {
85      Jnlst().Printf(J_DETAILED, J_MAIN, header.c_str());
86    }
87
88    // For now, just print the total NLP error for the restoration
89    // phase problem in the dual infeasibility column
90    Number inf_du =
91      IpCq().curr_dual_infeasibility(NORM_MAX);
92
93    Number mu = IpData().curr_mu();
94    Number dnrm = 0.;
95    if (IsValid(IpData().delta()) && IsValid(IpData().delta()->x()) && IsValid(IpData().delta()->s())) {
96      dnrm = Max(IpData().delta()->x()->Amax(), IpData().delta()->s()->Amax());
97    }
98
99    // Set  the trial  values  for  the original  Data  object to  the
100    // current restoration phase values
101    SmartPtr<const Vector> x = IpData().curr()->x();
102    const CompoundVector* cx =
103      dynamic_cast<const CompoundVector*>(GetRawPtr(x));
104    DBG_ASSERT(cx);
105
106    SmartPtr<IteratesVector> trial = orig_ip_data->trial()->MakeNewContainer();
107    trial->Set_x(*cx->GetComp(0));
108    trial->Set_s(*IpData().curr()->s());
109    orig_ip_data->set_trial(trial);
110
111    // Compute primal infeasibility
112    Number inf_pr = orig_ip_cq->trial_primal_infeasibility(NORM_MAX);
113    Number f = orig_ip_cq->unscaled_trial_f();
114
115    // Retrieve some information set in the different parts of the algorithm
116    char info_iter='r';
117
118    Number alpha_primal = IpData().info_alpha_primal();
119    char alpha_primal_char = IpData().info_alpha_primal_char();
120    Number alpha_dual = IpData().info_alpha_dual();
121    Number regu_x = IpData().info_regu_x();
122    char regu_x_buf[8];
123    char dashes[]="   - ";
124    char *regu_x_ptr;
125    if (regu_x==.0) {
126      regu_x_ptr = dashes;
127    }
128    else {
129      sprintf(regu_x_buf, "%5.1f", log10(regu_x));
130      regu_x_ptr = regu_x_buf;
131    }
132    Index ls_count = IpData().info_ls_count();
133    const std::string info_string = IpData().info_string();
134
135    Jnlst().Printf(J_ITERSUMMARY, J_MAIN,
136                   "%4d%c%14.7e %7.2e %7.2e %5.1f %7.2e %5s %7.2e %7.2e%c%3d",
137                   iter, info_iter, f, inf_pr, inf_du, log10(mu), dnrm, regu_x_ptr,
138                   alpha_dual, alpha_primal, alpha_primal_char,
139                   ls_count);
140    if (print_info_string_) {
141      Jnlst().Printf(J_ITERSUMMARY, J_MAIN, " %s", info_string.c_str());
142    }
143    else {
144      Jnlst().Printf(J_DETAILED, J_MAIN, " %s", info_string.c_str());
145    }
146    Jnlst().Printf(J_ITERSUMMARY, J_MAIN, "\n");
147
148    //////////////////////////////////////////////////////////////////////
149    //           Now if desired more detail on the iterates             //
150    //////////////////////////////////////////////////////////////////////
151
152    if (Jnlst().ProduceOutput(J_DETAILED, J_MAIN)) {
153      Jnlst().Printf(J_DETAILED, J_MAIN,
154                     "\n**************************************************\n");
155      Jnlst().Printf(J_DETAILED, J_MAIN,
156                     "*** Beginning Iteration %d from the following point:",
157                     IpData().iter_count());
158      Jnlst().Printf(J_DETAILED, J_MAIN,
159                     "\n**************************************************\n\n");
160
161      Jnlst().Printf(J_DETAILED, J_MAIN,
162                     "Primal infeasibility for restoration phase problem = %.16e\n",
163                     IpCq().curr_primal_infeasibility(NORM_MAX));
164      Jnlst().Printf(J_DETAILED, J_MAIN,
165                     "Dual infeasibility for restoration phase problem   = %.16e\n",
166                     IpCq().curr_dual_infeasibility(NORM_MAX));
167
168      Jnlst().Printf(J_DETAILED, J_MAIN,
169                     "||curr_x||_inf   = %.16e\n", IpData().curr()->x()->Amax());
170      Jnlst().Printf(J_DETAILED, J_MAIN,
171                     "||curr_s||_inf   = %.16e\n", IpData().curr()->s()->Amax());
172      Jnlst().Printf(J_DETAILED, J_MAIN,
173                     "||curr_y_c||_inf = %.16e\n", IpData().curr()->y_c()->Amax());
174      Jnlst().Printf(J_DETAILED, J_MAIN,
175                     "||curr_y_d||_inf = %.16e\n", IpData().curr()->y_d()->Amax());
176      Jnlst().Printf(J_DETAILED, J_MAIN,
177                     "||curr_z_L||_inf = %.16e\n", IpData().curr()->z_L()->Amax());
178      Jnlst().Printf(J_DETAILED, J_MAIN,
179                     "||curr_z_U||_inf = %.16e\n", IpData().curr()->z_U()->Amax());
180      Jnlst().Printf(J_DETAILED, J_MAIN,
181                     "||curr_v_L||_inf = %.16e\n", IpData().curr()->v_L()->Amax());
182      Jnlst().Printf(J_DETAILED, J_MAIN,
183                     "||curr_v_U||_inf = %.16e\n", IpData().curr()->v_U()->Amax());
184    }
185    if (Jnlst().ProduceOutput(J_MOREDETAILED, J_MAIN)) {
186      if (IsValid(IpData().delta())) {
187        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
188                       "\n||delta_x||_inf   = %.16e\n", IpData().delta()->x()->Amax());
189        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
190                       "||delta_s||_inf   = %.16e\n", IpData().delta()->s()->Amax());
191        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
192                       "||delta_y_c||_inf = %.16e\n", IpData().delta()->y_c()->Amax());
193        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
194                       "||delta_y_d||_inf = %.16e\n", IpData().delta()->y_d()->Amax());
195        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
196                       "||delta_z_L||_inf = %.16e\n", IpData().delta()->z_L()->Amax());
197        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
198                       "||delta_z_U||_inf = %.16e\n", IpData().delta()->z_U()->Amax());
199        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
200                       "||delta_v_L||_inf = %.16e\n", IpData().delta()->v_L()->Amax());
201        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
202                       "||delta_v_U||_inf = %.16e\n", IpData().delta()->v_U()->Amax());
203      }
204      else {
205        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
206                       "\nNo search direction has been computed yet.\n");
207      }
208    }
209    if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
210      IpData().curr()->x()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_x");
211      IpData().curr()->s()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_s");
212
213      IpData().curr()->y_c()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_y_c");
214      IpData().curr()->y_d()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_y_d");
215
216      IpCq().curr_slack_x_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_slack_x_L");
217      IpCq().curr_slack_x_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_slack_x_U");
218      IpData().curr()->z_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_z_L");
219      IpData().curr()->z_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_z_U");
220
221      IpCq().curr_slack_s_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_slack_s_L");
222      IpCq().curr_slack_s_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_slack_s_U");
223      IpData().curr()->v_L()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_v_L");
224      IpData().curr()->v_U()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_v_U");
225    }
226    if (Jnlst().ProduceOutput(J_MOREVECTOR, J_MAIN)) {
227      IpCq().curr_grad_lag_x()->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "curr_grad_lag_x");
228      IpCq().curr_grad_lag_s()->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "curr_grad_lag_s");
229      if (IsValid(IpData().delta())) {
230        IpData().delta()->Print(Jnlst(), J_MOREVECTOR, J_MAIN, "delta");
231      }
232    }
233
234    if (Jnlst().ProduceOutput(J_DETAILED, J_MAIN)) {
235      Jnlst().Printf(J_DETAILED, J_MAIN,
236                     "\n\n***Current NLP Values for Iteration (Restoration phase problem) %d:\n",
237                     IpData().iter_count());
238      Jnlst().Printf(J_DETAILED, J_MAIN, "\n                                   (scaled)                 (unscaled)\n");
239      Jnlst().Printf(J_DETAILED, J_MAIN, "Objective...............: %24.16e  %24.16e\n", IpCq().curr_f(), IpCq().unscaled_curr_f());
240      Jnlst().Printf(J_DETAILED, J_MAIN, "Dual infeasibility......: %24.16e  %24.16e\n", IpCq().curr_dual_infeasibility(NORM_MAX), IpCq().unscaled_curr_dual_infeasibility(NORM_MAX));
241      Jnlst().Printf(J_DETAILED, J_MAIN, "Constraint violation....: %24.16e  %24.16e\n", IpCq().curr_nlp_constraint_violation(NORM_MAX), IpCq().unscaled_curr_nlp_constraint_violation(NORM_MAX));
242      Jnlst().Printf(J_DETAILED, J_MAIN, "Complementarity.........: %24.16e  %24.16e\n", IpCq().curr_complementarity(0., NORM_MAX), IpCq().unscaled_curr_complementarity(0., NORM_MAX));
243      Jnlst().Printf(J_DETAILED, J_MAIN, "Overall NLP error.......: %24.16e  %24.16e\n\n", IpCq().curr_nlp_error(), IpCq().unscaled_curr_nlp_error());
244    }
245    if (Jnlst().ProduceOutput(J_VECTOR, J_MAIN)) {
246      IpCq().curr_grad_f()->Print(Jnlst(), J_VECTOR, J_MAIN, "grad_f");
247      IpCq().curr_c()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_c");
248      IpCq().curr_d()->Print(Jnlst(), J_VECTOR, J_MAIN, "curr_d");
249      IpCq().curr_d_minus_s()->Print(Jnlst(), J_VECTOR, J_MAIN,
250                                     "curr_d - curr_s");
251    }
252
253    if (Jnlst().ProduceOutput(J_MATRIX, J_MAIN)) {
254      IpCq().curr_jac_c()->Print(Jnlst(), J_MATRIX, J_MAIN, "jac_c");
255      IpCq().curr_jac_d()->Print(Jnlst(), J_MATRIX, J_MAIN, "jac_d");
256      IpCq().curr_exact_hessian()->Print(Jnlst(), J_MATRIX, J_MAIN, "h");
257    }
258
259    Jnlst().Printf(J_DETAILED, J_MAIN, "\n\n");
260  }
261
262} // namespace Ipopt
Note: See TracBrowser for help on using the repository browser.