1 // Copyright (C) 2008, 2011 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 "IpInexactPDTerminationTester.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 InexactPDTerminationTester()28 InexactPDTerminationTester::InexactPDTerminationTester() 29 { 30 DBG_START_METH("InexactPDTerminationTester::InexactPDTerminationTester",dbg_verbosity); 31 } 32 ~InexactPDTerminationTester()33 InexactPDTerminationTester::~InexactPDTerminationTester() 34 { 35 DBG_START_METH("InexactPDTerminationTester::~InexactPDTerminationTester()",dbg_verbosity); 36 } 37 RegisterOptions(SmartPtr<RegisteredOptions> roptions)38 void InexactPDTerminationTester::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 39 { 40 roptions->AddLowerBoundedNumberOption( 41 "tcc_psi", 42 "Psi factor in Tangential Component Condition.", 43 0.0, true, 44 1e-1, 45 ""); 46 roptions->AddLowerBoundedNumberOption( 47 "tcc_theta", 48 "theta factor in Tangential Component Condition.", 49 0.0, true, 50 1e-12, 51 ""); 52 roptions->AddLowerBoundedNumberOption( 53 "tcc_theta_mu_exponent", 54 "exponent for mu when multiplied with tcc_theta in Tangential Component Condition.", 55 0., false, 56 0., 57 ""); 58 roptions->AddLowerBoundedNumberOption( 59 "tcc_zeta", 60 "zeta factor in Tangential Component Condition.", 61 0.0, true, 62 1e-1, 63 ""); 64 roptions->AddLowerBoundedNumberOption( 65 "tt_kappa1", 66 "kappa1 factor in Termination Test 1 and 3.", 67 0.0, true, 68 1e-3, 69 ""); 70 roptions->AddLowerBoundedNumberOption( 71 "tt_kappa2", 72 "kappa2 factor in Termination Test 2.", 73 0.0, true, 74 1e-1, 75 ""); 76 roptions->AddLowerBoundedNumberOption( 77 "tt_eps2", 78 "eps2 factor in Termination Test 2.", 79 0.0, true, 80 1., 81 ""); 82 roptions->AddLowerBoundedNumberOption( 83 "tt_eps3", 84 "eps3 factor in Termination Test 3.", 85 0.0, true, 86 1.-1e-1, 87 ""); 88 roptions->AddLowerBoundedNumberOption( 89 "inexact_desired_pd_residual", 90 "Desired relative residual tolerance for iterative solver during primal-dual step computation.", 91 0.0, true, 1e-3, 92 ""); 93 roptions->AddLowerBoundedIntegerOption( 94 "inexact_desired_pd_residual_iter", 95 "Number of iterations willing to be spent in obtaining desired primal-dual ration.", 96 0, 1, 97 ""); 98 } 99 100 InitializeImpl(const OptionsList & options,const std::string & prefix)101 bool InexactPDTerminationTester::InitializeImpl(const OptionsList& options, 102 const std::string& prefix) 103 { 104 options.GetNumericValue("tcc_psi", tcc_psi_, prefix); 105 options.GetNumericValue("tcc_theta", tcc_theta_, prefix); 106 options.GetNumericValue("tcc_theta_mu_exponent", tcc_theta_mu_exponent_, prefix); 107 options.GetNumericValue("tcc_zeta", tcc_zeta_, prefix); 108 options.GetNumericValue("tt_kappa1", tt_kappa1_, prefix); 109 options.GetNumericValue("tt_kappa2", tt_kappa2_, prefix); 110 options.GetNumericValue("tt_eps2", tt_eps2_, prefix); 111 options.GetNumericValue("tt_eps3", tt_eps3_, prefix); 112 options.GetNumericValue("rho", rho_, prefix); 113 options.GetNumericValue("inexact_desired_pd_residual", 114 inexact_desired_pd_residual_, prefix); 115 options.GetIntegerValue("inexact_desired_pd_residual_iter", 116 inexact_desired_pd_residual_iter_, prefix); 117 118 std::string inexact_linear_system_scaling; 119 options.GetStringValue("inexact_linear_system_scaling", 120 inexact_linear_system_scaling, prefix); 121 if (inexact_linear_system_scaling=="slack-based") { 122 requires_scaling_ = true; 123 } 124 else { 125 requires_scaling_ = false; 126 } 127 128 return true; 129 } 130 InitializeSolve()131 bool InexactPDTerminationTester::InitializeSolve() 132 { 133 DBG_START_METH("InexactPDTerminationTester::InitializeSolve", 134 dbg_verbosity); 135 136 bool compute_normal = InexData().compute_normal(); 137 138 // compute the current infeasibility 139 c_norm_ = IpCq().curr_primal_infeasibility(NORM_2); 140 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 141 "TT: c_norm = %23.16e\n", c_norm_); 142 143 SmartPtr<const Vector> normal_x = InexData().normal_x(); 144 SmartPtr<const Vector> normal_s = InexData().normal_s(); 145 if (compute_normal) { 146 // calculate Jacobian times normal step (no scaling relevant in this space) 147 curr_Av_c_ = InexCq().curr_jac_times_normal_c(); 148 curr_Av_d_ = InexCq().curr_jac_times_normal_d(); 149 150 // compute the linearized infeasibility at the normal step 151 SmartPtr<const Vector> curr_c = IpCq().curr_c(); 152 SmartPtr<Vector> tmp1 = curr_c->MakeNew(); 153 tmp1->AddTwoVectors(1, *curr_c, 1., *curr_Av_c_, 0.); 154 SmartPtr<const Vector> curr_d_minus_s = IpCq().curr_d_minus_s(); 155 SmartPtr<Vector> tmp2 = curr_d_minus_s->MakeNew(); 156 tmp2->AddTwoVectors(1, *curr_d_minus_s, 1., *curr_Av_d_, 0.); 157 c_plus_Av_norm_ = IpCq().CalcNormOfType(NORM_2, *tmp1, *tmp2); 158 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 159 "TT: c_plus_Av_norm_ = %23.16e\n", c_plus_Av_norm_); 160 161 // compute norm of the normal step in the scaled space 162 v_norm_scaled_ = InexCq().slack_scaled_norm(*normal_x, *normal_s); 163 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 164 "TT: v_norm_scaled_ = %23.16e\n", v_norm_scaled_); 165 166 // compute Wv (Hessian times normal step) in the unscaled space 167 curr_Wv_x_ = InexCq().curr_W_times_vec_x(*normal_x); 168 curr_Wv_s_ = InexCq().curr_W_times_vec_s(*normal_s); 169 } 170 else { 171 curr_Av_c_ = NULL; 172 curr_Av_d_ = NULL; 173 c_plus_Av_norm_ = -1.; 174 v_norm_scaled_ = -1.; 175 curr_Wv_x_ = NULL; 176 curr_Wv_s_ = NULL; 177 } 178 179 // store the previous gradient and Jacobian information 180 SmartPtr<const Vector> last_grad_barrier_obj_x = curr_grad_barrier_obj_x_; 181 SmartPtr<const Vector> last_grad_barrier_obj_s = curr_grad_barrier_obj_s_; 182 SmartPtr<const Matrix> last_jac_c = curr_jac_c_; 183 SmartPtr<const Matrix> last_jac_d = curr_jac_d_; 184 SmartPtr<const Vector> last_scaling_slacks = curr_scaling_slacks_; 185 last_Av_norm_ = curr_Av_norm_; 186 187 // get the current gradient and Jacobian information 188 curr_grad_barrier_obj_x_ = IpCq().curr_grad_barrier_obj_x(); 189 curr_grad_barrier_obj_s_ = IpCq().curr_grad_barrier_obj_s(); // (unscaled) 190 curr_jac_c_ = IpCq().curr_jac_c(); 191 curr_jac_d_ = IpCq().curr_jac_d(); 192 curr_scaling_slacks_ = InexCq().curr_scaling_slacks(); 193 194 // calculate \nabla phi(x_{k}) + A_{k}^Ty_k (in scaled space) 195 SmartPtr<const Vector> curr_jac_cT_times_curr_y_c = 196 IpCq().curr_jac_cT_times_curr_y_c(); 197 SmartPtr<const Vector> curr_jac_cT_times_curr_y_d = 198 IpCq().curr_jac_dT_times_curr_y_d(); 199 curr_nabla_phi_plus_ATy_x_ = curr_grad_barrier_obj_x_->MakeNewCopy(); 200 curr_nabla_phi_plus_ATy_x_->AddTwoVectors(1., *curr_jac_cT_times_curr_y_c, 201 1., *curr_jac_cT_times_curr_y_d, 1.); 202 curr_nabla_phi_plus_ATy_s_ = curr_grad_barrier_obj_s_->MakeNew(); 203 curr_nabla_phi_plus_ATy_s_->AddTwoVectors(1., *curr_grad_barrier_obj_s_, 204 -1., *IpData().curr()->y_d(), 0.); 205 curr_nabla_phi_plus_ATy_s_->ElementWiseMultiply(*curr_scaling_slacks_); 206 207 // calculate norms appearing in termination tests 208 curr_tt2_norm_ = IpCq().CalcNormOfType(NORM_2, *curr_nabla_phi_plus_ATy_x_, 209 *curr_nabla_phi_plus_ATy_s_); 210 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 211 "TT: curr_tt2_norm_ = %23.16e\n", curr_tt2_norm_); 212 if (compute_normal) { 213 curr_Av_norm_ = IpCq().CalcNormOfType(NORM_2, *curr_Av_c_, *curr_Av_d_); 214 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 215 "TT: curr_Av_norm_ = %23.16e\n", curr_Av_norm_); 216 curr_tt1_norm_ = sqrt(pow(curr_tt2_norm_, 2) + pow(curr_Av_norm_, 2)); 217 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 218 "TT: curr_tt1_norm_ = %23.16e\n", curr_tt1_norm_); 219 } 220 else { 221 curr_Av_norm_ = -1.; 222 curr_tt1_norm_ = sqrt(pow(curr_tt2_norm_, 2) + pow(c_norm_, 2)); 223 } 224 225 if (compute_normal && IsValid(last_grad_barrier_obj_x)) { 226 // calculate \nabla phi(x_{k-1}) + A_{k-1}^Ty_k (in scaled space) 227 SmartPtr<Vector> last_nabla_phi_plus_ATy_x = 228 last_grad_barrier_obj_x->MakeNewCopy(); 229 last_jac_c->TransMultVector(1., *IpData().curr()->y_c(), 230 1., *last_nabla_phi_plus_ATy_x); 231 last_jac_d->TransMultVector(1., *IpData().curr()->y_d(), 232 1., *last_nabla_phi_plus_ATy_x); 233 SmartPtr<Vector> last_nabla_phi_plus_ATy_s = 234 last_grad_barrier_obj_s->MakeNew(); 235 last_nabla_phi_plus_ATy_s->AddTwoVectors(1., *last_grad_barrier_obj_s, 236 -1., *IpData().curr()->y_d(), 0.); 237 last_nabla_phi_plus_ATy_s->ElementWiseMultiply(*last_scaling_slacks); 238 last_tt1_norm_ = 239 IpCq().CalcNormOfType(NORM_2, *last_nabla_phi_plus_ATy_x, 240 *last_nabla_phi_plus_ATy_s); 241 last_tt1_norm_ = sqrt(pow(last_tt1_norm_, 2) + pow(last_Av_norm_, 2)); 242 } 243 else { 244 last_tt1_norm_= 1e100; 245 } 246 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 247 "TT: last_tt1_norm_ = %23.16e\n", last_tt1_norm_); 248 249 // check if we need to test termination test 2 250 Number ATc_norm = InexCq().curr_scaled_Ac_norm(); 251 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 252 "TT: current ATc norm = %23.16e\n", ATc_norm); 253 try_tt2_ = (ATc_norm <= tt_eps2_*curr_tt2_norm_); 254 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 255 "TT: will%s try termination test 2.\n", try_tt2_ ? "" : " not"); 256 257 return true; 258 } 259 260 InexactPDTerminationTester::ETerminationTest 261 InexactPDTerminationTester:: TestTermination(Index ndim,const Number * sol,const Number * resid,Index iter,Number norm2_rhs)262 TestTermination(Index ndim, const Number* sol, const Number* resid, 263 Index iter, Number norm2_rhs) 264 { 265 DBG_START_METH("InexactPDTerminationTester::TestTermination", 266 dbg_verbosity); 267 268 bool compute_normal = InexData().compute_normal(); 269 270 last_iter_ = iter; 271 272 ETerminationTest retval = CONTINUE; 273 274 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 275 "Starting PD Termination Tester for iteration %d.\n", iter); 276 /* 277 if (iter%5 != 4) { 278 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 279 "TT: immediately leaving tester for iteration %d.\n", iter); 280 return retval; 281 } 282 */ 283 Number norm2_resid = IpBlasDnrm2(ndim, resid, 1); 284 Number test_ratio = norm2_resid/norm2_rhs; // Min(norm2_resid/norm2_rhs, norm2_resid); 285 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 286 "TT: test ratio %e (norm2_rhs = %e norm2_resid = %e).\n", 287 test_ratio, norm2_rhs, norm2_resid); 288 if (iter < inexact_desired_pd_residual_iter_ && 289 test_ratio > inexact_desired_pd_residual_) { 290 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 291 "TT: immediately leaving tester with test ratio %e (norm2_rhs = %e norm2_resid = %e).\n", 292 test_ratio, norm2_rhs, norm2_resid); 293 return retval; 294 } 295 296 SmartPtr<const Vector> sol_x; 297 SmartPtr<const Vector> sol_s; 298 SmartPtr<const Vector> sol_c; 299 SmartPtr<const Vector> sol_d; 300 GetVectors(ndim, sol, sol_x, sol_s, sol_c, sol_d); 301 302 DBG_PRINT_VECTOR(2, "sol_x", *sol_x); 303 DBG_PRINT_VECTOR(2, "sol_s", *sol_s); 304 DBG_PRINT_VECTOR(2, "sol_c", *sol_c); 305 DBG_PRINT_VECTOR(2, "sol_d", *sol_d); 306 307 SmartPtr<const Vector> resid_x; 308 SmartPtr<const Vector> resid_s; 309 SmartPtr<const Vector> resid_c; 310 SmartPtr<const Vector> resid_d; 311 GetVectors(ndim, resid, resid_x, resid_s, resid_c, resid_d); 312 313 if (requires_scaling_) { 314 SmartPtr<const Vector> scaling_vec = curr_scaling_slacks_; 315 SmartPtr<Vector> tmp = sol_s->MakeNewCopy(); 316 tmp->ElementWiseMultiply(*scaling_vec); 317 sol_s = ConstPtr(tmp); 318 tmp = resid_s->MakeNewCopy(); 319 tmp->ElementWiseDivide(*scaling_vec); 320 resid_s = ConstPtr(tmp); 321 } 322 323 //// Set algorithm 324 if (!compute_normal) { 325 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "RUNNING TERMINATION TESTS FOR INEXACT NEWTON (which means u=d)\n"); 326 } 327 else { 328 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "RUNNING TERMINATION TESTS FOR INEXACT NEWTON - TRUST REGION\n"); 329 } 330 331 // Get the tangential step and its scaled norm 332 SmartPtr<const Vector> tangential_x; 333 SmartPtr<const Vector> tangential_s; 334 if (compute_normal) { 335 SmartPtr<const Vector> normal_x = InexData().normal_x(); 336 SmartPtr<Vector> tmp = sol_x->MakeNew(); 337 tmp->AddTwoVectors(1., *sol_x, -1, *normal_x, 0.); 338 tangential_x = ConstPtr(tmp); 339 SmartPtr<const Vector> normal_s = InexData().normal_s(); 340 tmp = sol_s->MakeNew(); 341 tmp->AddTwoVectors(1., *sol_s, -1, *normal_s, 0.); 342 tangential_s = ConstPtr(tmp); 343 } 344 else { 345 tangential_x = sol_x; 346 tangential_s = sol_s; 347 } 348 Number u_norm_scaled = 349 InexCq().slack_scaled_norm(*tangential_x, *tangential_s); 350 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 351 "TT: u_norm_scaled = %23.16e\n", u_norm_scaled); 352 353 // Compute u^TWu 354 DBG_PRINT_VECTOR(2, "tangential_x", *tangential_x); 355 DBG_PRINT_VECTOR(2, "tangential_s", *tangential_s); 356 SmartPtr<const Vector> Wu_x = InexCq().curr_W_times_vec_x(*tangential_x); 357 SmartPtr<const Vector> Wu_s = InexCq().curr_W_times_vec_s(*tangential_s); 358 DBG_PRINT_VECTOR(2, "Wu_x", *Wu_x); 359 DBG_PRINT_VECTOR(2, "Wu_s", *Wu_s); 360 Number uWu = Wu_x->Dot(*tangential_x) + Wu_s->Dot(*tangential_s); 361 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 362 "TT: uWu = %23.16e\n", uWu); 363 364 // Compute norm of c + Ad 365 SmartPtr<Vector> c_plus_Ad_c = IpCq().curr_c()->MakeNewCopy(); 366 curr_jac_c_->MultVector(1., *sol_x, 1., *c_plus_Ad_c); 367 SmartPtr<const Vector> curr_d_minus_s = IpCq().curr_d_minus_s(); 368 SmartPtr<Vector> c_plus_Ad_d = curr_d_minus_s->MakeNew(); 369 c_plus_Ad_d->AddTwoVectors(1., *curr_d_minus_s, -1., *sol_s, 0.); 370 curr_jac_d_->MultVector(1., *sol_x, 1., *c_plus_Ad_d); 371 Number c_plus_Ad_norm = 372 IpCq().CalcNormOfType(NORM_2, *c_plus_Ad_c, *c_plus_Ad_d); 373 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 374 "TT: c_plus_Ad_norm = %23.16e\n", c_plus_Ad_norm); 375 376 // Compute norm of scaled residual rho 377 SmartPtr<Vector> tmp = resid_s->MakeNewCopy(); 378 tmp->ElementWiseMultiply(*curr_scaling_slacks_); 379 Number rho_norm = IpCq().CalcNormOfType(NORM_2, *resid_x, *tmp); 380 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 381 "TT: rho_norm = %23.16e\n", rho_norm); 382 tmp = NULL; 383 384 // TODO: AW wants to discuss with Frank 385 Number Upsilon = -1.; 386 Number Nu = -1.; 387 if (!compute_normal) { 388 #if 0 389 // Compute Nu = ||A*u||^2/||A||^2 390 SmartPtr<const Vector> curr_Au_c = IpCq().curr_jac_c_times_vec(*tangential_x); 391 SmartPtr<Vector> curr_Au_d = sol_s->MakeNew(); 392 curr_Au_d->AddTwoVectors(1., *IpCq().curr_jac_d_times_vec(*tangential_x), -1., *tangential_s, 0.); 393 Number Nu = IpCq().CalcNormOfType(NORM_2, *curr_Au_c, *curr_Au_d); 394 Nu = pow(Nu,2); 395 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: ||A*u||^2 = %23.16e\n", Nu); 396 Number A_norm2 = InexCq().curr_scaled_A_norm2(); 397 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: ||A||^2 = %23.16e\n", A_norm2); 398 #endif 399 Nu = 0;//Nu/A_norm2; 400 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: Nu = ||A*u||^2/||A||^2 = %23.16e\n", Nu); 401 402 // Compute Upsilon = ||u||^2 - Nu 403 Upsilon = u_norm_scaled*u_norm_scaled - Nu; 404 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: Upsilon = ||u||^2 - ||A*u||^2/||A||^2 = %23.16e\n", Upsilon); 405 } 406 407 // Base value, something on the order of square root of machine epsilon; TODO: find a better base value 408 Number BasVal = Max(IpData().curr()->x()->Amax(), IpData().curr()->s()->Amax()); 409 410 // Check tangential component condition, part 1 411 Number lhs; 412 Number rhs; 413 if (!compute_normal) { 414 lhs = Upsilon; 415 rhs = pow(tcc_psi_,2)*Nu; 416 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 417 "TCC1 testing Upsilon(=%23.16e) <= (tcc_psi_^2)*Nu(=%23.16e) --> ", lhs, rhs); 418 } 419 else { 420 lhs = u_norm_scaled; 421 rhs = tcc_psi_*v_norm_scaled_; 422 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 423 "TCC1 testing u_norm_scaled(=%23.16e) <= tcc_psi_*v_norm_scaled(=%23.16e) --> ", lhs, rhs); 424 } 425 bool tcc1 = Compare_le(lhs, rhs, BasVal); 426 if (tcc1) { 427 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 428 } 429 else { 430 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 431 } 432 433 // Check tangential component condition, part 2a 434 const Number mu = IpData().curr_mu(); 435 rhs = 0.5*uWu; 436 if (!compute_normal) { 437 lhs = tcc_theta_*pow(mu,tcc_theta_mu_exponent_)*Upsilon; 438 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 439 "TCC2a testing 0.5*uWu(=%23.16e) >= tcc_theta_*pow(mu,tcc_theta_mu_exponent_)*Upsilon(=%23.16e) -->", rhs, lhs); 440 } 441 else { 442 lhs = tcc_theta_*pow(mu,tcc_theta_mu_exponent_)*pow(u_norm_scaled, 2); 443 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 444 "TCC2a testing 0.5*uWu(=%23.16e) >= tcc_theta_*pow(mu,tcc_theta_mu_exponent_)*u_norm^2(=%23.16e) -->", rhs, lhs); 445 } 446 bool tcc2a = Compare_le(lhs, rhs, BasVal); 447 if (tcc2a) { 448 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 449 } 450 else { 451 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 452 } 453 454 // Check tangential component condition, part 2b (only in MIPS) 455 bool tcc = tcc1; 456 if (!tcc && tcc2a) { 457 if (!compute_normal) { 458 tcc = tcc2a; 459 } 460 else { 461 lhs = 0.5*uWu + curr_grad_barrier_obj_x_->Dot(*tangential_x) + curr_grad_barrier_obj_s_->Dot(*tangential_s) + curr_Wv_x_->Dot(*tangential_x) + curr_Wv_s_->Dot(*tangential_s); 462 rhs = tcc_zeta_*v_norm_scaled_; 463 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 464 "TCC2b testing (grad_barr^Tu + v^TWu + 0.5*uWu)(=%23.16e) <= tcc_zeta_*v_norm(=%23.16e) -->", lhs, rhs); 465 tcc = Compare_le(lhs, rhs, BasVal); 466 if (tcc) { 467 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 468 } 469 else { 470 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 471 } 472 } 473 } 474 475 // Check tangential component condition 476 if (tcc) { 477 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Tangential Component Condition satisfied\n"); 478 } 479 else { 480 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Tangential Component Condition violated\n"); 481 } 482 483 // Check termination test 1, residual condition 484 bool tt1 = tcc; 485 bool tt1_kappa1 = tcc; 486 if (!compute_normal) { 487 // Compute scaled norm of entire residual in case there is no step 488 // decomposition. In that case, c_plus_Ad_norm should indeed be 489 // the same as what resid_c and resid_d woulod give (TODO: 490 // check?!?) 491 Number resid_norm = sqrt(pow(rho_norm, 2) + pow(c_plus_Ad_norm, 2)); 492 lhs = resid_norm; 493 rhs = tt_kappa1_*curr_tt1_norm_; 494 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 495 "TT1 testing resid_norm(=%23.16e) <= tt_kappa1_*curr_tt1_norm_(=%23.16e) --> ", lhs, rhs); 496 tt1 = Compare_le(lhs, rhs, BasVal); 497 } 498 else if (tt1) { 499 lhs = rho_norm; 500 rhs = tt_kappa1_*Min(curr_tt1_norm_, last_tt1_norm_); 501 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 502 "TT1 testing rho_norm(=%23.16e) <= kappa1*min(curr_tt1_norm_, last_tt1_norm_)(=%23.16e) -->", lhs, rhs); 503 tt1_kappa1 = Compare_le(lhs, rhs, BasVal); 504 tt1 = tt1_kappa1; 505 } 506 if (tt1) { 507 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 508 } 509 else { 510 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 511 } 512 513 // Check termination test 1, model reduction condition 514 bool model_reduction = false; 515 if (!compute_normal || tt1) { 516 Number curr_nu = InexData().curr_nu(); 517 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: curr_nu = %23.16e\n", curr_nu); 518 Number delta_m = -(curr_grad_barrier_obj_x_->Dot(*sol_x) + curr_grad_barrier_obj_s_->Dot(*sol_s)); 519 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: -grad_barr^Td = %23.16e\n", delta_m); 520 delta_m += curr_nu*(c_norm_ - c_plus_Ad_norm); 521 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "TT: delta_m = %23.16e\n", delta_m); 522 rhs = delta_m; 523 Number sigma = rho_*tt_eps3_; 524 if (!compute_normal) { 525 lhs = Max(0.5*uWu, tcc_theta_*Upsilon) + sigma*curr_nu*Max(c_norm_, c_plus_Ad_norm - c_norm_); 526 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 527 "MRC testing delta_m(=%23.16e) >= max(0.5*uWu,tcc_theta_*Upsilon) + sigma*nu*max(c_norm_, c_plus_Ad_norm - c_norm_)(=%23.16e) -->", rhs, lhs); 528 model_reduction = Compare_le(lhs, rhs, BasVal); 529 if (model_reduction) { 530 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 531 } 532 else { 533 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 534 } 535 if (tt1) { 536 tt1 = model_reduction; 537 } 538 } 539 else { 540 lhs = Max(0.5*uWu, tcc_theta_*pow(u_norm_scaled, 2)) + sigma*curr_nu*(c_norm_ - c_plus_Av_norm_); 541 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 542 "MRC testing delta_m(=%23.16e) >= max(0.5*uWu,tcc_theta_*u_norm^2) + sigma*nu*(c_norm_ - c_plus_Av_norm_)(=%23.16e) -->", rhs, lhs); 543 tt1 = Compare_le(lhs, rhs, BasVal); 544 if (tt1) { 545 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 546 } 547 else { 548 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 549 } 550 } 551 } 552 553 // Check termination test 1 554 if (tt1) { 555 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 1 satisfied.\n"); 556 return TEST_1_SATISFIED; 557 } 558 else { 559 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 1 not satisfied.\n"); 560 } 561 562 // Check termination test 3, residual condition 563 bool tt3 = tcc; 564 if (tt3) { 565 if (!compute_normal) { 566 lhs = rho_norm; 567 rhs = tt_kappa1_*c_norm_; 568 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 569 "TT3 testing rho_norm(=%23.16e) <= tt_kappa1_*c_norm_(=%23.16e) -->", lhs, rhs); 570 tt3 = Compare_le(lhs, rhs, BasVal); 571 } 572 else { 573 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 574 "TT3 with residual condition from TT1 -->"); 575 tt3 = tt1_kappa1; 576 } 577 if (tt3) { 578 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 579 } 580 else { 581 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 582 } 583 } 584 585 // Check termination test 3, linearized feasibility condition 586 if (tt3) { 587 if (!compute_normal) { 588 lhs = c_plus_Ad_norm; 589 rhs = tt_kappa1_*c_norm_; 590 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 591 "TT3 testing c_plus_Ad_norm(=%23.16e) <= tt_kappa1_*c_norm(=%23.16e) -->", lhs, rhs); 592 tt3 = Compare_le(lhs, rhs, BasVal); 593 } 594 else { 595 lhs = tt_eps3_*(c_norm_ - c_plus_Av_norm_); 596 rhs = c_norm_ - c_plus_Ad_norm; 597 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 598 "TT3 testing (c_norm_-c_plus_Ad_norm)(=%23.16e) >= eps3*(c_norm_-c_plus_Av_norm_)(=%23.16e) -->", rhs, lhs); 599 tt3 = Compare_le(lhs, rhs, BasVal); 600 } 601 if (tt3) { 602 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 603 } 604 else { 605 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 606 } 607 } 608 609 // Check termination test 3 610 if (tt3) { 611 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 3 satisfied.\n"); 612 return TEST_3_SATISFIED; 613 } 614 else { 615 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 3 not satisfied.\n"); 616 } 617 618 // Check termination test 2 619 bool tt2 = try_tt2_; 620 if (tt2) { 621 DBG_PRINT_VECTOR(2, "curr_nabla_phi_plus_ATy_x_", *curr_nabla_phi_plus_ATy_x_); 622 DBG_PRINT_VECTOR(2, "curr_nabla_phi_plus_ATy_s_", *curr_nabla_phi_plus_ATy_s_); 623 DBG_PRINT_VECTOR(2, "sol_c", *sol_c); 624 DBG_PRINT_VECTOR(2, "sol_d", *sol_d); 625 SmartPtr<Vector> sol_d_scaled = sol_d->MakeNewCopy(); 626 sol_d_scaled->ElementWiseMultiply(*curr_scaling_slacks_); 627 DBG_PRINT_VECTOR(2, "sol_d_scaled", *sol_d_scaled); 628 SmartPtr<Vector> nabla_phi_plus_ATydelta_x = curr_nabla_phi_plus_ATy_x_->MakeNewCopy(); 629 curr_jac_c_->TransMultVector(1., *sol_c, 1., *curr_nabla_phi_plus_ATy_x_); 630 curr_jac_d_->TransMultVector(1., *sol_d, 1., *curr_nabla_phi_plus_ATy_x_); 631 SmartPtr<Vector> nabla_phi_plus_ATydelta_s = curr_nabla_phi_plus_ATy_s_->MakeNew(); 632 nabla_phi_plus_ATydelta_s->AddTwoVectors(1., *curr_nabla_phi_plus_ATy_s_, -1., *sol_d_scaled, 0.); 633 Number nabla_phi_plus_ATydelta_norm = IpCq().CalcNormOfType(NORM_2, *nabla_phi_plus_ATydelta_x, *nabla_phi_plus_ATydelta_s); 634 lhs = nabla_phi_plus_ATydelta_norm; 635 if (!compute_normal) { 636 rhs = tt_kappa2_*curr_tt2_norm_; 637 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 638 "TT2 testing ||gamma+A^T(y+delta)||(=%23.16e) <= kappa2*curr_tt2_norm_(=%23.16e) -->", lhs, rhs); 639 } 640 else { 641 rhs = tt_kappa2_*Min(curr_tt2_norm_, last_tt1_norm_); 642 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 643 "TT2 testing ||gamma+A^T(y+delta)||(=%23.16e) <= kappa2*min(curr_tt2_norm_, last_tt1_norm_)(=%23.16e) -->", lhs, rhs); 644 } 645 tt2 = Compare_le(lhs, rhs, BasVal); 646 if (tt2) { 647 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "satisfied\n"); 648 } 649 else { 650 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, "violated\n"); 651 } 652 } 653 654 // Check termination test 2 655 if (tt2) { 656 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 2 satisfied.\n"); 657 return TEST_2_SATISFIED; 658 } 659 else if (try_tt2_) { 660 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Termination Test 2 not satisfied.\n"); 661 } 662 663 // Check if the Hessian should be modified 664 if (tcc1 || tcc2a) {// || (!compute_normal && model_reduction)) { 665 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Hessian Modification not requested.\n"); 666 } 667 else { 668 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, "Hessian Modification requested.\n"); 669 return MODIFY_HESSIAN; 670 } 671 672 return retval; 673 } 674 675 void Clear()676 InexactPDTerminationTester::Clear() 677 { 678 DBG_START_METH("InexactPDTerminationTester::Clear", 679 dbg_verbosity); 680 681 curr_Av_c_ = NULL; 682 curr_Av_d_ = NULL; 683 curr_Wv_x_ = NULL; 684 curr_Wv_s_ = NULL; 685 curr_nabla_phi_plus_ATy_x_ = NULL; 686 curr_nabla_phi_plus_ATy_s_ = NULL; 687 } 688 689 } // namespace Ipopt 690