1 // Copyright (C) 2008 International Business Machines and others.
2 // All Rights Reserved.
3 // This code is published under the Eclipse Public License.
4 //
5 // $Id$
6 //
7 // Authors:  Andreas Waechter            IBM    2008-09-19
8 
9 #include "IpInexactNormalTerminationTester.hpp"
10 #include "IpBlas.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 
22 namespace Ipopt
23 {
24 #if COIN_IPOPT_VERBOSITY > 0
25   static const Index dbg_verbosity = 0;
26 #endif
27 
InexactNormalTerminationTester()28   InexactNormalTerminationTester::InexactNormalTerminationTester()
29   {
30     DBG_START_METH("InexactNormalTerminationTester::InexactNormalTerminationTester",dbg_verbosity);
31   }
32 
~InexactNormalTerminationTester()33   InexactNormalTerminationTester::~InexactNormalTerminationTester()
34   {
35     DBG_START_METH("InexactNormalTerminationTester::~InexactNormalTerminationTester()",dbg_verbosity);
36   }
37 
RegisterOptions(SmartPtr<RegisteredOptions> roptions)38   void InexactNormalTerminationTester::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
39   {
40     roptions->AddLowerBoundedNumberOption(
41       "inexact_normal_tol",
42       "Desired relative residual tolerance for iterative solver during normal step computation.",
43       0.0, true, 1e-3,
44       "");
45     roptions->AddLowerBoundedIntegerOption(
46       "inexact_normal_max_iter",
47       "Maximal number of iterative solver iterations during normal step computation",
48       0, 200,
49       "");
50   }
51 
52 
InitializeImpl(const OptionsList & options,const std::string & prefix)53   bool InexactNormalTerminationTester::InitializeImpl(const OptionsList& options,
54       const std::string& prefix)
55   {
56     options.GetNumericValue("inexact_normal_tol", inexact_normal_tol_, prefix);
57     options.GetIntegerValue("inexact_normal_max_iter",
58                             inexact_normal_max_iter_, prefix);
59 
60     std::string inexact_linear_system_scaling;
61     options.GetStringValue("inexact_linear_system_scaling",
62                            inexact_linear_system_scaling, prefix);
63     if (inexact_linear_system_scaling=="slack-based") {
64       requires_scaling_ = true;
65     }
66     else {
67       requires_scaling_ = false;
68     }
69 
70     c_Avc_norm_cauchy_ = -1;
71     return true;
72   }
73 
InitializeSolve()74   bool InexactNormalTerminationTester::InitializeSolve()
75   {
76     DBG_START_METH("InexactNormalTerminationTester::InitializeSolve",
77                    dbg_verbosity);
78 
79 
80     return true;
81   }
82 
83   InexactNormalTerminationTester::ETerminationTest
84   InexactNormalTerminationTester::
TestTermination(Index ndim,const Number * sol,const Number * resid,Index iter,Number norm2_rhs)85   TestTermination(Index ndim, const Number* sol, const Number* resid,
86                   Index iter, Number norm2_rhs)
87   {
88     DBG_START_METH("InexactNormalTerminationTester::TestTermination",
89                    dbg_verbosity);
90 
91     last_iter_ = iter;
92 
93     DBG_ASSERT(c_Avc_norm_cauchy_ >= 0.);
94 
95     ETerminationTest retval = CONTINUE;
96 
97     double norm2_resid = IpBlasDnrm2(ndim, resid, 1);
98     Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
99                    "TTNormal: iter = %d ||resid|| = %23.16e ||rhs|| = %23.16e\n",
100                    iter, norm2_resid, norm2_rhs);
101 
102     if (iter > inexact_normal_max_iter_) {
103       IpData().Append_info_string("N+");
104       retval = OTHER_SATISFIED;
105       return retval;
106     }
107 
108     //if (Min(norm2_resid/norm2_rhs,norm2_resid) > inexact_normal_tol_) {
109     if (norm2_resid/norm2_rhs > inexact_normal_tol_ && norm2_resid > 1e-10) {
110       return retval;
111     }
112 
113     // compute ||c+Av|| for current iterative solution v
114     // note that the sol_x and sol_s have the wrong sign
115     SmartPtr<const Vector> sol_x;
116     SmartPtr<const Vector> sol_s;
117     SmartPtr<const Vector> sol_c;
118     SmartPtr<const Vector> sol_d;
119     GetVectors(ndim, sol, sol_x, sol_s, sol_c, sol_d);
120 
121     if (requires_scaling_) {
122       SmartPtr<const Vector> scaling_vec = InexCq().curr_scaling_slacks();
123       SmartPtr<Vector> tmp = sol_s->MakeNewCopy();
124       tmp->ElementWiseMultiply(*scaling_vec);
125       sol_s = ConstPtr(tmp);
126     }
127 
128     SmartPtr<Vector> inf_c = IpCq().curr_c()->MakeNewCopy();
129     IpCq().curr_jac_c()->MultVector(-1., *sol_x, 1., *inf_c);
130     SmartPtr<const Vector> curr_d_minus_s = IpCq().curr_d_minus_s();
131     SmartPtr<Vector> inf_d = curr_d_minus_s->MakeNew();
132     inf_d->AddTwoVectors(1., *curr_d_minus_s, 1., *sol_s, 0.);
133     IpCq().curr_jac_d()->MultVector(-1., *sol_x, 1., *inf_d);
134 
135     Number trial_c_Av_norm = IpCq().CalcNormOfType(NORM_2, *inf_c, *inf_d);
136 
137     Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
138                    "TTNormal: c_Avc_norm_cauchy = %23.16e trial_c_Av_norm = %23.16e\n", c_Avc_norm_cauchy_, trial_c_Av_norm);
139     Number BasVal = Max(1., IpCq().curr_primal_infeasibility(NORM_2));
140 
141     if (Compare_le(trial_c_Av_norm, c_Avc_norm_cauchy_, BasVal)) {
142       retval = OTHER_SATISFIED;
143     }
144 
145     return retval;
146   }
147 
148   void
Clear()149   InexactNormalTerminationTester::Clear()
150   {
151     DBG_START_METH("InexactNormalTerminationTester::Clear",
152                    dbg_verbosity);
153   }
154 
155 } // namespace Ipopt
156