1 // Copyright (C) 2004, 2007 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: Carl Laird, Andreas Waechter IBM 2004-08-13 8 9 #include "IpPDFullSpaceSolver.hpp" 10 #include "IpDebug.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 25 #if COIN_IPOPT_VERBOSITY > 0 26 static const Index dbg_verbosity = 0; 27 #endif 28 PDFullSpaceSolver(AugSystemSolver & augSysSolver,PDPerturbationHandler & perturbHandler)29 PDFullSpaceSolver::PDFullSpaceSolver(AugSystemSolver& augSysSolver, 30 PDPerturbationHandler& perturbHandler) 31 : 32 PDSystemSolver(), 33 augSysSolver_(&augSysSolver), 34 perturbHandler_(&perturbHandler), 35 dummy_cache_(1) 36 { 37 DBG_START_METH("PDFullSpaceSolver::PDFullSpaceSolver",dbg_verbosity); 38 } 39 ~PDFullSpaceSolver()40 PDFullSpaceSolver::~PDFullSpaceSolver() 41 { 42 DBG_START_METH("PDFullSpaceSolver::~PDFullSpaceSolver()",dbg_verbosity); 43 } 44 RegisterOptions(SmartPtr<RegisteredOptions> roptions)45 void PDFullSpaceSolver::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 46 { 47 roptions->AddLowerBoundedIntegerOption( 48 "min_refinement_steps", 49 "Minimum number of iterative refinement steps per linear system solve.", 50 0, 1, 51 "Iterative refinement (on the full unsymmetric system) is performed for " 52 "each right hand side. This option determines the minimum number " 53 "of iterative refinements (i.e. at least \"min_refinement_steps\" " 54 "iterative refinement steps are enforced per right hand side.)"); 55 roptions->AddLowerBoundedIntegerOption( 56 "max_refinement_steps", 57 "Maximum number of iterative refinement steps per linear system solve.", 58 0, 10, 59 "Iterative refinement (on the full unsymmetric system) is performed for " 60 "each right hand side. This option determines the maximum number " 61 "of iterative refinement steps."); 62 roptions->AddLowerBoundedNumberOption( 63 "residual_ratio_max", 64 "Iterative refinement tolerance", 65 0.0, true, 1e-10, 66 "Iterative refinement is performed until the residual test ratio is " 67 "less than this tolerance (or until \"max_refinement_steps\" refinement " 68 "steps are performed)."); 69 roptions->AddLowerBoundedNumberOption( 70 "residual_ratio_singular", 71 "Threshold for declaring linear system singular after failed iterative refinement.", 72 0.0, true, 1e-5, 73 "If the residual test ratio is larger than this value after failed " 74 "iterative refinement, the algorithm pretends that the linear system is " 75 "singular."); 76 // ToDo Think about following option - are the correct norms used? 77 roptions->AddLowerBoundedNumberOption( 78 "residual_improvement_factor", 79 "Minimal required reduction of residual test ratio in iterative refinement.", 80 0.0, true, 0.999999999, 81 "If the improvement of the residual test ratio made by one iterative " 82 "refinement step is not better than this factor, iterative refinement " 83 "is aborted."); 84 roptions->AddLowerBoundedNumberOption( 85 "neg_curv_test_tol", 86 "Tolerance for heuristic to ignore wrong inertia.", 87 0.0, false, 0.0, 88 "If nonzero, incorrect inertia in the augmented system is ignored, and " 89 "Ipopt tests if the direction is a direction of positive curvature. This " 90 "tolerance is alpha_n in the paper by Zavala and Chiang (2014) and it " 91 "determines when the direction is considered to be sufficiently positive. " 92 "A value in the range of [1e-12, 1e-11] is recommended."); 93 roptions->AddStringOption2( 94 "neg_curv_test_reg", 95 "Whether to do the curvature test with the primal regularization (see Zavala and Chiang, 2014).", 96 "yes", 97 "yes", "use primal regularization with the inertia-free curvature test", 98 "no", "use original IPOPT approach, in which the primal regularization is ignored", 99 ""); 100 } 101 102 InitializeImpl(const OptionsList & options,const std::string & prefix)103 bool PDFullSpaceSolver::InitializeImpl(const OptionsList& options, 104 const std::string& prefix) 105 { 106 // Check for the algorithm options 107 options.GetIntegerValue("min_refinement_steps", min_refinement_steps_, prefix); 108 options.GetIntegerValue("max_refinement_steps", max_refinement_steps_, prefix); 109 ASSERT_EXCEPTION(max_refinement_steps_ >= min_refinement_steps_, OPTION_INVALID, 110 "Option \"max_refinement_steps\": This value must be larger than or equal to min_refinement_steps (default 1)"); 111 112 options.GetNumericValue("residual_ratio_max", residual_ratio_max_, prefix); 113 options.GetNumericValue("residual_ratio_singular", residual_ratio_singular_, prefix); 114 ASSERT_EXCEPTION(residual_ratio_singular_ >= residual_ratio_max_, OPTION_INVALID, 115 "Option \"residual_ratio_singular\": This value must be not smaller than residual_ratio_max."); 116 options.GetNumericValue("residual_improvement_factor", residual_improvement_factor_, prefix); 117 options.GetNumericValue("neg_curv_test_tol", neg_curv_test_tol_, prefix); 118 options.GetBoolValue("neg_curv_test_reg", neg_curv_test_reg_, prefix); 119 120 // Reset internal flags and data 121 augsys_improved_ = false; 122 123 if (!augSysSolver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(), 124 options, prefix)) { 125 return false; 126 } 127 128 return perturbHandler_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(), 129 options, prefix); 130 } 131 Solve(Number alpha,Number beta,const IteratesVector & rhs,IteratesVector & res,bool allow_inexact,bool improve_solution)132 bool PDFullSpaceSolver::Solve(Number alpha, 133 Number beta, 134 const IteratesVector& rhs, 135 IteratesVector& res, 136 bool allow_inexact, 137 bool improve_solution /* = false */) 138 { 139 DBG_START_METH("PDFullSpaceSolver::Solve",dbg_verbosity); 140 DBG_ASSERT(!allow_inexact || !improve_solution); 141 DBG_ASSERT(!improve_solution || beta==0.); 142 143 // Timing of PDSystem solver starts here 144 IpData().TimingStats().PDSystemSolverTotal().Start(); 145 146 DBG_PRINT_VECTOR(2, "rhs_x", *rhs.x()); 147 DBG_PRINT_VECTOR(2, "rhs_s", *rhs.s()); 148 DBG_PRINT_VECTOR(2, "rhs_c", *rhs.y_c()); 149 DBG_PRINT_VECTOR(2, "rhs_d", *rhs.y_d()); 150 DBG_PRINT_VECTOR(2, "rhs_zL", *rhs.z_L()); 151 DBG_PRINT_VECTOR(2, "rhs_zU", *rhs.z_U()); 152 DBG_PRINT_VECTOR(2, "rhs_vL", *rhs.v_L()); 153 DBG_PRINT_VECTOR(2, "rhs_vU", *rhs.v_U()); 154 DBG_PRINT_VECTOR(2, "res_x in", *res.x()); 155 DBG_PRINT_VECTOR(2, "res_s in", *res.s()); 156 DBG_PRINT_VECTOR(2, "res_c in", *res.y_c()); 157 DBG_PRINT_VECTOR(2, "res_d in", *res.y_d()); 158 DBG_PRINT_VECTOR(2, "res_zL in", *res.z_L()); 159 DBG_PRINT_VECTOR(2, "res_zU in", *res.z_U()); 160 DBG_PRINT_VECTOR(2, "res_vL in", *res.v_L()); 161 DBG_PRINT_VECTOR(2, "res_vU in", *res.v_U()); 162 163 // if beta is nonzero, keep a copy of the incoming values in res_ */ 164 SmartPtr<IteratesVector> copy_res; 165 if (beta != 0.) { 166 copy_res = res.MakeNewIteratesVectorCopy(); 167 } 168 169 // Receive data about matrix 170 SmartPtr<const Vector> x = IpData().curr()->x(); 171 SmartPtr<const Vector> s = IpData().curr()->s(); 172 SmartPtr<const SymMatrix> W = IpData().W(); 173 SmartPtr<const Matrix> J_c = IpCq().curr_jac_c(); 174 SmartPtr<const Matrix> J_d = IpCq().curr_jac_d(); 175 SmartPtr<const Matrix> Px_L = IpNLP().Px_L(); 176 SmartPtr<const Matrix> Px_U = IpNLP().Px_U(); 177 SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L(); 178 SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U(); 179 SmartPtr<const Vector> z_L = IpData().curr()->z_L(); 180 SmartPtr<const Vector> z_U = IpData().curr()->z_U(); 181 SmartPtr<const Vector> v_L = IpData().curr()->v_L(); 182 SmartPtr<const Vector> v_U = IpData().curr()->v_U(); 183 SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L(); 184 SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U(); 185 SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L(); 186 SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U(); 187 SmartPtr<const Vector> sigma_x = IpCq().curr_sigma_x(); 188 SmartPtr<const Vector> sigma_s = IpCq().curr_sigma_s(); 189 DBG_PRINT_VECTOR(2, "Sigma_x", *sigma_x); 190 DBG_PRINT_VECTOR(2, "Sigma_s", *sigma_s); 191 192 bool done = false; 193 // The following flag is set to true, if we asked the linear 194 // solver to improve the quality of the solution in 195 // the next solve 196 bool resolve_with_better_quality = false; 197 // the following flag is set to true, if iterative refinement 198 // failed and we want to try if a modified system is able to 199 // remedy that problem by pretending the matrix is singular 200 bool pretend_singular = false; 201 bool pretend_singular_last_time = false; 202 203 // Beginning of loop for solving the system (including all 204 // modifications for the linear system to ensure good solution 205 // quality) 206 while (!done) { 207 208 // if improve_solution is true, we are given already a solution 209 // from the calling function, so we can skip the first solve 210 bool solve_retval = true; 211 if (!improve_solution) { 212 solve_retval = 213 SolveOnce(resolve_with_better_quality, pretend_singular, 214 *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, 215 *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, 216 *sigma_x, *sigma_s, 1., 0., rhs, res); 217 resolve_with_better_quality = false; 218 pretend_singular = false; 219 } 220 improve_solution = false; 221 222 if (!solve_retval) { 223 // If system seems not to be solvable, we return with false 224 // and let the calling routine deal with it. 225 IpData().TimingStats().PDSystemSolverTotal().End(); 226 return false; 227 } 228 229 if (allow_inexact) { 230 // no safety checks required 231 if (Jnlst().ProduceOutput(J_MOREDETAILED, J_LINEAR_ALGEBRA)) { 232 SmartPtr<IteratesVector> resid = res.MakeNewIteratesVector(true); 233 ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, 234 *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, 235 *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, 236 alpha, beta, rhs, res, *resid); 237 } 238 break; 239 } 240 241 // Get space for the residual 242 SmartPtr<IteratesVector> resid = res.MakeNewIteratesVector(true); 243 244 // ToDo don't to that after max refinement? 245 ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, 246 *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, 247 *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, 248 alpha, beta, rhs, res, *resid); 249 250 Number residual_ratio = 251 ComputeResidualRatio(rhs, res, *resid); 252 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 253 "residual_ratio = %e\n", residual_ratio); 254 Number residual_ratio_old = residual_ratio; 255 256 // Beginning of loop for iterative refinement 257 Index num_iter_ref = 0; 258 bool quit_refinement = false; 259 while (!allow_inexact && !quit_refinement && 260 (num_iter_ref < min_refinement_steps_ || 261 residual_ratio > residual_ratio_max_) ) { 262 263 // To the next back solve 264 solve_retval = 265 SolveOnce(resolve_with_better_quality, false, 266 *W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, *z_L, *z_U, 267 *v_L, *v_U, *slack_x_L, *slack_x_U, *slack_s_L, *slack_s_U, 268 *sigma_x, *sigma_s, -1., 1., *resid, res); 269 ASSERT_EXCEPTION(solve_retval, INTERNAL_ABORT, 270 "SolveOnce returns false during iterative refinement."); 271 272 ComputeResiduals(*W, *J_c, *J_d, *Px_L, *Px_U, *Pd_L, *Pd_U, 273 *z_L, *z_U, *v_L, *v_U, *slack_x_L, *slack_x_U, 274 *slack_s_L, *slack_s_U, *sigma_x, *sigma_s, 275 alpha, beta, rhs, res, *resid); 276 277 residual_ratio = 278 ComputeResidualRatio(rhs, res, *resid); 279 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 280 "residual_ratio = %e\n", residual_ratio); 281 282 num_iter_ref++; 283 // Check if we have to give up on iterative refinement 284 if (residual_ratio > residual_ratio_max_ && 285 num_iter_ref>min_refinement_steps_ && 286 (num_iter_ref>max_refinement_steps_ || 287 residual_ratio>residual_improvement_factor_*residual_ratio_old)) { 288 289 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 290 "Iterative refinement failed with residual_ratio = %e\n", residual_ratio); 291 quit_refinement = true; 292 293 // Pretend singularity only once - if it didn't help, we 294 // have to live with what we got so far 295 resolve_with_better_quality = false; 296 DBG_PRINT((1, "pretend_singular = %d\n", pretend_singular)); 297 if (!pretend_singular_last_time) { 298 // First try if we can ask the augmented system solver to 299 // improve the quality of the solution (only if that hasn't 300 // been done before for this linear system) 301 if (!augsys_improved_) { 302 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 303 "Asking augmented system solver to improve quality of its solutions.\n"); 304 augsys_improved_ = augSysSolver_->IncreaseQuality(); 305 if (augsys_improved_) { 306 IpData().Append_info_string("q"); 307 resolve_with_better_quality = true; 308 } 309 else { 310 // solver said it cannot improve quality, so let 311 // possibly conclude that the current modification is 312 // singular 313 pretend_singular = true; 314 } 315 } 316 else { 317 // we had already asked the solver before to improve the 318 // quality of the solution, so let's now pretend that the 319 // modification is possibly singular 320 pretend_singular = true; 321 } 322 pretend_singular_last_time = pretend_singular; 323 if (pretend_singular) { 324 // let's only conclude that the current linear system 325 // including modifications is singular, if the residual is 326 // quite bad 327 if (residual_ratio < residual_ratio_singular_) { 328 pretend_singular = false; 329 IpData().Append_info_string("S"); 330 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 331 "Just accept current solution.\n"); 332 } 333 else { 334 IpData().Append_info_string("s"); 335 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 336 "Pretend that the current system (including modifications) is singular.\n"); 337 } 338 } 339 } 340 else { 341 pretend_singular = false; 342 DBG_PRINT((1,"Resetting pretend_singular to false.\n")); 343 } 344 } 345 346 residual_ratio_old = residual_ratio; 347 } // End of loop for iterative refinement 348 349 done = !(resolve_with_better_quality) && !(pretend_singular); 350 351 } // End of loop for solving the linear system (incl. modifications) 352 353 // Finally let's assemble the res result vectors 354 if (alpha != 0.) { 355 res.Scal(alpha); 356 } 357 358 if (beta != 0.) { 359 res.Axpy(beta, *copy_res); 360 } 361 362 DBG_PRINT_VECTOR(2, "res_x", *res.x()); 363 DBG_PRINT_VECTOR(2, "res_s", *res.s()); 364 DBG_PRINT_VECTOR(2, "res_c", *res.y_c()); 365 DBG_PRINT_VECTOR(2, "res_d", *res.y_d()); 366 DBG_PRINT_VECTOR(2, "res_zL", *res.z_L()); 367 DBG_PRINT_VECTOR(2, "res_zU", *res.z_U()); 368 DBG_PRINT_VECTOR(2, "res_vL", *res.v_L()); 369 DBG_PRINT_VECTOR(2, "res_vU", *res.v_U()); 370 371 IpData().TimingStats().PDSystemSolverTotal().End(); 372 373 return true; 374 } 375 SolveOnce(bool resolve_with_better_quality,bool pretend_singular,const SymMatrix & W,const Matrix & J_c,const Matrix & J_d,const Matrix & Px_L,const Matrix & Px_U,const Matrix & Pd_L,const Matrix & Pd_U,const Vector & z_L,const Vector & z_U,const Vector & v_L,const Vector & v_U,const Vector & slack_x_L,const Vector & slack_x_U,const Vector & slack_s_L,const Vector & slack_s_U,const Vector & sigma_x,const Vector & sigma_s,Number alpha,Number beta,const IteratesVector & rhs,IteratesVector & res)376 bool PDFullSpaceSolver::SolveOnce(bool resolve_with_better_quality, 377 bool pretend_singular, 378 const SymMatrix& W, 379 const Matrix& J_c, 380 const Matrix& J_d, 381 const Matrix& Px_L, 382 const Matrix& Px_U, 383 const Matrix& Pd_L, 384 const Matrix& Pd_U, 385 const Vector& z_L, 386 const Vector& z_U, 387 const Vector& v_L, 388 const Vector& v_U, 389 const Vector& slack_x_L, 390 const Vector& slack_x_U, 391 const Vector& slack_s_L, 392 const Vector& slack_s_U, 393 const Vector& sigma_x, 394 const Vector& sigma_s, 395 Number alpha, 396 Number beta, 397 const IteratesVector& rhs, 398 IteratesVector& res) 399 { 400 // TO DO LIST: 401 // 402 // 1. decide for reasonable return codes (e.g. fatal error, too 403 // ill-conditioned...) 404 // 2. Make constants parameters that can be set from the outside 405 // 3. Get Information out of Ipopt structures 406 // 4. add heuristic for structurally singular problems 407 // 5. see if it makes sense to distinguish delta_x and delta_s, 408 // or delta_c and delta_d 409 // 6. increase pivot tolerance if number of get evals so too small 410 DBG_START_METH("PDFullSpaceSolver::SolveOnce",dbg_verbosity); 411 412 IpData().TimingStats().PDSystemSolverSolveOnce().Start(); 413 414 // Compute the right hand side for the augmented system formulation 415 SmartPtr<Vector> augRhs_x = rhs.x()->MakeNewCopy(); 416 Px_L.AddMSinvZ(1.0, slack_x_L, *rhs.z_L(), *augRhs_x); 417 Px_U.AddMSinvZ(-1.0, slack_x_U, *rhs.z_U(), *augRhs_x); 418 419 SmartPtr<Vector> augRhs_s = rhs.s()->MakeNewCopy(); 420 Pd_L.AddMSinvZ(1.0, slack_s_L, *rhs.v_L(), *augRhs_s); 421 Pd_U.AddMSinvZ(-1.0, slack_s_U, *rhs.v_U(), *augRhs_s); 422 423 // Get space into which we can put the solution of the augmented system 424 SmartPtr<IteratesVector> sol = res.MakeNewIteratesVector(true); 425 426 // Now check whether any data has changed 427 std::vector<const TaggedObject*> deps(13); 428 deps[0] = &W; 429 deps[1] = &J_c; 430 deps[2] = &J_d; 431 deps[3] = &z_L; 432 deps[4] = &z_U; 433 deps[5] = &v_L; 434 deps[6] = &v_U; 435 deps[7] = &slack_x_L; 436 deps[8] = &slack_x_U; 437 deps[9] = &slack_s_L; 438 deps[10] = &slack_s_U; 439 deps[11] = &sigma_x; 440 deps[12] = &sigma_s; 441 void* dummy; 442 bool uptodate = dummy_cache_.GetCachedResult(dummy, deps); 443 if (!uptodate) { 444 dummy_cache_.AddCachedResult(dummy, deps); 445 augsys_improved_ = false; 446 } 447 // improve_current_solution can only be true, if that system has 448 // been solved before 449 DBG_ASSERT((!resolve_with_better_quality && !pretend_singular) || uptodate); 450 451 ESymSolverStatus retval; 452 if (uptodate && !pretend_singular) { 453 454 // Get the perturbation values 455 Number delta_x; 456 Number delta_s; 457 Number delta_c; 458 Number delta_d; 459 perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d); 460 461 // No need to go through the pain of finding the appropriate 462 // values for the deltas, because the matrix hasn't changed since 463 // the last call. So, just call the Solve Method 464 // 465 // Note: resolve_with_better_quality is true, then the Solve 466 // method has already asked the augSysSolver to increase the 467 // quality at the end solve, and we are now getting the solution 468 // with that better quality 469 retval = augSysSolver_->Solve(&W, 1.0, &sigma_x, delta_x, 470 &sigma_s, delta_s, &J_c, NULL, 471 delta_c, &J_d, NULL, delta_d, 472 *augRhs_x, *augRhs_s, *rhs.y_c(), *rhs.y_d(), 473 *sol->x_NonConst(), *sol->s_NonConst(), 474 *sol->y_c_NonConst(), *sol->y_d_NonConst(), 475 false, 0); 476 if (retval!=SYMSOLVER_SUCCESS) { 477 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 478 return false; 479 } 480 } 481 else { 482 const Index numberOfEVals=rhs.y_c()->Dim()+rhs.y_d()->Dim(); 483 // counter for the number of trial evaluations 484 // (ToDo is not at the correct place) 485 Index count = 0; 486 487 // Get the very first perturbation values from the perturbation 488 // Handler 489 Number delta_x; 490 Number delta_s; 491 Number delta_c; 492 Number delta_d; 493 perturbHandler_->ConsiderNewSystem(delta_x, delta_s, delta_c, delta_d); 494 495 retval = SYMSOLVER_SINGULAR; 496 bool fail = false; 497 498 while (retval!= SYMSOLVER_SUCCESS && !fail) { 499 500 if (pretend_singular) { 501 retval = SYMSOLVER_SINGULAR; 502 pretend_singular = false; 503 } 504 else { 505 count++; 506 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 507 "Solving system with delta_x=%e delta_s=%e\n delta_c=%e delta_d=%e\n", 508 delta_x, delta_s, delta_c, delta_d); 509 bool check_inertia = true; 510 if (neg_curv_test_tol_ > 0.) { 511 check_inertia = false; 512 } 513 retval = augSysSolver_->Solve(&W, 1.0, &sigma_x, delta_x, 514 &sigma_s, delta_s, &J_c, NULL, 515 delta_c, &J_d, NULL, delta_d, 516 *augRhs_x, *augRhs_s, *rhs.y_c(), *rhs.y_d(), 517 *sol->x_NonConst(), *sol->s_NonConst(), 518 *sol->y_c_NonConst(), *sol->y_d_NonConst(), check_inertia, numberOfEVals); 519 } 520 if (retval==SYMSOLVER_FATAL_ERROR) return false; 521 if (retval==SYMSOLVER_SINGULAR && 522 (rhs.y_c()->Dim()+rhs.y_d()->Dim() > 0) ) { 523 524 // Get new perturbation factors from the perturbation 525 // handlers for the singular case 526 bool pert_return = perturbHandler_->PerturbForSingularity(delta_x, delta_s, 527 delta_c, delta_d); 528 if (!pert_return) { 529 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 530 "PerturbForSingularity can't be done\n"); 531 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 532 return false; 533 } 534 } 535 else if (retval==SYMSOLVER_WRONG_INERTIA && 536 augSysSolver_->NumberOfNegEVals() < numberOfEVals) { 537 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 538 "Number of negative eigenvalues too small!\n"); 539 // If the number of negative eigenvalues is too small, then 540 // we first try to remedy this by asking for better quality 541 // solution (e.g. increasing pivot tolerance), and if that 542 // doesn't help, we assume that the system is singular 543 bool assume_singular = true; 544 if (!augsys_improved_) { 545 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 546 "Asking augmented system solver to improve quality of its solutions.\n"); 547 augsys_improved_ = augSysSolver_->IncreaseQuality(); 548 if (augsys_improved_) { 549 IpData().Append_info_string("q"); 550 assume_singular = false; 551 } 552 else { 553 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 554 "Quality could not be improved\n"); 555 } 556 } 557 if (assume_singular) { 558 bool pert_return = 559 perturbHandler_->PerturbForSingularity(delta_x, delta_s, 560 delta_c, delta_d); 561 if (!pert_return) { 562 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 563 "PerturbForSingularity can't be done for assume singular.\n"); 564 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 565 return false; 566 } 567 IpData().Append_info_string("a"); 568 } 569 } 570 else if (retval==SYMSOLVER_WRONG_INERTIA || 571 retval==SYMSOLVER_SINGULAR) { 572 // Get new perturbation factors from the perturbation 573 // handlers for the case of wrong inertia 574 bool pert_return = perturbHandler_->PerturbForWrongInertia(delta_x, delta_s, 575 delta_c, delta_d); 576 if (!pert_return) { 577 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 578 "PerturbForWrongInertia can't be done for wrong interia or singular.\n"); 579 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 580 return false; 581 } 582 } 583 else if (neg_curv_test_tol_ > 0.) { 584 DBG_ASSERT(augSysSolver_->ProvidesInertia()); 585 // we now check if the inertia is possible wrong 586 Index neg_values = augSysSolver_->NumberOfNegEVals(); 587 if (neg_values != numberOfEVals) { 588 // check if we have a direction of sufficient positive curvature 589 SmartPtr<Vector> x_tmp = sol->x()->MakeNew(); 590 W.MultVector(1., *sol->x(), 0., *x_tmp); 591 Number xWx = x_tmp->Dot(*sol->x()); 592 x_tmp->Copy(*sol->x()); 593 x_tmp->ElementWiseMultiply(sigma_x); 594 xWx += x_tmp->Dot(*sol->x()); 595 SmartPtr<Vector> s_tmp = sol->s()->MakeNewCopy(); 596 s_tmp->ElementWiseMultiply(sigma_s); 597 xWx += s_tmp->Dot(*sol->s()); 598 if (neg_curv_test_reg_) { 599 x_tmp->Copy(*sol->x()); 600 x_tmp->Scal(delta_x); 601 xWx += x_tmp->Dot(*sol->x()); 602 603 s_tmp->Copy(*sol->s()); 604 s_tmp->Scal(delta_s); 605 xWx += s_tmp->Dot(*sol->s()); 606 } 607 Number xs_nrmsq = pow(sol->x()->Nrm2(),2) + pow(sol->s()->Nrm2(),2); 608 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 609 "In inertia heuristic: xWx = %e xx = %e\n", 610 xWx, xs_nrmsq); 611 if (xWx < neg_curv_test_tol_*xs_nrmsq) { 612 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 613 " -> Redo with modified matrix.\n"); 614 bool pert_return = perturbHandler_->PerturbForWrongInertia(delta_x, delta_s, 615 delta_c, delta_d); 616 if (!pert_return) { 617 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 618 "PerturbForWrongInertia can't be done for inertia heuristic.\n"); 619 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 620 return false; 621 } 622 retval = SYMSOLVER_WRONG_INERTIA; 623 } 624 } 625 } 626 } // while (retval!=SYMSOLVER_SUCCESS && !fail) { 627 628 // Some output 629 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 630 "Number of trial factorizations performed: %d\n", 631 count); 632 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 633 "Perturbation parameters: delta_x=%e delta_s=%e\n delta_c=%e delta_d=%e\n", 634 delta_x, delta_s, delta_c, delta_d); 635 // Set the perturbation values in the Data object 636 IpData().setPDPert(delta_x, delta_s, delta_c, delta_d); 637 } 638 639 // Compute the remaining sol Vectors 640 Px_L.SinvBlrmZMTdBr(-1., slack_x_L, *rhs.z_L(), z_L, *sol->x(), *sol->z_L_NonConst()); 641 Px_U.SinvBlrmZMTdBr(1., slack_x_U, *rhs.z_U(), z_U, *sol->x(), *sol->z_U_NonConst()); 642 Pd_L.SinvBlrmZMTdBr(-1., slack_s_L, *rhs.v_L(), v_L, *sol->s(), *sol->v_L_NonConst()); 643 Pd_U.SinvBlrmZMTdBr(1., slack_s_U, *rhs.v_U(), v_U, *sol->s(), *sol->v_U_NonConst()); 644 645 // Finally let's assemble the res result vectors 646 res.AddOneVector(alpha, *sol, beta); 647 648 IpData().TimingStats().PDSystemSolverSolveOnce().End(); 649 650 return true; 651 } 652 ComputeResiduals(const SymMatrix & W,const Matrix & J_c,const Matrix & J_d,const Matrix & Px_L,const Matrix & Px_U,const Matrix & Pd_L,const Matrix & Pd_U,const Vector & z_L,const Vector & z_U,const Vector & v_L,const Vector & v_U,const Vector & slack_x_L,const Vector & slack_x_U,const Vector & slack_s_L,const Vector & slack_s_U,const Vector & sigma_x,const Vector & sigma_s,Number alpha,Number beta,const IteratesVector & rhs,const IteratesVector & res,IteratesVector & resid)653 void PDFullSpaceSolver::ComputeResiduals( 654 const SymMatrix& W, 655 const Matrix& J_c, 656 const Matrix& J_d, 657 const Matrix& Px_L, 658 const Matrix& Px_U, 659 const Matrix& Pd_L, 660 const Matrix& Pd_U, 661 const Vector& z_L, 662 const Vector& z_U, 663 const Vector& v_L, 664 const Vector& v_U, 665 const Vector& slack_x_L, 666 const Vector& slack_x_U, 667 const Vector& slack_s_L, 668 const Vector& slack_s_U, 669 const Vector& sigma_x, 670 const Vector& sigma_s, 671 Number alpha, 672 Number beta, 673 const IteratesVector& rhs, 674 const IteratesVector& res, 675 IteratesVector& resid) 676 { 677 DBG_START_METH("PDFullSpaceSolver::ComputeResiduals", dbg_verbosity); 678 679 DBG_PRINT_VECTOR(2, "res", res); 680 IpData().TimingStats().ComputeResiduals().Start(); 681 682 // Get the current sizes of the perturbation factors 683 Number delta_x; 684 Number delta_s; 685 Number delta_c; 686 Number delta_d; 687 perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d); 688 689 SmartPtr<Vector> tmp; 690 691 // x 692 W.MultVector(1., *res.x(), 0., *resid.x_NonConst()); 693 J_c.TransMultVector(1., *res.y_c(), 1., *resid.x_NonConst()); 694 J_d.TransMultVector(1., *res.y_d(), 1., *resid.x_NonConst()); 695 Px_L.MultVector(-1., *res.z_L(), 1., *resid.x_NonConst()); 696 Px_U.MultVector(1., *res.z_U(), 1., *resid.x_NonConst()); 697 resid.x_NonConst()->AddTwoVectors(delta_x, *res.x(), -1., *rhs.x(), 1.); 698 699 // s 700 Pd_U.MultVector(1., *res.v_U(), 0., *resid.s_NonConst()); 701 Pd_L.MultVector(-1., *res.v_L(), 1., *resid.s_NonConst()); 702 resid.s_NonConst()->AddTwoVectors(-1., *res.y_d(), -1., *rhs.s(), 1.); 703 if (delta_s!=0.) { 704 resid.s_NonConst()->Axpy(delta_s, *res.s()); 705 } 706 707 // c 708 J_c.MultVector(1., *res.x(), 0., *resid.y_c_NonConst()); 709 resid.y_c_NonConst()->AddTwoVectors(-delta_c, *res.y_c(), -1., *rhs.y_c(), 1.); 710 711 // d 712 J_d.MultVector(1., *res.x(), 0., *resid.y_d_NonConst()); 713 resid.y_d_NonConst()->AddTwoVectors(-1., *res.s(), -1., *rhs.y_d(), 1.); 714 if (delta_d!=0.) { 715 resid.y_d_NonConst()->Axpy(-delta_d, *res.y_d()); 716 } 717 718 // zL 719 resid.z_L_NonConst()->Copy(*res.z_L()); 720 resid.z_L_NonConst()->ElementWiseMultiply(slack_x_L); 721 tmp = z_L.MakeNew(); 722 Px_L.TransMultVector(1., *res.x(), 0., *tmp); 723 tmp->ElementWiseMultiply(z_L); 724 resid.z_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.z_L(), 1.); 725 726 // zU 727 resid.z_U_NonConst()->Copy(*res.z_U()); 728 resid.z_U_NonConst()->ElementWiseMultiply(slack_x_U); 729 tmp = z_U.MakeNew(); 730 Px_U.TransMultVector(1., *res.x(), 0., *tmp); 731 tmp->ElementWiseMultiply(z_U); 732 resid.z_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.z_U(), 1.); 733 734 // vL 735 resid.v_L_NonConst()->Copy(*res.v_L()); 736 resid.v_L_NonConst()->ElementWiseMultiply(slack_s_L); 737 tmp = v_L.MakeNew(); 738 Pd_L.TransMultVector(1., *res.s(), 0., *tmp); 739 tmp->ElementWiseMultiply(v_L); 740 resid.v_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.v_L(), 1.); 741 742 // vU 743 resid.v_U_NonConst()->Copy(*res.v_U()); 744 resid.v_U_NonConst()->ElementWiseMultiply(slack_s_U); 745 tmp = v_U.MakeNew(); 746 Pd_U.TransMultVector(1., *res.s(), 0., *tmp); 747 tmp->ElementWiseMultiply(v_U); 748 resid.v_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.v_U(), 1.); 749 750 DBG_PRINT_VECTOR(2, "resid", resid); 751 752 if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) { 753 resid.Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "resid"); 754 } 755 756 if (Jnlst().ProduceOutput(J_MOREDETAILED, J_LINEAR_ALGEBRA)) { 757 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 758 "max-norm resid_x %e\n", resid.x()->Amax()); 759 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 760 "max-norm resid_s %e\n", resid.s()->Amax()); 761 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 762 "max-norm resid_c %e\n", resid.y_c()->Amax()); 763 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 764 "max-norm resid_d %e\n", resid.y_d()->Amax()); 765 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 766 "max-norm resid_zL %e\n", resid.z_L()->Amax()); 767 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 768 "max-norm resid_zU %e\n", resid.z_U()->Amax()); 769 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 770 "max-norm resid_vL %e\n", resid.v_L()->Amax()); 771 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 772 "max-norm resid_vU %e\n", resid.v_U()->Amax()); 773 } 774 IpData().TimingStats().ComputeResiduals().End(); 775 } 776 ComputeResidualRatio(const IteratesVector & rhs,const IteratesVector & res,const IteratesVector & resid)777 Number PDFullSpaceSolver::ComputeResidualRatio(const IteratesVector& rhs, 778 const IteratesVector& res, 779 const IteratesVector& resid) 780 { 781 DBG_START_METH("PDFullSpaceSolver::ComputeResidualRatio", dbg_verbosity); 782 783 Number nrm_rhs = rhs.Amax(); 784 Number nrm_res = res.Amax(); 785 Number nrm_resid = resid.Amax(); 786 Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA, 787 "nrm_rhs = %8.2e nrm_sol = %8.2e nrm_resid = %8.2e\n", 788 nrm_rhs, nrm_res, nrm_resid); 789 790 if (nrm_rhs+nrm_res == 0.) { 791 return nrm_resid; // this should be zero 792 } 793 else { 794 // ToDo: determine how to include norm of matrix, and what 795 // safeguard to use against incredibly large solution vectors 796 Number max_cond = 1e6; 797 return nrm_resid/(Min(nrm_res, max_cond*nrm_rhs)+nrm_rhs); 798 } 799 } 800 801 } // namespace Ipopt 802