1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the class library */ 4 /* SoPlex --- the Sequential object-oriented simPlex. */ 5 /* */ 6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 #include <iostream> 26 27 #include "soplex/array.h" 28 #include "soplex/dataarray.h" 29 #include "soplex/sorter.h" 30 #include "soplex/spxout.h" 31 #include <sstream> 32 #include <iostream> 33 #include <fstream> 34 #include <memory> 35 36 37 //rows 38 #define SOPLEX_FREE_LHS_RHS 1 39 #define SOPLEX_FREE_CONSTRAINT 1 40 #define SOPLEX_EMPTY_CONSTRAINT 1 41 #define SOPLEX_ROW_SINGLETON 1 42 #define SOPLEX_AGGREGATE_VARS 1 43 #define SOPLEX_FORCE_CONSTRAINT 1 44 //cols 45 #define SOPLEX_FREE_BOUNDS 1 46 #define SOPLEX_COMPTY_COLUMN 1 47 #define SOPLEX_FIX_VARIABLE 1 48 #define SOPLEX_FREE_ZERO_OBJ_VARIABLE 1 49 #define SOPLEX_ZERO_OBJ_COL_SINGLETON 1 50 #define SOPLEX_DOUBLETON_EQUATION 1 51 #define SOPLEX_FREE_COL_SINGLETON 1 52 //dual 53 #define SOPLEX_DOMINATED_COLUMN 1 54 #define SOPLEX_WEAKLY_DOMINATED_COLUMN 1 55 #define SOPLEX_MULTI_AGGREGATE 1 56 //other 57 #define SOPLEX_TRIVIAL_HEURISTICS 1 58 #define SOPLEX_PSEUDOOBJ 1 59 60 61 #define SOPLEX_EXTREMES 1 62 #define SOPLEX_ROWS_SPXMAINSM 1 63 #define SOPLEX_COLS_SPXMAINSM 1 64 #define SOPLEX_DUAL_SPXMAINSM 1 65 ///@todo check: with this simplification step, the unsimplified basis seems to be slightly suboptimal for some instances 66 #define SOPLEX_DUPLICATE_ROWS 1 67 #define SOPLEX_DUPLICATE_COLS 1 68 69 70 #ifndef NDEBUG 71 #define SOPLEX_CHECK_BASIS_DIM 72 #endif // NDEBUG 73 74 namespace soplex 75 { 76 77 template <class R> 78 bool SPxMainSM<R>::PostStep::checkBasisDim(DataArray<typename SPxSolverBase<R>::VarStatus> rows, 79 DataArray<typename SPxSolverBase<R>::VarStatus> cols) const 80 { 81 int numBasis = 0; 82 83 for(int rs = 0; rs < nRows; ++rs) 84 { 85 if(rows[rs] == SPxSolverBase<R>::BASIC) 86 numBasis++; 87 } 88 89 for(int cs = 0; cs < nCols; ++cs) 90 { 91 if(cols[cs] == SPxSolverBase<R>::BASIC) 92 numBasis++; 93 } 94 95 assert(numBasis == nRows); 96 return numBasis == nRows; 97 } 98 99 template <class R> 100 void SPxMainSM<R>::RowObjPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 101 VectorBase<R>&, 102 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 103 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 104 { 105 s[m_i] = s[m_i] - x[m_j]; 106 107 assert(rStatus[m_i] != SPxSolverBase<R>::UNDEFINED); 108 assert(cStatus[m_j] != SPxSolverBase<R>::UNDEFINED); 109 assert(rStatus[m_i] != SPxSolverBase<R>::BASIC || cStatus[m_j] != SPxSolverBase<R>::BASIC); 110 111 SPxOut::debug(this, "RowObjPS: removing slack column {} ({}) for row {} ({}).\n", m_j, cStatus[m_j], 112 m_i, rStatus[m_i]); 113 114 if(rStatus[m_i] != SPxSolverBase<R>::BASIC) 115 { 116 switch(cStatus[m_j]) 117 { 118 case SPxSolverBase<R>::ON_UPPER: 119 rStatus[m_i] = SPxSolverBase<R>::ON_LOWER; 120 break; 121 122 case SPxSolverBase<R>::ON_LOWER: 123 rStatus[m_i] = SPxSolverBase<R>::ON_UPPER; 124 break; 125 126 default: 127 rStatus[m_i] = cStatus[m_j]; 128 } 129 130 // otherwise checkBasisDim() may fail 131 cStatus[m_j] = SPxSolverBase<R>::ZERO; 132 } 133 134 #ifdef SOPLEX_CHECK_BASIS_DIM 135 136 if(!this->checkBasisDim(rStatus, cStatus)) 137 { 138 assert(false); 139 throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step."); 140 } 141 142 #endif 143 } 144 145 template <class R> 146 void SPxMainSM<R>::FreeConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 147 VectorBase<R>&, 148 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 149 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 150 { 151 // correcting the change of idx by deletion of the row: 152 if(m_i != m_old_i) 153 { 154 s[m_old_i] = s[m_i]; 155 y[m_old_i] = y[m_i]; 156 rStatus[m_old_i] = rStatus[m_i]; 157 } 158 159 // primal: 160 R slack = 0.0; 161 162 for(int k = 0; k < m_row.size(); ++k) 163 slack += m_row.value(k) * x[m_row.index(k)]; 164 165 s[m_i] = slack; 166 167 // dual: 168 y[m_i] = m_row_obj; 169 170 // basis: 171 rStatus[m_i] = SPxSolverBase<R>::BASIC; 172 173 #ifdef SOPLEX_CHECK_BASIS_DIM 174 175 if(!this->checkBasisDim(rStatus, cStatus)) 176 { 177 throw SPxInternalCodeException("XMAISM15 Dimension doesn't match after this step."); 178 } 179 180 #endif 181 } 182 183 template <class R> 184 void SPxMainSM<R>::EmptyConstraintPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s, 185 VectorBase<R>&, 186 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 187 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 188 { 189 // correcting the change of idx by deletion of the row: 190 if(m_i != m_old_i) 191 { 192 s[m_old_i] = s[m_i]; 193 y[m_old_i] = y[m_i]; 194 rStatus[m_old_i] = rStatus[m_i]; 195 } 196 197 // primal: 198 s[m_i] = 0.0; 199 200 // dual: 201 y[m_i] = m_row_obj; 202 203 // basis: 204 rStatus[m_i] = SPxSolverBase<R>::BASIC; 205 206 #ifdef SOPLEX_CHECK_BASIS_DIM 207 208 if(!this->checkBasisDim(rStatus, cStatus)) 209 { 210 throw SPxInternalCodeException("XMAISM16 Dimension doesn't match after this step."); 211 } 212 213 #endif 214 } 215 216 template <class R> 217 void SPxMainSM<R>::RowSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 218 VectorBase<R>& r, 219 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 220 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 221 { 222 // correcting the change of idx by deletion of the row: 223 if(m_i != m_old_i) 224 { 225 y[m_old_i] = y[m_i]; 226 s[m_old_i] = s[m_i]; 227 rStatus[m_old_i] = rStatus[m_i]; 228 } 229 230 R aij = m_col[m_i]; 231 assert(aij != 0.0); 232 233 // primal: 234 s[m_i] = aij * x[m_j]; 235 236 // dual & basis: 237 R val = m_obj; 238 239 for(int k = 0; k < m_col.size(); ++k) 240 { 241 if(m_col.index(k) != m_i) 242 val -= m_col.value(k) * y[m_col.index(k)]; 243 } 244 245 R newLo = (aij > 0) ? m_lhs / aij : m_rhs / aij; // implicit lhs 246 R newUp = (aij > 0) ? m_rhs / aij : m_lhs / aij; // implicit rhs 247 248 switch(cStatus[m_j]) 249 { 250 case SPxSolverBase<R>::FIXED: 251 if(newLo <= m_oldLo && newUp >= m_oldUp) 252 { 253 // this row is totally redundant, has not changed bound of xj 254 rStatus[m_i] = SPxSolverBase<R>::BASIC; 255 y[m_i] = m_row_obj; 256 } 257 else if(EQrel(newLo, newUp, this->feastol())) 258 { 259 // row is in the type aij * xj = b 260 assert(EQrel(newLo, x[m_j], this->feastol())); 261 262 if(EQrel(m_oldLo, m_oldUp, this->feastol())) 263 { 264 // xj has been fixed in other row 265 rStatus[m_i] = SPxSolverBase<R>::BASIC; 266 y[m_i] = m_row_obj; 267 } 268 else if((EQrel(m_oldLo, x[m_j], this->feastol()) && r[m_j] <= -this->feastol()) 269 || (EQrel(m_oldUp, x[m_j], this->feastol()) && r[m_j] >= this->feastol()) 270 || (!EQrel(m_oldLo, x[m_j], this->feastol()) && !(EQrel(m_oldUp, x[m_j], this->feastol())))) 271 { 272 // if x_j on lower but reduced cost is negative, or x_j on upper but reduced cost is positive, or x_j not on bound: basic 273 rStatus[m_i] = (EQrel(m_lhs, x[m_j] * aij, 274 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 275 cStatus[m_j] = SPxSolverBase<R>::BASIC; 276 y[m_i] = val / aij; 277 r[m_j] = 0.0; 278 } 279 else 280 { 281 // set x_j on one of the bound 282 assert(EQrel(m_oldLo, x[m_j], this->feastol()) || EQrel(m_oldUp, x[m_j], this->feastol())); 283 284 cStatus[m_j] = EQrel(m_oldLo, x[m_j], 285 this->feastol()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 286 rStatus[m_i] = SPxSolverBase<R>::BASIC; 287 y[m_i] = m_row_obj; 288 r[m_j] = val; 289 } 290 } 291 else if(EQrel(newLo, m_oldUp, this->feastol())) 292 { 293 // row is in the type xj >= b/aij, try to set xj on upper 294 if(r[m_j] >= this->feastol()) 295 { 296 // the reduced cost is positive, xj should in the basic 297 assert(EQrel(m_rhs / aij, x[m_j], this->feastol()) || EQrel(m_lhs / aij, x[m_j], this->feastol())); 298 299 rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j], 300 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 301 cStatus[m_j] = SPxSolverBase<R>::BASIC; 302 y[m_i] = val / aij; 303 r[m_j] = 0.0; 304 } 305 else 306 { 307 assert(EQrel(m_oldUp, x[m_j], this->feastol())); 308 309 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 310 rStatus[m_i] = SPxSolverBase<R>::BASIC; 311 y[m_i] = m_row_obj; 312 r[m_j] = val; 313 } 314 } 315 else if(EQrel(newUp, m_oldLo, this->feastol())) 316 { 317 // row is in the type xj <= b/aij, try to set xj on lower 318 if(r[m_j] <= -this->feastol()) 319 { 320 // the reduced cost is negative, xj should in the basic 321 assert(EQrel(m_rhs / aij, x[m_j], this->feastol()) || EQrel(m_lhs / aij, x[m_j], this->feastol())); 322 323 rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j], 324 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 325 cStatus[m_j] = SPxSolverBase<R>::BASIC; 326 y[m_i] = val / aij; 327 r[m_j] = 0.0; 328 } 329 else 330 { 331 assert(EQrel(m_oldLo, x[m_j], this->feastol())); 332 333 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 334 rStatus[m_i] = SPxSolverBase<R>::BASIC; 335 y[m_i] = m_row_obj; 336 r[m_j] = val; 337 } 338 } 339 else 340 { 341 // the variable is set to FIXED by other constraints, i.e., this singleton row is redundant 342 rStatus[m_i] = SPxSolverBase<R>::BASIC; 343 y[m_i] = m_row_obj; 344 } 345 346 break; 347 348 case SPxSolverBase<R>::BASIC: 349 rStatus[m_i] = SPxSolverBase<R>::BASIC; 350 y[m_i] = m_row_obj; 351 r[m_j] = 0.0; 352 break; 353 354 case SPxSolverBase<R>::ON_LOWER: 355 if(EQrel(m_oldLo, x[m_j], this->feastol())) // xj may stay on lower 356 { 357 rStatus[m_i] = SPxSolverBase<R>::BASIC; 358 y[m_i] = m_row_obj; 359 r[m_j] = val; 360 } 361 else // if reduced costs are negative or old lower bound not equal to xj, we need to change xj into the basis 362 { 363 assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->feastol()) 364 || EQrel(m_lhs / aij, x[m_j], this->feastol())); 365 366 cStatus[m_j] = SPxSolverBase<R>::BASIC; 367 rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j], 368 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 369 y[m_i] = val / aij; 370 r[m_j] = 0.0; 371 } 372 373 break; 374 375 case SPxSolverBase<R>::ON_UPPER: 376 if(EQrel(m_oldUp, x[m_j], this->feastol())) // xj may stay on upper 377 { 378 rStatus[m_i] = SPxSolverBase<R>::BASIC; 379 y[m_i] = m_row_obj; 380 r[m_j] = val; 381 } 382 else // if reduced costs are positive or old upper bound not equal to xj, we need to change xj into the basis 383 { 384 assert(!isOptimal || EQrel(m_rhs / aij, x[m_j], this->feastol()) 385 || EQrel(m_lhs / aij, x[m_j], this->feastol())); 386 387 cStatus[m_j] = SPxSolverBase<R>::BASIC; 388 rStatus[m_i] = (EQrel(m_lhs / aij, x[m_j], 389 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 390 y[m_i] = val / aij; 391 r[m_j] = 0.0; 392 } 393 394 break; 395 396 case SPxSolverBase<R>::ZERO: 397 rStatus[m_i] = SPxSolverBase<R>::BASIC; 398 y[m_i] = m_row_obj; 399 r[m_j] = val; 400 break; 401 402 default: 403 break; 404 } 405 406 #ifdef SOPLEX_CHECK_BASIS_DIM 407 408 if(!this->checkBasisDim(rStatus, cStatus)) 409 { 410 throw SPxInternalCodeException("XMAISM17 Dimension doesn't match after this step."); 411 } 412 413 #endif 414 } 415 416 template <class R> 417 void SPxMainSM<R>::ForceConstraintPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 418 VectorBase<R>& r, 419 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 420 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 421 { 422 // correcting the change of idx by deletion of the row: 423 if(m_i != m_old_i) 424 { 425 s[m_old_i] = s[m_i]; 426 y[m_old_i] = y[m_i]; 427 rStatus[m_old_i] = rStatus[m_i]; 428 } 429 430 // primal: 431 s[m_i] = m_lRhs; 432 433 // basis: 434 int cBasisCandidate = -1; 435 R maxViolation = -1.0; 436 int bas_k = -1; 437 438 for(int k = 0; k < m_row.size(); ++k) 439 { 440 int cIdx = m_row.index(k); 441 R aij = m_row.value(k); 442 R oldLo = m_oldLowers[k]; 443 R oldUp = m_oldUppers[k]; 444 445 switch(cStatus[cIdx]) 446 { 447 case SPxSolverBase<R>::FIXED: 448 if(m_fixed[k]) 449 { 450 assert(EQrel(oldLo, x[cIdx], this->feastol()) || EQrel(oldUp, x[cIdx], this->feastol())); 451 452 R violation = spxAbs(r[cIdx] / aij); 453 454 cStatus[cIdx] = EQrel(oldLo, x[cIdx], 455 this->feastol()) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 456 457 if(violation > maxViolation && ((EQrel(oldLo, x[cIdx], this->feastol()) 458 && r[cIdx] < -this->feastol()) 459 || (EQrel(oldUp, x[cIdx], this->feastol()) && r[cIdx] > this->feastol()))) 460 { 461 maxViolation = violation; 462 cBasisCandidate = cIdx; 463 bas_k = k; 464 } 465 } // do nothing, if the old bounds are equal, i.e. variable has been not fixed in this row 466 467 break; 468 469 case SPxSolverBase<R>::ON_LOWER: 470 case SPxSolverBase<R>::ON_UPPER: 471 case SPxSolverBase<R>::BASIC: 472 break; 473 474 default: 475 break; 476 } 477 } 478 479 // dual and basis : 480 if(cBasisCandidate >= 0) // one of the variable in the row should in the basis 481 { 482 assert(EQrel(m_lRhs, m_rhs, this->feastol()) || EQrel(m_lRhs, m_lhs, this->feastol())); 483 assert(bas_k >= 0); 484 assert(cBasisCandidate == m_row.index(bas_k)); 485 486 cStatus[cBasisCandidate] = SPxSolverBase<R>::BASIC; 487 rStatus[m_i] = (EQrel(m_lRhs, m_lhs, 488 this->feastol())) ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER; 489 490 R aij = m_row.value(bas_k); 491 R multiplier = r[cBasisCandidate] / aij; 492 r[cBasisCandidate] = 0.0; 493 494 for(int k = 0; k < m_row.size(); ++k) // update the reduced cost 495 { 496 if(k == bas_k) 497 { 498 continue; 499 } 500 501 r[m_row.index(k)] -= m_row.value(k) * multiplier; 502 } 503 504 // compute the value of new dual variable (because we have a new row) 505 R val = m_objs[bas_k]; 506 DSVectorBase<R> basis_col = m_cols[bas_k]; 507 508 for(int k = 0; k < basis_col.size(); ++k) 509 { 510 if(basis_col.index(k) != m_i) 511 val -= basis_col.value(k) * y[basis_col.index(k)]; 512 } 513 514 y[m_i] = val / aij; 515 } 516 else // slack in the basis 517 { 518 rStatus[m_i] = SPxSolverBase<R>::BASIC; 519 y[m_i] = m_rowobj; 520 } 521 522 #ifdef SOPLEX_CHECK_BASIS_DIM 523 524 if(!this->checkBasisDim(rStatus, cStatus)) 525 { 526 throw SPxInternalCodeException("XMAISM18 Dimension doesn't match after this step."); 527 } 528 529 #endif 530 } 531 532 template <class R> 533 void SPxMainSM<R>::FixVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 534 VectorBase<R>& r, 535 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 536 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 537 { 538 // update the index mapping; if m_correctIdx is false, we assume that this has happened already 539 if(m_correctIdx) 540 { 541 x[m_old_j] = x[m_j]; 542 r[m_old_j] = r[m_j]; 543 cStatus[m_old_j] = cStatus[m_j]; 544 } 545 546 // primal: 547 x[m_j] = m_val; 548 549 for(int k = 0; k < m_col.size(); ++k) 550 s[m_col.index(k)] += m_col.value(k) * x[m_j]; 551 552 // dual: 553 R val = m_obj; 554 555 for(int k = 0; k < m_col.size(); ++k) 556 val -= m_col.value(k) * y[m_col.index(k)]; 557 558 r[m_j] = val; 559 560 // basis: 561 if(m_lower == m_upper) 562 { 563 assert(EQrel(m_lower, m_val, this->epsilon())); 564 565 cStatus[m_j] = SPxSolverBase<R>::FIXED; 566 } 567 else 568 { 569 assert(EQrel(m_val, m_lower, this->epsilon()) || EQrel(m_val, m_upper, this->epsilon()) 570 || m_val == 0.0); 571 572 cStatus[m_j] = EQrel(m_val, m_lower, this->epsilon()) ? SPxSolverBase<R>::ON_LOWER : (EQrel(m_val, 573 m_upper, this->epsilon()) ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ZERO); 574 } 575 576 #ifdef SOPLEX_CHECK_BASIS_DIM 577 578 if(m_correctIdx) 579 { 580 if(!this->checkBasisDim(rStatus, cStatus)) 581 { 582 throw SPxInternalCodeException("XMAISM19 Dimension doesn't match after this step."); 583 } 584 } 585 586 #endif 587 } 588 589 template <class R> 590 void SPxMainSM<R>::FixBoundsPS::execute(VectorBase<R>&, VectorBase<R>&, VectorBase<R>&, 591 VectorBase<R>&, 592 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 593 DataArray<typename SPxSolverBase<R>::VarStatus>&, bool isOptimal) const 594 { 595 // basis: 596 cStatus[m_j] = m_status; 597 } 598 599 template <class R> 600 void SPxMainSM<R>::FreeZeroObjVariablePS::execute(VectorBase<R>& x, VectorBase<R>& y, 601 VectorBase<R>& s, VectorBase<R>& r, 602 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 603 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 604 { 605 // correcting the change of idx by deletion of the column and corresponding rows: 606 if(m_j != m_old_j) 607 { 608 x[m_old_j] = x[m_j]; 609 r[m_old_j] = r[m_j]; 610 cStatus[m_old_j] = cStatus[m_j]; 611 } 612 613 int rIdx = m_old_i - m_col.size() + 1; 614 615 for(int k = 0; k < m_col.size(); ++k) 616 { 617 int rIdx_new = m_col.index(k); 618 s[rIdx] = s[rIdx_new]; 619 y[rIdx] = y[rIdx_new]; 620 rStatus[rIdx] = rStatus[rIdx_new]; 621 rIdx++; 622 } 623 624 // primal: 625 int domIdx = -1; 626 DSVectorBase<R> slack(m_col.size()); 627 628 if(m_loFree) 629 { 630 R minRowUp = R(infinity); 631 632 for(int k = 0; k < m_rows.size(); ++k) 633 { 634 R val = 0.0; 635 const SVectorBase<R>& row = m_rows[k]; 636 637 for(int l = 0; l < row.size(); ++l) 638 { 639 if(row.index(l) != m_j) 640 val += row.value(l) * x[row.index(l)]; 641 } 642 643 R scale = maxAbs(m_lRhs[k], val); 644 645 if(scale < 1.0) 646 scale = 1.0; 647 648 R z = (m_lRhs[k] / scale) - (val / scale); 649 650 if(isZero(z, this->epsilon())) 651 z = 0.0; 652 653 R up = z * scale / row[m_j]; 654 slack.add(k, val); 655 656 if(up < minRowUp) 657 { 658 minRowUp = up; 659 domIdx = k; 660 } 661 } 662 663 if(m_bnd < minRowUp) 664 { 665 x[m_j] = m_bnd; 666 domIdx = -1; 667 } 668 else 669 x[m_j] = minRowUp; 670 } 671 else 672 { 673 R maxRowLo = R(-infinity); 674 675 for(int k = 0; k < m_rows.size(); ++k) 676 { 677 R val = 0.0; 678 const SVectorBase<R>& row = m_rows[k]; 679 680 for(int l = 0; l < row.size(); ++l) 681 { 682 if(row.index(l) != m_j) 683 val += row.value(l) * x[row.index(l)]; 684 } 685 686 R scale = maxAbs(m_lRhs[k], val); 687 688 if(scale < 1.0) 689 scale = 1.0; 690 691 R z = (m_lRhs[k] / scale) - (val / scale); 692 693 if(isZero(z, this->epsilon())) 694 z = 0.0; 695 696 R lo = z * scale / row[m_j]; 697 slack.add(k, val); 698 699 if(lo > maxRowLo) 700 { 701 maxRowLo = lo; 702 domIdx = k; 703 } 704 } 705 706 if(m_bnd > maxRowLo) 707 { 708 x[m_j] = m_bnd; 709 domIdx = -1; 710 } 711 else 712 x[m_j] = maxRowLo; 713 } 714 715 for(int k = 0; k < m_col.size(); ++k) 716 s[m_col.index(k)] = slack[k] + m_col.value(k) * x[m_j]; 717 718 // dual: 719 r[m_j] = 0.0; 720 721 for(int k = 0; k < m_col.size(); ++k) 722 { 723 int idx = m_col.index(k); 724 y[idx] = m_rowObj[idx]; 725 } 726 727 // basis: 728 for(int k = 0; k < m_col.size(); ++k) 729 { 730 if(k != domIdx) 731 rStatus[m_col.index(k)] = SPxSolverBase<R>::BASIC; 732 733 else 734 { 735 cStatus[m_j] = SPxSolverBase<R>::BASIC; 736 737 if(m_loFree) 738 rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_UPPER : 739 SPxSolverBase<R>::ON_LOWER; 740 else 741 rStatus[m_col.index(k)] = (m_col.value(k) > 0) ? SPxSolverBase<R>::ON_LOWER : 742 SPxSolverBase<R>::ON_UPPER; 743 } 744 } 745 746 if(domIdx == -1) 747 { 748 if(m_loFree) 749 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 750 else 751 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 752 } 753 754 #ifdef SOPLEX_CHECK_BASIS_DIM 755 756 if(!this->checkBasisDim(rStatus, cStatus)) 757 { 758 throw SPxInternalCodeException("XMAISM20 Dimension doesn't match after this step."); 759 } 760 761 #endif 762 } 763 764 template <class R> 765 void SPxMainSM<R>::ZeroObjColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, 766 VectorBase<R>& s, VectorBase<R>& r, 767 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 768 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 769 { 770 // correcting the change of idx by deletion of the column and corresponding rows: 771 if(m_j != m_old_j) 772 { 773 x[m_old_j] = x[m_j]; 774 r[m_old_j] = r[m_j]; 775 cStatus[m_old_j] = cStatus[m_j]; 776 } 777 778 // primal & basis: 779 R aij = m_row[m_j]; 780 781 if(isZero(s[m_i], R(1e-6))) 782 s[m_i] = 0.0; 783 else if(s[m_i] >= R(infinity)) 784 // this is a fix for a highly ill conditioned instance that is "solved" in presolving (ilaser0 from MINLP, mittelmann) 785 throw SPxException("Simplifier: infinite activities - aborting unsimplification"); 786 787 R scale1 = maxAbs(m_lhs, s[m_i]); 788 R scale2 = maxAbs(m_rhs, s[m_i]); 789 790 if(scale1 < 1.0) 791 scale1 = 1.0; 792 793 if(scale2 < 1.0) 794 scale2 = 1.0; 795 796 R z1 = (m_lhs / scale1) - (s[m_i] / scale1); 797 R z2 = (m_rhs / scale2) - (s[m_i] / scale2); 798 799 if(isZero(z1, this->epsilon())) 800 z1 = 0.0; 801 802 if(isZero(z2, this->epsilon())) 803 z2 = 0.0; 804 805 R lo = (aij > 0) ? z1 * scale1 / aij : z2 * scale2 / aij; 806 R up = (aij > 0) ? z2 * scale2 / aij : z1 * scale1 / aij; 807 808 if(isZero(lo, this->feastol())) 809 lo = 0.0; 810 811 if(isZero(up, this->feastol())) 812 up = 0.0; 813 814 assert(LErel(lo, up, this->epsilon())); 815 SOPLEX_ASSERT_WARN("WMAISM01", isNotZero(aij, R(1.0 / R(infinity)))); 816 817 if(rStatus[m_i] == SPxSolverBase<R>::ON_LOWER) 818 { 819 if(m_lower <= R(-infinity) && m_upper >= R(infinity)) 820 { 821 x[m_j] = 0.0; 822 cStatus[m_j] = SPxSolverBase<R>::ZERO; 823 } 824 else if(m_lower == m_upper) 825 { 826 x[m_j] = m_lower; 827 cStatus[m_j] = SPxSolverBase<R>::FIXED; 828 } 829 else if(aij > 0) 830 { 831 x[m_j] = m_upper; 832 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 833 } 834 else if(aij < 0) 835 { 836 x[m_j] = m_lower; 837 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 838 } 839 else 840 throw SPxInternalCodeException("XMAISM01 This should never happen."); 841 } 842 else if(rStatus[m_i] == SPxSolverBase<R>::ON_UPPER) 843 { 844 if(m_lower <= R(-infinity) && m_upper >= R(infinity)) 845 { 846 x[m_j] = 0.0; 847 cStatus[m_j] = SPxSolverBase<R>::ZERO; 848 } 849 else if(m_lower == m_upper) 850 { 851 x[m_j] = m_lower; 852 cStatus[m_j] = SPxSolverBase<R>::FIXED; 853 } 854 else if(aij > 0) 855 { 856 x[m_j] = m_lower; 857 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 858 } 859 else if(aij < 0) 860 { 861 x[m_j] = m_upper; 862 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 863 } 864 else 865 throw SPxInternalCodeException("XMAISM02 This should never happen."); 866 } 867 else if(rStatus[m_i] == SPxSolverBase<R>::FIXED) 868 { 869 if(m_lower <= R(-infinity) && m_upper >= R(infinity)) 870 { 871 x[m_j] = 0.0; 872 cStatus[m_j] = SPxSolverBase<R>::ZERO; 873 } 874 else 875 { 876 assert(EQrel(m_lower, m_upper, this->feastol())); 877 878 x[m_j] = (m_lower + m_upper) / 2.0; 879 cStatus[m_j] = SPxSolverBase<R>::FIXED; 880 } 881 } 882 else if(rStatus[m_i] == SPxSolverBase<R>::BASIC) 883 { 884 if(GErel(m_lower, lo, this->feastol()) && m_lower > R(-infinity)) 885 { 886 x[m_j] = m_lower; 887 cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 888 } 889 else if(LErel(m_upper, up, this->feastol()) && m_upper < R(infinity)) 890 { 891 x[m_j] = m_upper; 892 cStatus[m_j] = (m_lower == m_upper) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 893 } 894 else if(lo > R(-infinity)) 895 { 896 // make m_i non-basic and m_j basic 897 x[m_j] = lo; 898 cStatus[m_j] = SPxSolverBase<R>::BASIC; 899 rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_LOWER : SPxSolverBase<R>::ON_UPPER); 900 } 901 else if(up < R(infinity)) 902 { 903 // make m_i non-basic and m_j basic 904 x[m_j] = up; 905 cStatus[m_j] = SPxSolverBase<R>::BASIC; 906 rStatus[m_i] = (aij > 0 ? SPxSolverBase<R>::ON_UPPER : SPxSolverBase<R>::ON_LOWER); 907 } 908 else 909 throw SPxInternalCodeException("XMAISM03 This should never happen."); 910 } 911 else 912 throw SPxInternalCodeException("XMAISM04 This should never happen."); 913 914 s[m_i] += aij * x[m_j]; 915 916 // dual: 917 r[m_j] = -1.0 * aij * y[m_i]; 918 919 assert(!isOptimal || (cStatus[m_j] != SPxSolverBase<R>::BASIC || isZero(r[m_j], this->feastol()))); 920 921 #ifdef SOPLEX_CHECK_BASIS_DIM 922 923 if(!this->checkBasisDim(rStatus, cStatus)) 924 { 925 throw SPxInternalCodeException("XMAISM21 Dimension doesn't match after this step."); 926 } 927 928 #endif 929 } 930 931 template <class R> 932 void SPxMainSM<R>::FreeColSingletonPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 933 VectorBase<R>& r, 934 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 935 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 936 { 937 938 // correcting the change of idx by deletion of the row: 939 if(m_i != m_old_i) 940 { 941 s[m_old_i] = s[m_i]; 942 y[m_old_i] = y[m_i]; 943 rStatus[m_old_i] = rStatus[m_i]; 944 } 945 946 // correcting the change of idx by deletion of the column: 947 if(m_j != m_old_j) 948 { 949 x[m_old_j] = x[m_j]; 950 r[m_old_j] = r[m_j]; 951 cStatus[m_old_j] = cStatus[m_j]; 952 } 953 954 // primal: 955 R val = 0.0; 956 R aij = m_row[m_j]; 957 958 for(int k = 0; k < m_row.size(); ++k) 959 { 960 if(m_row.index(k) != m_j) 961 val += m_row.value(k) * x[m_row.index(k)]; 962 } 963 964 R scale = maxAbs(m_lRhs, val); 965 966 if(scale < 1.0) 967 scale = 1.0; 968 969 R z = (m_lRhs / scale) - (val / scale); 970 971 if(isZero(z, this->epsilon())) 972 z = 0.0; 973 974 x[m_j] = z * scale / aij; 975 s[m_i] = m_lRhs; 976 977 // dual: 978 y[m_i] = m_obj / aij; 979 r[m_j] = 0.0; 980 981 // basis: 982 cStatus[m_j] = SPxSolverBase<R>::BASIC; 983 984 if(m_eqCons) 985 rStatus[m_i] = SPxSolverBase<R>::FIXED; 986 else if(m_onLhs) 987 rStatus[m_i] = SPxSolverBase<R>::ON_LOWER; 988 else 989 rStatus[m_i] = SPxSolverBase<R>::ON_UPPER; 990 991 #ifdef SOPLEX_CHECK_BASIS_DIM 992 993 if(!this->checkBasisDim(rStatus, cStatus)) 994 { 995 throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step."); 996 } 997 998 #endif 999 } 1000 1001 template <class R> 1002 void SPxMainSM<R>::DoubletonEquationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>&, 1003 VectorBase<R>& r, 1004 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1005 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1006 { 1007 // dual: 1008 if((cStatus[m_k] != SPxSolverBase<R>::BASIC) && 1009 ((cStatus[m_k] == SPxSolverBase<R>::ON_LOWER && m_strictLo) || 1010 (cStatus[m_k] == SPxSolverBase<R>::ON_UPPER && m_strictUp) || 1011 (cStatus[m_k] == SPxSolverBase<R>::FIXED && 1012 ((m_maxSense && ((r[m_j] > 0 && m_strictUp) || (r[m_j] < 0 && m_strictLo))) || 1013 (!m_maxSense && ((r[m_j] > 0 && m_strictLo) || (r[m_j] < 0 && m_strictUp))))))) 1014 { 1015 R val = m_kObj; 1016 R aik = m_col[m_i]; 1017 1018 for(int _k = 0; _k < m_col.size(); ++_k) 1019 { 1020 if(m_col.index(_k) != m_i) 1021 val -= m_col.value(_k) * y[m_col.index(_k)]; 1022 } 1023 1024 y[m_i] = val / aik; 1025 r[m_k] = 0.0; 1026 1027 r[m_j] = m_jObj - val * m_aij / aik; 1028 1029 SOPLEX_ASSERT_WARN("WMAISM73", isNotZero(m_aij * aik, this->epsilon())); 1030 1031 // basis: 1032 if(m_jFixed) 1033 cStatus[m_j] = SPxSolverBase<R>::FIXED; 1034 else 1035 { 1036 if(GT(r[m_j], (R) 0, this->epsilon()) || (isZero(r[m_j], this->epsilon()) 1037 && EQ(x[m_j], m_Lo_j, this->epsilon()))) 1038 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 1039 else 1040 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 1041 } 1042 1043 cStatus[m_k] = SPxSolverBase<R>::BASIC; 1044 } 1045 1046 #ifdef SOPLEX_CHECK_BASIS_DIM 1047 1048 if(!this->checkBasisDim(rStatus, cStatus)) 1049 { 1050 throw SPxInternalCodeException("XMAISM23 Dimension doesn't match after this step."); 1051 } 1052 1053 #endif 1054 } 1055 1056 template <class R> 1057 void SPxMainSM<R>::DuplicateRowsPS::execute(VectorBase<R>&, VectorBase<R>& y, VectorBase<R>& s, 1058 VectorBase<R>&, 1059 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1060 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1061 { 1062 // correcting the change of idx by deletion of the duplicated rows: 1063 if(m_isLast) 1064 { 1065 for(int i = m_perm.size() - 1; i >= 0; --i) 1066 { 1067 if(m_perm[i] >= 0) 1068 { 1069 int rIdx_new = m_perm[i]; 1070 int rIdx = i; 1071 s[rIdx] = s[rIdx_new]; 1072 y[rIdx] = y[rIdx_new]; 1073 rStatus[rIdx] = rStatus[rIdx_new]; 1074 } 1075 } 1076 } 1077 1078 // primal: 1079 for(int k = 0; k < m_scale.size(); ++k) 1080 { 1081 if(m_scale.index(k) != m_i) 1082 s[m_scale.index(k)] = s[m_i] / m_scale.value(k); 1083 } 1084 1085 // dual & basis: 1086 bool haveSetBasis = false; 1087 1088 for(int k = 0; k < m_scale.size(); ++k) 1089 { 1090 int i = m_scale.index(k); 1091 1092 if(rStatus[m_i] == SPxSolverBase<R>::BASIC || (haveSetBasis && i != m_i)) 1093 // if the row with tightest lower and upper bound in the basic, every duplicate row should in basic 1094 // or basis status of row m_i has been set, this row should be in basis 1095 { 1096 y[i] = m_rowObj.value(k); 1097 rStatus[i] = SPxSolverBase<R>::BASIC; 1098 continue; 1099 } 1100 1101 SOPLEX_ASSERT_WARN("WMAISM02", isNotZero(m_scale.value(k), this->epsilon())); 1102 1103 if(rStatus[m_i] == SPxSolverBase<R>::FIXED && (i == m_maxLhsIdx || i == m_minRhsIdx)) 1104 { 1105 // this row leads to the tightest lower or upper bound, slack should not be in the basis 1106 y[i] = y[m_i] * m_scale.value(k); 1107 y[m_i] = m_i_rowObj; 1108 1109 if(m_isLhsEqualRhs[k]) 1110 { 1111 rStatus[i] = SPxSolverBase<R>::FIXED; 1112 } 1113 else if(i == m_maxLhsIdx) 1114 { 1115 rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER : 1116 SPxSolverBase<R>::ON_UPPER; 1117 } 1118 else 1119 { 1120 assert(i == m_minRhsIdx); 1121 1122 rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER : 1123 SPxSolverBase<R>::ON_LOWER; 1124 } 1125 1126 if(i != m_i) 1127 rStatus[m_i] = SPxSolverBase<R>::BASIC; 1128 1129 haveSetBasis = true; 1130 } 1131 else if(i == m_maxLhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_LOWER) 1132 { 1133 // this row leads to the tightest lower bound, slack should not be in the basis 1134 y[i] = y[m_i] * m_scale.value(k); 1135 y[m_i] = m_i_rowObj; 1136 1137 rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_LOWER : 1138 SPxSolverBase<R>::ON_UPPER; 1139 1140 if(i != m_i) 1141 rStatus[m_i] = SPxSolverBase<R>::BASIC; 1142 1143 haveSetBasis = true; 1144 } 1145 else if(i == m_minRhsIdx && rStatus[m_i] == SPxSolverBase<R>::ON_UPPER) 1146 { 1147 // this row leads to the tightest upper bound, slack should not be in the basis 1148 y[i] = y[m_i] * m_scale.value(k); 1149 y[m_i] = m_i_rowObj; 1150 1151 rStatus[i] = m_scale.value(k) * m_scale.value(0) > 0 ? SPxSolverBase<R>::ON_UPPER : 1152 SPxSolverBase<R>::ON_LOWER; 1153 1154 if(i != m_i) 1155 rStatus[m_i] = SPxSolverBase<R>::BASIC; 1156 1157 haveSetBasis = true; 1158 } 1159 else if(i != m_i) 1160 { 1161 // this row does not lead to the tightest lower or upper bound, slack should be in the basis 1162 y[i] = m_rowObj.value(k); 1163 rStatus[i] = SPxSolverBase<R>::BASIC; 1164 } 1165 } 1166 1167 #ifdef SOPLEX_CHECK_BASIS_DIM 1168 1169 if(m_isFirst && !this->checkBasisDim(rStatus, cStatus)) 1170 { 1171 throw SPxInternalCodeException("XMAISM24 Dimension doesn't match after this step."); 1172 } 1173 1174 #endif 1175 1176 // nothing to do for the reduced cost values 1177 } 1178 1179 template <class R> 1180 void SPxMainSM<R>::DuplicateColsPS::execute(VectorBase<R>& x, 1181 VectorBase<R>&, 1182 VectorBase<R>&, 1183 VectorBase<R>& r, 1184 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1185 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1186 { 1187 1188 if(m_isFirst) 1189 { 1190 #ifdef SOPLEX_CHECK_BASIS_DIM 1191 1192 if(!this->checkBasisDim(rStatus, cStatus)) 1193 { 1194 throw SPxInternalCodeException("XMAISM25 Dimension doesn't match after this step."); 1195 } 1196 1197 #endif 1198 return; 1199 } 1200 1201 1202 // correcting the change of idx by deletion of the columns: 1203 if(m_isLast) 1204 { 1205 for(int i = m_perm.size() - 1; i >= 0; --i) 1206 { 1207 if(m_perm[i] >= 0) 1208 { 1209 int cIdx_new = m_perm[i]; 1210 int cIdx = i; 1211 x[cIdx] = x[cIdx_new]; 1212 r[cIdx] = r[cIdx_new]; 1213 cStatus[cIdx] = cStatus[cIdx_new]; 1214 } 1215 } 1216 1217 return; 1218 } 1219 1220 // primal & basis: 1221 SOPLEX_ASSERT_WARN("WMAISM03", isNotZero(m_scale, this->epsilon())); 1222 1223 if(cStatus[m_k] == SPxSolverBase<R>::ON_LOWER) 1224 { 1225 x[m_k] = m_loK; 1226 1227 if(m_scale > 0) 1228 { 1229 x[m_j] = m_loJ; 1230 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1231 } 1232 else 1233 { 1234 x[m_j] = m_upJ; 1235 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1236 } 1237 } 1238 else if(cStatus[m_k] == SPxSolverBase<R>::ON_UPPER) 1239 { 1240 x[m_k] = m_upK; 1241 1242 if(m_scale > 0) 1243 { 1244 x[m_j] = m_upJ; 1245 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1246 } 1247 else 1248 { 1249 x[m_j] = m_loJ; 1250 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1251 } 1252 } 1253 else if(cStatus[m_k] == SPxSolverBase<R>::FIXED) 1254 { 1255 // => x[m_k] and x[m_j] are also fixed before the corresponding preprocessing step 1256 x[m_j] = m_loJ; 1257 cStatus[m_j] = SPxSolverBase<R>::FIXED; 1258 } 1259 else if(cStatus[m_k] == SPxSolverBase<R>::ZERO) 1260 { 1261 /* we only aggregate duplicate columns if 0 is contained in their bounds, so we can handle this case properly */ 1262 assert(isZero(x[m_k], this->epsilon())); 1263 assert(LErel(m_loJ, R(0.0), this->epsilon())); 1264 assert(GErel(m_upJ, R(0.0), this->epsilon())); 1265 assert(LErel(m_loK, R(0.0), this->epsilon())); 1266 assert(GErel(m_upK, R(0.0), this->epsilon())); 1267 1268 if(isZero(m_loK, this->epsilon()) && isZero(m_upK, this->epsilon()) && m_loK == m_upK) 1269 cStatus[m_k] = SPxSolverBase<R>::FIXED; 1270 else if(isZero(m_loK, this->epsilon())) 1271 cStatus[m_k] = SPxSolverBase<R>::ON_LOWER; 1272 else if(isZero(m_upK, this->epsilon())) 1273 cStatus[m_k] = SPxSolverBase<R>::ON_UPPER; 1274 else if(LErel(m_loK, R(0.0), this->epsilon()) && GErel(m_upK, R(0.0), this->epsilon())) 1275 cStatus[m_k] = SPxSolverBase<R>::ZERO; 1276 else 1277 throw SPxInternalCodeException("XMAISM05 This should never happen."); 1278 1279 x[m_j] = 0.0; 1280 1281 if(isZero(m_loJ, this->epsilon()) && isZero(m_upJ, this->epsilon()) && m_loJ == m_upJ) 1282 cStatus[m_j] = SPxSolverBase<R>::FIXED; 1283 else if(isZero(m_loJ, this->epsilon())) 1284 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 1285 else if(isZero(m_upJ, this->epsilon())) 1286 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 1287 else if(LErel(m_loJ, R(0.0), this->epsilon()) && GErel(m_upJ, R(0.0), this->epsilon())) 1288 cStatus[m_j] = SPxSolverBase<R>::ZERO; 1289 else 1290 throw SPxInternalCodeException("XMAISM06 This should never happen."); 1291 } 1292 else if(cStatus[m_k] == SPxSolverBase<R>::BASIC) 1293 { 1294 R scale1 = maxAbs(x[m_k], m_loK); 1295 R scale2 = maxAbs(x[m_k], m_upK); 1296 1297 if(scale1 < 1.0) 1298 scale1 = 1.0; 1299 1300 if(scale2 < 1.0) 1301 scale2 = 1.0; 1302 1303 R z1 = (x[m_k] / scale1) - (m_loK / scale1); 1304 R z2 = (x[m_k] / scale2) - (m_upK / scale2); 1305 1306 if(isZero(z1, this->epsilon())) 1307 z1 = 0.0; 1308 1309 if(isZero(z2, this->epsilon())) 1310 z2 = 0.0; 1311 1312 if(m_loJ <= R(-infinity) && m_upJ >= R(infinity) && m_loK <= R(-infinity) && m_upK >= R(infinity)) 1313 { 1314 cStatus[m_j] = SPxSolverBase<R>::ZERO; 1315 x[m_j] = 0.0; 1316 } 1317 else if(m_scale > 0.0) 1318 { 1319 if(GErel(x[m_k], m_upK + m_scale * m_upJ, this->epsilon())) 1320 { 1321 assert(m_upJ < R(infinity)); 1322 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1323 x[m_j] = m_upJ; 1324 x[m_k] -= m_scale * x[m_j]; 1325 } 1326 else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_upJ < R(infinity)) 1327 { 1328 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1329 x[m_j] = m_upJ; 1330 x[m_k] -= m_scale * x[m_j]; 1331 } 1332 else if(GErel(x[m_k], m_upK + m_scale * m_loJ, this->epsilon()) && m_upK < R(infinity)) 1333 { 1334 cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1335 x[m_k] = m_upK; 1336 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1337 x[m_j] = z2 * scale2 / m_scale; 1338 } 1339 else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loJ > R(-infinity)) 1340 { 1341 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1342 x[m_j] = m_loJ; 1343 x[m_k] -= m_scale * x[m_j]; 1344 } 1345 else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loK > R(-infinity)) 1346 { 1347 cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1348 x[m_k] = m_loK; 1349 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1350 x[m_j] = z1 * scale1 / m_scale; 1351 } 1352 else if(LTrel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon())) 1353 { 1354 assert(m_loJ > R(-infinity)); 1355 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1356 x[m_j] = m_loJ; 1357 x[m_k] -= m_scale * x[m_j]; 1358 } 1359 else 1360 { 1361 throw SPxInternalCodeException("XMAISM08 This should never happen."); 1362 } 1363 } 1364 else 1365 { 1366 assert(m_scale < 0.0); 1367 1368 if(GErel(x[m_k], m_upK + m_scale * m_loJ, this->epsilon())) 1369 { 1370 assert(m_loJ > R(-infinity)); 1371 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1372 x[m_j] = m_loJ; 1373 x[m_k] -= m_scale * x[m_j]; 1374 } 1375 else if(GErel(x[m_k], m_loK + m_scale * m_loJ, this->epsilon()) && m_loJ > R(-infinity)) 1376 { 1377 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1378 x[m_j] = m_loJ; 1379 x[m_k] -= m_scale * x[m_j]; 1380 } 1381 else if(GErel(x[m_k], m_upK + m_scale * m_upJ, this->epsilon()) && m_upK < R(infinity)) 1382 { 1383 cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1384 x[m_k] = m_upK; 1385 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1386 x[m_j] = z2 * scale2 / m_scale; 1387 } 1388 else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_upJ < R(infinity)) 1389 { 1390 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1391 x[m_j] = m_upJ; 1392 x[m_k] -= m_scale * x[m_j]; 1393 } 1394 else if(GErel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon()) && m_loK > R(-infinity)) 1395 { 1396 cStatus[m_k] = (m_loK == m_upK) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_LOWER; 1397 x[m_k] = m_loK; 1398 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1399 x[m_j] = z1 * scale1 / m_scale; 1400 } 1401 else if(LTrel(x[m_k], m_loK + m_scale * m_upJ, this->epsilon())) 1402 { 1403 assert(m_upJ < R(infinity)); 1404 cStatus[m_j] = (m_loJ == m_upJ) ? SPxSolverBase<R>::FIXED : SPxSolverBase<R>::ON_UPPER; 1405 x[m_j] = m_upJ; 1406 x[m_k] -= m_scale * x[m_j]; 1407 } 1408 else 1409 { 1410 throw SPxInternalCodeException("XMAISM09 This should never happen."); 1411 } 1412 } 1413 } 1414 1415 // dual: 1416 r[m_j] = m_scale * r[m_k]; 1417 } 1418 1419 template <class R> 1420 void SPxMainSM<R>::AggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 1421 VectorBase<R>& r, 1422 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1423 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1424 { 1425 // correcting the change of idx by deletion of the row: 1426 if(m_i != m_old_i) 1427 { 1428 s[m_old_i] = s[m_i]; 1429 y[m_old_i] = y[m_i]; 1430 rStatus[m_old_i] = rStatus[m_i]; 1431 } 1432 1433 // correcting the change of idx by deletion of the column: 1434 if(m_j != m_old_j) 1435 { 1436 x[m_old_j] = x[m_j]; 1437 r[m_old_j] = r[m_j]; 1438 cStatus[m_old_j] = cStatus[m_j]; 1439 } 1440 1441 // primal: 1442 R val = 0.0; 1443 R aij = m_row[m_j]; 1444 int active_idx = -1; 1445 1446 assert(m_row.size() == 2); 1447 1448 for(int k = 0; k < 2; ++k) 1449 { 1450 if(m_row.index(k) != m_j) 1451 { 1452 active_idx = m_row.index(k); 1453 val = m_row.value(k) * x[active_idx]; 1454 } 1455 } 1456 1457 assert(active_idx >= 0); 1458 1459 R scale = maxAbs(m_rhs, val); 1460 1461 if(scale < 1.0) 1462 scale = 1.0; 1463 1464 R z = (m_rhs / scale) - (val / scale); 1465 1466 if(isZero(z, this->epsilon())) 1467 z = 0.0; 1468 1469 x[m_j] = z * scale / aij; 1470 s[m_i] = m_rhs; 1471 1472 if(isOptimal && (LT(x[m_j], m_lower, this->feastol()) || GT(x[m_j], m_upper, this->feastol()))) 1473 { 1474 SPX_MSG_ERROR(std::cerr << "EMAISM: numerical violation after disaggregating variable" << std::endl; 1475 ) 1476 } 1477 1478 // dual: 1479 R dualVal = 0.0; 1480 1481 for(int k = 0; k < m_col.size(); ++k) 1482 { 1483 if(m_col.index(k) != m_i) 1484 dualVal += m_col.value(k) * y[m_col.index(k)]; 1485 } 1486 1487 z = m_obj - dualVal; 1488 1489 y[m_i] = z / aij; 1490 r[m_j] = 0.0; 1491 1492 // basis: 1493 if(((cStatus[active_idx] == SPxSolverBase<R>::ON_UPPER 1494 || cStatus[active_idx] == SPxSolverBase<R>::FIXED) 1495 && NE(x[active_idx], m_oldupper, this->feastol())) || 1496 ((cStatus[active_idx] == SPxSolverBase<R>::ON_LOWER 1497 || cStatus[active_idx] == SPxSolverBase<R>::FIXED) 1498 && NE(x[active_idx], m_oldlower, this->feastol()))) 1499 { 1500 cStatus[active_idx] = SPxSolverBase<R>::BASIC; 1501 r[active_idx] = 0.0; 1502 assert(NE(m_upper, m_lower, this->epsilon())); 1503 1504 if(EQ(x[m_j], m_upper, this->feastol())) 1505 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 1506 else if(EQ(x[m_j], m_lower, this->feastol())) 1507 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 1508 else if(m_upper >= R(infinity) && m_lower <= R(-infinity)) 1509 cStatus[m_j] = SPxSolverBase<R>::ZERO; 1510 else 1511 throw SPxInternalCodeException("XMAISM unexpected basis status in aggregation unsimplifier."); 1512 } 1513 else 1514 { 1515 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1516 } 1517 1518 // sides may not be equal and we always only consider the rhs during aggregation, so set ON_UPPER 1519 // (in theory and with exact arithmetic setting it to FIXED would be correct) 1520 rStatus[m_i] = SPxSolverBase<R>::ON_UPPER; 1521 1522 #ifdef SOPLEX_CHECK_BASIS_DIM 1523 1524 if(!this->checkBasisDim(rStatus, cStatus)) 1525 { 1526 throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step."); 1527 } 1528 1529 #endif 1530 } 1531 1532 template <class R> 1533 void SPxMainSM<R>::MultiAggregationPS::execute(VectorBase<R>& x, VectorBase<R>& y, VectorBase<R>& s, 1534 VectorBase<R>& r, 1535 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1536 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1537 { 1538 1539 // correcting the change of idx by deletion of the row: 1540 if(m_i != m_old_i) 1541 { 1542 s[m_old_i] = s[m_i]; 1543 y[m_old_i] = y[m_i]; 1544 rStatus[m_old_i] = rStatus[m_i]; 1545 } 1546 1547 // correcting the change of idx by deletion of the column: 1548 if(m_j != m_old_j) 1549 { 1550 x[m_old_j] = x[m_j]; 1551 r[m_old_j] = r[m_j]; 1552 cStatus[m_old_j] = cStatus[m_j]; 1553 } 1554 1555 // primal: 1556 R val = 0.0; 1557 R aij = m_row[m_j]; 1558 1559 for(int k = 0; k < m_row.size(); ++k) 1560 { 1561 if(m_row.index(k) != m_j) 1562 val += m_row.value(k) * x[m_row.index(k)]; 1563 } 1564 1565 R scale = maxAbs(m_const, val); 1566 1567 if(scale < 1.0) 1568 scale = 1.0; 1569 1570 R z = (m_const / scale) - (val / scale); 1571 1572 if(isZero(z, this->epsilon())) 1573 z = 0.0; 1574 1575 x[m_j] = z * scale / aij; 1576 s[m_i] = 0.0; 1577 1578 #ifndef NDEBUG 1579 1580 if(isOptimal && (LT(x[m_j], m_lower, this->feastol()) || GT(x[m_j], m_upper, this->feastol()))) 1581 SPX_MSG_ERROR(std::cerr << "numerical violation in original space due to MultiAggregation\n";) 1582 #endif 1583 1584 // dual: 1585 R dualVal = 0.0; 1586 1587 for(int k = 0; k < m_col.size(); ++k) 1588 { 1589 if(m_col.index(k) != m_i) 1590 dualVal += m_col.value(k) * y[m_col.index(k)]; 1591 } 1592 1593 z = m_obj - dualVal; 1594 1595 y[m_i] = z / aij; 1596 r[m_j] = 0.0; 1597 1598 // basis: 1599 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1600 1601 if(m_eqCons) 1602 rStatus[m_i] = SPxSolverBase<R>::FIXED; 1603 else if(m_onLhs) 1604 rStatus[m_i] = SPxSolverBase<R>::ON_LOWER; 1605 else 1606 rStatus[m_i] = SPxSolverBase<R>::ON_UPPER; 1607 1608 #ifdef SOPLEX_CHECK_BASIS_DIM 1609 1610 if(!this->checkBasisDim(rStatus, cStatus)) 1611 { 1612 throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step."); 1613 } 1614 1615 #endif 1616 } 1617 1618 template <class R> 1619 void SPxMainSM<R>::TightenBoundsPS::execute(VectorBase<R>& x, VectorBase<R>&, VectorBase<R>&, 1620 VectorBase<R>&, 1621 DataArray<typename SPxSolverBase<R>::VarStatus>& cStatus, 1622 DataArray<typename SPxSolverBase<R>::VarStatus>& rStatus, bool isOptimal) const 1623 { 1624 // basis: 1625 switch(cStatus[m_j]) 1626 { 1627 case SPxSolverBase<R>::FIXED: 1628 if(LT(x[m_j], m_origupper, this->feastol()) && GT(x[m_j], m_origlower, this->feastol())) 1629 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1630 else if(LT(x[m_j], m_origupper, this->feastol())) 1631 cStatus[m_j] = SPxSolverBase<R>::ON_LOWER; 1632 else if(GT(x[m_j], m_origlower, this->feastol())) 1633 cStatus[m_j] = SPxSolverBase<R>::ON_UPPER; 1634 1635 break; 1636 1637 case SPxSolverBase<R>::ON_LOWER: 1638 if(GT(x[m_j], m_origlower, this->feastol())) 1639 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1640 1641 break; 1642 1643 case SPxSolverBase<R>::ON_UPPER: 1644 if(LT(x[m_j], m_origupper, this->feastol())) 1645 cStatus[m_j] = SPxSolverBase<R>::BASIC; 1646 1647 break; 1648 1649 default: 1650 break; 1651 } 1652 1653 #ifdef SOPLEX_CHECK_BASIS_DIM 1654 1655 if(!this->checkBasisDim(rStatus, cStatus)) 1656 { 1657 throw SPxInternalCodeException("XMAISM22 Dimension doesn't match after this step."); 1658 } 1659 1660 #endif 1661 } 1662 1663 template <class R> 1664 void SPxMainSM<R>::handleRowObjectives(SPxLPBase<R>& lp) 1665 { 1666 for(int i = lp.nRows() - 1; i >= 0; --i) 1667 { 1668 if(lp.maxRowObj(i) != 0.0) 1669 { 1670 std::shared_ptr<PostStep> ptr(new RowObjPS(lp, i, lp.nCols(), this->_tolerances)); 1671 m_hist.append(ptr); 1672 lp.addCol(lp.rowObj(i), -lp.rhs(i), UnitVectorBase<R>(i), -lp.lhs(i)); 1673 lp.changeRange(i, R(0.0), R(0.0)); 1674 lp.changeRowObj(i, R(0.0)); 1675 m_addedcols++; 1676 } 1677 } 1678 } 1679 1680 template <class R> 1681 void SPxMainSM<R>::handleExtremes(SPxLPBase<R>& lp) 1682 { 1683 1684 // This method handles extreme value of the given LP by 1685 // 1686 // 1. setting numbers of very small absolute values to zero and 1687 // 2. setting numbers of very large absolute values to R(-infinity) or +R(infinity), respectively. 1688 1689 R maxVal = R(infinity) / 5.0; 1690 R tol = feastol() * 1e-2; 1691 tol = (tol < this->epsZero()) ? this->epsZero() : tol; 1692 int remRows = 0; 1693 int remNzos = 0; 1694 int chgBnds = 0; 1695 int chgLRhs = 0; 1696 int objCnt = 0; 1697 1698 for(int i = lp.nRows() - 1; i >= 0; --i) 1699 { 1700 // lhs 1701 R lhs = lp.lhs(i); 1702 1703 if(lhs != 0.0 && isZero(lhs, this->epsZero())) 1704 { 1705 lp.changeLhs(i, R(0.0)); 1706 ++chgLRhs; 1707 } 1708 else if(lhs > R(-infinity) && lhs < -maxVal) 1709 { 1710 lp.changeLhs(i, R(-infinity)); 1711 ++chgLRhs; 1712 } 1713 else if(lhs < R(infinity) && lhs > maxVal) 1714 { 1715 lp.changeLhs(i, R(infinity)); 1716 ++chgLRhs; 1717 } 1718 1719 // rhs 1720 R rhs = lp.rhs(i); 1721 1722 if(rhs != 0.0 && isZero(rhs, this->epsZero())) 1723 { 1724 lp.changeRhs(i, R(0.0)); 1725 ++chgLRhs; 1726 } 1727 else if(rhs > R(-infinity) && rhs < -maxVal) 1728 { 1729 lp.changeRhs(i, R(-infinity)); 1730 ++chgLRhs; 1731 } 1732 else if(rhs < R(infinity) && rhs > maxVal) 1733 { 1734 lp.changeRhs(i, R(infinity)); 1735 ++chgLRhs; 1736 } 1737 1738 if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity)) 1739 { 1740 std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances)); 1741 m_hist.append(ptr); 1742 1743 removeRow(lp, i); 1744 ++remRows; 1745 1746 ++m_stat[FREE_ROW]; 1747 } 1748 } 1749 1750 for(int j = 0; j < lp.nCols(); ++j) 1751 { 1752 // lower bound 1753 R lo = lp.lower(j); 1754 1755 if(lo != 0.0 && isZero(lo, this->epsZero())) 1756 { 1757 lp.changeLower(j, R(0.0)); 1758 ++chgBnds; 1759 } 1760 else if(lo > R(-infinity) && lo < -maxVal) 1761 { 1762 lp.changeLower(j, R(-infinity)); 1763 ++chgBnds; 1764 } 1765 else if(lo < R(infinity) && lo > maxVal) 1766 { 1767 lp.changeLower(j, R(infinity)); 1768 ++chgBnds; 1769 } 1770 1771 // upper bound 1772 R up = lp.upper(j); 1773 1774 if(up != 0.0 && isZero(up, this->epsZero())) 1775 { 1776 lp.changeUpper(j, R(0.0)); 1777 ++chgBnds; 1778 } 1779 else if(up > R(-infinity) && up < -maxVal) 1780 { 1781 lp.changeUpper(j, R(-infinity)); 1782 ++chgBnds; 1783 } 1784 else if(up < R(infinity) && up > maxVal) 1785 { 1786 lp.changeUpper(j, R(infinity)); 1787 ++chgBnds; 1788 } 1789 1790 // fixed columns will be eliminated later 1791 if(NE(lo, up, this->tolerances()->epsilon())) 1792 { 1793 lo = spxAbs(lo); 1794 up = spxAbs(up); 1795 1796 R absBnd = (lo > up) ? lo : up; 1797 1798 if(absBnd < 1.0) 1799 absBnd = 1.0; 1800 1801 // non-zeros 1802 SVectorBase<R>& col = lp.colVector_w(j); 1803 int i = 0; 1804 1805 while(i < col.size()) 1806 { 1807 R aij = spxAbs(col.value(i)); 1808 1809 if(isZero(aij * absBnd, tol)) 1810 { 1811 SVectorBase<R>& row = lp.rowVector_w(col.index(i)); 1812 int row_j = row.pos(j); 1813 1814 // this changes col.size() 1815 if(row_j >= 0) 1816 row.remove(row_j); 1817 1818 col.remove(i); 1819 1820 SPxOut::debug(this, "IMAISM04 aij={} removed, absBnd={} \n", aij, absBnd); 1821 ++remNzos; 1822 } 1823 else 1824 { 1825 if(aij > maxVal) 1826 { 1827 SPX_MSG_WARNING((*this->spxout), 1828 (*this->spxout) << "WMAISM05 Warning! Big matrix coefficient " << aij 1829 << std::endl); 1830 } 1831 else if(isZero(aij, tol)) 1832 { 1833 SPX_MSG_WARNING((*this->spxout), 1834 (*this->spxout) << "WMAISM06 Warning! Tiny matrix coefficient " << aij 1835 << std::endl); 1836 } 1837 1838 ++i; 1839 } 1840 } 1841 } 1842 1843 // objective 1844 R obj = lp.obj(j); 1845 1846 if(obj != 0.0 && isZero(obj, this->epsZero())) 1847 { 1848 lp.changeObj(j, R(0.0)); 1849 ++objCnt; 1850 } 1851 else if(obj > R(-infinity) && obj < -maxVal) 1852 { 1853 lp.changeObj(j, R(-infinity)); 1854 ++objCnt; 1855 } 1856 else if(obj < R(infinity) && obj > maxVal) 1857 { 1858 lp.changeObj(j, R(infinity)); 1859 ++objCnt; 1860 } 1861 } 1862 1863 if(remRows + remNzos + chgLRhs + chgBnds + objCnt > 0) 1864 { 1865 this->m_remRows += remRows; 1866 this->m_remNzos += remNzos; 1867 this->m_chgLRhs += chgLRhs; 1868 this->m_chgBnds += chgBnds; 1869 1870 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (extremes) removed " 1871 << remRows << " rows, " 1872 << remNzos << " non-zeros, " 1873 << chgBnds << " col bounds, " 1874 << chgLRhs << " row bounds, " 1875 << objCnt << " objective coefficients" << std::endl;) 1876 } 1877 1878 assert(lp.isConsistent()); 1879 } 1880 1881 /// computes the minimum and maximum residual activity for a given variable 1882 template <class R> 1883 void SPxMainSM<R>::computeMinMaxResidualActivity(SPxLPBase<R>& lp, int rowNumber, int colNumber, 1884 R& minAct, R& maxAct) 1885 { 1886 const SVectorBase<R>& row = lp.rowVector(rowNumber); 1887 bool minNegInfinite = false; 1888 bool maxInfinite = false; 1889 1890 minAct = 0; // this is the minimum value that the aggregation can attain 1891 maxAct = 0; // this is the maximum value that the aggregation can attain 1892 1893 for(int l = 0; l < row.size(); ++l) 1894 { 1895 if(colNumber < 0 || row.index(l) != colNumber) 1896 { 1897 // computing the minimum activity of the aggregated variables 1898 if(GT(row.value(l), R(0.0), this->tolerances()->epsilon())) 1899 { 1900 if(lp.lower(row.index(l)) <= R(-infinity)) 1901 minNegInfinite = true; 1902 else 1903 minAct += row.value(l) * lp.lower(row.index(l)); 1904 } 1905 else if(LT(row.value(l), R(0.0), this->tolerances()->epsilon())) 1906 { 1907 if(lp.upper(row.index(l)) >= R(infinity)) 1908 minNegInfinite = true; 1909 else 1910 minAct += row.value(l) * lp.upper(row.index(l)); 1911 } 1912 1913 // computing the maximum activity of the aggregated variables 1914 if(GT(row.value(l), R(0.0), this->tolerances()->epsilon())) 1915 { 1916 if(lp.upper(row.index(l)) >= R(infinity)) 1917 maxInfinite = true; 1918 else 1919 maxAct += row.value(l) * lp.upper(row.index(l)); 1920 } 1921 else if(LT(row.value(l), R(0.0), this->tolerances()->epsilon())) 1922 { 1923 if(lp.lower(row.index(l)) <= R(-infinity)) 1924 maxInfinite = true; 1925 else 1926 maxAct += row.value(l) * lp.lower(row.index(l)); 1927 } 1928 } 1929 } 1930 1931 // if an infinite value exists for the minimum activity, then that it taken 1932 if(minNegInfinite) 1933 minAct = R(-infinity); 1934 1935 // if an -infinite value exists for the maximum activity, then that value is taken 1936 if(maxInfinite) 1937 maxAct = R(infinity); 1938 } 1939 1940 1941 /// calculate min/max value for the multi aggregated variables 1942 template <class R> 1943 void SPxMainSM<R>::computeMinMaxValues(SPxLPBase<R>& lp, R side, R val, R minRes, R maxRes, 1944 R& minVal, R& maxVal) 1945 { 1946 minVal = 0; 1947 maxVal = 0; 1948 1949 if(LT(val, R(0.0), this->tolerances()->epsilon())) 1950 { 1951 if(minRes <= R(-infinity)) 1952 minVal = R(-infinity); 1953 else 1954 minVal = (side - minRes) / val; 1955 1956 if(maxRes >= R(infinity)) 1957 maxVal = R(infinity); 1958 else 1959 maxVal = (side - maxRes) / val; 1960 } 1961 else if(GT(val, R(0.0), this->tolerances()->epsilon())) 1962 { 1963 if(maxRes >= R(infinity)) 1964 minVal = R(-infinity); 1965 else 1966 minVal = (side - maxRes) / val; 1967 1968 if(minRes <= R(-infinity)) 1969 maxVal = R(infinity); 1970 else 1971 maxVal = (side - minRes) / val; 1972 } 1973 } 1974 1975 1976 /// tries to find good lower bound solutions by applying some trivial heuristics 1977 template <class R> 1978 void SPxMainSM<R>::trivialHeuristic(SPxLPBase<R>& lp) 1979 { 1980 VectorBase<R> zerosol(lp.nCols()); // the zero solution VectorBase<R> 1981 VectorBase<R> lowersol(lp.nCols()); // the lower bound solution VectorBase<R> 1982 VectorBase<R> uppersol(lp.nCols()); // the upper bound solution VectorBase<R> 1983 VectorBase<R> locksol(lp.nCols()); // the locks solution VectorBase<R> 1984 1985 VectorBase<R> upLocks(lp.nCols()); 1986 VectorBase<R> downLocks(lp.nCols()); 1987 1988 R zeroObj = this->m_objoffset; 1989 R lowerObj = this->m_objoffset; 1990 R upperObj = this->m_objoffset; 1991 R lockObj = this->m_objoffset; 1992 1993 bool zerovalid = true; 1994 1995 R largeValue = R(infinity); 1996 1997 if(LT(R(1.0 / feastol()), R(infinity), this->tolerances()->epsilon())) 1998 largeValue = 1.0 / feastol(); 1999 2000 2001 2002 for(int j = lp.nCols() - 1; j >= 0; --j) 2003 { 2004 upLocks[j] = 0; 2005 downLocks[j] = 0; 2006 2007 // computing the locks on the variables 2008 const SVectorBase<R>& col = lp.colVector(j); 2009 2010 for(int k = 0; k < col.size(); ++k) 2011 { 2012 R aij = col.value(k); 2013 2014 SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity)))); 2015 2016 if(GT(lp.lhs(col.index(k)), R(-infinity), this->tolerances()->epsilon()) 2017 && LT(lp.rhs(col.index(k)), R(infinity), this->tolerances()->epsilon())) 2018 { 2019 upLocks[j]++; 2020 downLocks[j]++; 2021 } 2022 else if(GT(lp.lhs(col.index(k)), R(-infinity), this->tolerances()->epsilon())) 2023 { 2024 if(aij > 0) 2025 downLocks[j]++; 2026 else if(aij < 0) 2027 upLocks[j]++; 2028 } 2029 else if(LT(lp.rhs(col.index(k)), R(infinity), this->tolerances()->epsilon())) 2030 { 2031 if(aij > 0) 2032 upLocks[j]++; 2033 else if(aij < 0) 2034 downLocks[j]++; 2035 } 2036 } 2037 2038 R lower = lp.lower(j); 2039 R upper = lp.upper(j); 2040 2041 if(LE(lower, R(-infinity), this->tolerances()->epsilon())) 2042 lower = SOPLEX_MIN(-largeValue, upper); 2043 2044 if(GE(upper, R(infinity), this->tolerances()->epsilon())) 2045 upper = SOPLEX_MAX(lp.lower(j), largeValue); 2046 2047 if(zerovalid) 2048 { 2049 if(LE(lower, R(0.0), feastol()) && GE(upper, R(0.0), feastol())) 2050 zerosol[j] = 0.0; 2051 else 2052 zerovalid = false; 2053 } 2054 2055 lowersol[j] = lower; 2056 uppersol[j] = upper; 2057 2058 if(downLocks[j] > upLocks[j]) 2059 locksol[j] = upper; 2060 else if(downLocks[j] < upLocks[j]) 2061 locksol[j] = lower; 2062 else 2063 locksol[j] = (lower + upper) / 2.0; 2064 2065 lowerObj += lp.maxObj(j) * lowersol[j]; 2066 upperObj += lp.maxObj(j) * uppersol[j]; 2067 lockObj += lp.maxObj(j) * locksol[j]; 2068 } 2069 2070 // trying the lower bound solution 2071 if(checkSolution(lp, lowersol)) 2072 { 2073 if(lowerObj > m_cutoffbound) 2074 m_cutoffbound = lowerObj; 2075 } 2076 2077 // trying the upper bound solution 2078 if(checkSolution(lp, uppersol)) 2079 { 2080 if(upperObj > m_cutoffbound) 2081 m_cutoffbound = upperObj; 2082 } 2083 2084 // trying the zero solution 2085 if(zerovalid && checkSolution(lp, zerosol)) 2086 { 2087 if(zeroObj > m_cutoffbound) 2088 m_cutoffbound = zeroObj; 2089 } 2090 2091 // trying the lock solution 2092 if(checkSolution(lp, locksol)) 2093 { 2094 if(lockObj > m_cutoffbound) 2095 m_cutoffbound = lockObj; 2096 } 2097 } 2098 2099 2100 2101 /// checks a solution for feasibility 2102 template <class R> 2103 bool SPxMainSM<R>::checkSolution(SPxLPBase<R>& lp, VectorBase<R> sol) 2104 { 2105 for(int i = lp.nRows() - 1; i >= 0; --i) 2106 { 2107 const SVectorBase<R>& row = lp.rowVector(i); 2108 R activity = 0; 2109 2110 for(int k = 0; k < row.size(); k++) 2111 activity += row.value(k) * sol[row.index(k)]; 2112 2113 if(!GE(activity, lp.lhs(i), feastol()) || !LE(activity, lp.rhs(i), feastol())) 2114 return false; 2115 } 2116 2117 return true; 2118 } 2119 2120 2121 /// tightens variable bounds by propagating the pseudo objective function value. 2122 template <class R> 2123 void SPxMainSM<R>::propagatePseudoobj(SPxLPBase<R>& lp) 2124 { 2125 R pseudoObj = this->m_objoffset; 2126 2127 for(int j = lp.nCols() - 1; j >= 0; --j) 2128 { 2129 R val = lp.maxObj(j); 2130 2131 if(val < 0) 2132 { 2133 if(lp.lower(j) <= R(-infinity)) 2134 return; 2135 2136 pseudoObj += val * lp.lower(j); 2137 } 2138 else if(val > 0) 2139 { 2140 if(lp.upper(j) >= R(-infinity)) 2141 return; 2142 2143 pseudoObj += val * lp.upper(j); 2144 } 2145 } 2146 2147 if(GT(m_cutoffbound, R(-infinity), this->tolerances()->epsilon()) 2148 && LT(m_cutoffbound, R(infinity), this->tolerances()->epsilon())) 2149 { 2150 if(pseudoObj > m_pseudoobj) 2151 m_pseudoobj = pseudoObj; 2152 2153 for(int j = lp.nCols() - 1; j >= 0; --j) 2154 { 2155 R objval = lp.maxObj(j); 2156 2157 if(EQ(objval, R(0.0), this->tolerances()->epsilon())) 2158 continue; 2159 2160 if(objval < 0.0) 2161 { 2162 R newbound = lp.lower(j) + (m_cutoffbound - m_pseudoobj) / objval; 2163 2164 if(LT(newbound, lp.upper(j), this->tolerances()->epsilon())) 2165 { 2166 std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j), 2167 this->_tolerances)); 2168 m_hist.append(ptr); 2169 lp.changeUpper(j, newbound); 2170 } 2171 } 2172 else if(objval > 0.0) 2173 { 2174 R newbound = lp.upper(j) + (m_cutoffbound - m_pseudoobj) / objval; 2175 2176 if(GT(newbound, lp.lower(j), this->tolerances()->epsilon())) 2177 { 2178 std::shared_ptr<PostStep> ptr(new TightenBoundsPS(lp, j, lp.upper(j), lp.lower(j), 2179 this->_tolerances)); 2180 m_hist.append(ptr); 2181 lp.changeLower(j, newbound); 2182 } 2183 } 2184 } 2185 } 2186 } 2187 2188 2189 2190 template <class R> 2191 typename SPxSimplifier<R>::Result SPxMainSM<R>::removeEmpty(SPxLPBase<R>& lp) 2192 { 2193 2194 // This method removes empty rows and columns from the LP. 2195 2196 int remRows = 0; 2197 int remCols = 0; 2198 2199 for(int i = lp.nRows() - 1; i >= 0; --i) 2200 { 2201 const SVectorBase<R>& row = lp.rowVector(i); 2202 2203 if(row.size() == 0) 2204 { 2205 SPxOut::debug(this, "IMAISM07 row {}: empty ->", i); 2206 2207 if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol())) 2208 { 2209 SPxOut::debug(this, " infeasible lhs={} rhs={}\n", lp.lhs(i), lp.rhs(i)); 2210 return this->INFEASIBLE; 2211 } 2212 2213 SPxOut::debug(this, " removed\n"); 2214 2215 std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i, this->_tolerances)); 2216 m_hist.append(ptr); 2217 2218 ++remRows; 2219 removeRow(lp, i); 2220 2221 ++m_stat[EMPTY_ROW]; 2222 } 2223 } 2224 2225 for(int j = lp.nCols() - 1; j >= 0; --j) 2226 { 2227 const SVectorBase<R>& col = lp.colVector(j); 2228 2229 if(col.size() == 0) 2230 { 2231 SPxOut::debug(this, "IMAISM08 col {}: empty -> maxObj={} lower={} upper={}", j, lp.maxObj(j), 2232 lp.lower(j), lp.upper(j)); 2233 2234 R val; 2235 2236 if(GT(lp.maxObj(j), R(0.0), this->epsZero())) 2237 { 2238 if(lp.upper(j) >= R(infinity)) 2239 { 2240 SPxOut::debug(this, " unbounded\n"); 2241 return this->UNBOUNDED; 2242 } 2243 2244 val = lp.upper(j); 2245 } 2246 else if(LT(lp.maxObj(j), R(0.0), this->epsZero())) 2247 { 2248 if(lp.lower(j) <= R(-infinity)) 2249 { 2250 SPxOut::debug(this, " unbounded\n"); 2251 return this->UNBOUNDED; 2252 } 2253 2254 val = lp.lower(j); 2255 } 2256 else 2257 { 2258 SOPLEX_ASSERT_WARN("WMAISM09", isZero(lp.maxObj(j), this->epsZero())); 2259 2260 // any value within the bounds is ok 2261 if(lp.lower(j) > R(-infinity)) 2262 val = lp.lower(j); 2263 else if(lp.upper(j) < R(infinity)) 2264 val = lp.upper(j); 2265 else 2266 val = 0.0; 2267 } 2268 2269 SPxOut::debug(this, " removed\n"); 2270 2271 std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val, this->_tolerances)); 2272 std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val, this->_tolerances)); 2273 m_hist.append(ptr1); 2274 m_hist.append(ptr2); 2275 2276 ++remCols; 2277 removeCol(lp, j); 2278 2279 ++m_stat[EMPTY_COL]; 2280 } 2281 } 2282 2283 if(remRows + remCols > 0) 2284 { 2285 this->m_remRows += remRows; 2286 this->m_remCols += remCols; 2287 2288 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (empty rows/colums) removed " 2289 << remRows << " rows, " 2290 << remCols << " cols" 2291 << std::endl;) 2292 2293 } 2294 2295 return this->OKAY; 2296 } 2297 2298 template <class R> 2299 typename SPxSimplifier<R>::Result SPxMainSM<R>::removeRowSingleton(SPxLPBase<R>& lp, 2300 const SVectorBase<R>& row, int& i) 2301 { 2302 assert(row.size() == 1); 2303 2304 R aij = row.value(0); 2305 int j = row.index(0); 2306 R lo = R(-infinity); 2307 R up = R(infinity); 2308 2309 SPxOut::debug(this, "IMAISM22 row {}: singleton -> val={} lhs={} rhs={}", i, aij, lp.lhs(i), 2310 lp.rhs(i)); 2311 2312 if(GT(aij, R(0.0), this->epsZero())) // aij > 0 2313 { 2314 lo = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) / aij; 2315 up = (lp.rhs(i) >= R(infinity)) ? R(infinity) : lp.rhs(i) / aij; 2316 } 2317 else if(LT(aij, R(0.0), this->epsZero())) // aij < 0 2318 { 2319 lo = (lp.rhs(i) >= R(infinity)) ? R(-infinity) : lp.rhs(i) / aij; 2320 up = (lp.lhs(i) <= R(-infinity)) ? R(infinity) : lp.lhs(i) / aij; 2321 } 2322 else if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol())) 2323 { 2324 // aij == 0, rhs < 0 or lhs > 0 2325 SPxOut::debug(this, " infeasible\n"); 2326 return this->INFEASIBLE; 2327 } 2328 2329 if(isZero(lo, this->epsZero())) 2330 lo = 0.0; 2331 2332 if(isZero(up, this->epsZero())) 2333 up = 0.0; 2334 2335 SPxOut::debug(this, " removed, lower={} ({}) upper={} ({})\n", lo, lp.lower(j), up, lp.upper(j)); 2336 2337 bool stricterUp = false; 2338 bool stricterLo = false; 2339 2340 R oldLo = lp.lower(j); 2341 R oldUp = lp.upper(j); 2342 2343 if(LTrel(up, lp.upper(j), feastol())) 2344 { 2345 lp.changeUpper(j, up); 2346 stricterUp = true; 2347 } 2348 2349 if(GTrel(lo, lp.lower(j), feastol())) 2350 { 2351 lp.changeLower(j, lo); 2352 stricterLo = true; 2353 } 2354 2355 std::shared_ptr<PostStep> ptr(new RowSingletonPS(lp, i, j, stricterLo, stricterUp, 2356 lp.lower(j), 2357 lp.upper(j), oldLo, oldUp, this->_tolerances)); 2358 m_hist.append(ptr); 2359 2360 removeRow(lp, i); 2361 2362 this->m_remRows++; 2363 this->m_remNzos++; 2364 ++m_stat[SINGLETON_ROW]; 2365 2366 return this->OKAY; 2367 } 2368 2369 /// aggregate variable x_j to x_j = (rhs - aik * x_k) / aij from row i: aij * x_j + aik * x_k = rhs 2370 template <class R> 2371 typename SPxSimplifier<R>::Result SPxMainSM<R>::aggregateVars(SPxLPBase<R>& lp, 2372 const SVectorBase<R>& row, int& i) 2373 { 2374 assert(row.size() == 2); 2375 assert(EQrel(lp.lhs(i), lp.rhs(i), feastol())); 2376 2377 R rhs = lp.rhs(i); 2378 assert(rhs < R(infinity) && rhs > R(-infinity)); 2379 2380 int j = row.index(0); 2381 int k = row.index(1); 2382 R aij = row.value(0); 2383 R aik = row.value(1); 2384 R lower_j = lp.lower(j); 2385 R upper_j = lp.upper(j); 2386 R lower_k = lp.lower(k); 2387 R upper_k = lp.upper(k); 2388 2389 // fixed variables should be removed by simplifyCols() 2390 if(EQrel(lower_j, upper_j, feastol()) || EQrel(lower_k, upper_k, feastol())) 2391 return this->OKAY; 2392 2393 assert(isNotZero(aij, this->epsZero()) && isNotZero(aik, this->epsZero())); 2394 2395 SPxOut::debug(this, "IMAISM23 row {}: doubleton equation -> {} x_{} + {} x_{} = {}", i, aij, j, aik, 2396 k, rhs); 2397 2398 // determine which variable can be aggregated without requiring bound tightening of the other variable 2399 R new_lo_j; 2400 R new_up_j; 2401 R new_lo_k; 2402 R new_up_k; 2403 2404 if(aij * aik < 0.0) 2405 { 2406 // orientation persists 2407 new_lo_j = (upper_k >= R(infinity)) ? R(-infinity) : (rhs - aik * upper_k) / aij; 2408 new_up_j = (lower_k <= R(-infinity)) ? R(infinity) : (rhs - aik * lower_k) / aij; 2409 new_lo_k = (upper_j >= R(infinity)) ? R(-infinity) : (rhs - aij * upper_j) / aik; 2410 new_up_k = (lower_j <= R(-infinity)) ? R(infinity) : (rhs - aij * lower_j) / aik; 2411 } 2412 else if(aij * aik > 0.0) 2413 { 2414 // orientation is reversed 2415 new_lo_j = (lower_k <= R(-infinity)) ? R(-infinity) : (rhs - aik * lower_k) / aij; 2416 new_up_j = (upper_k >= R(infinity)) ? R(infinity) : (rhs - aik * upper_k) / aij; 2417 new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : (rhs - aij * lower_j) / aik; 2418 new_up_k = (upper_j >= R(infinity)) ? R(infinity) : (rhs - aij * upper_j) / aik; 2419 } 2420 else 2421 throw SPxInternalCodeException("XMAISM12 This should never happen."); 2422 2423 bool flip_jk = false; 2424 2425 if(new_lo_j <= R(-infinity) && new_up_j >= R(infinity)) 2426 { 2427 // no bound tightening on x_j when x_k is aggregated 2428 flip_jk = true; 2429 } 2430 else if(new_lo_k <= R(-infinity) && new_up_k >= R(infinity)) 2431 { 2432 // no bound tightening on x_k when x_j is aggregated 2433 flip_jk = false; 2434 } 2435 else if(LE(new_lo_j, lower_j, this->tolerances()->epsilon()) 2436 && GE(new_up_j, upper_j, this->tolerances()->epsilon())) 2437 { 2438 if(LE(new_lo_k, lower_k, this->tolerances()->epsilon()) 2439 && GE(new_up_k, upper_k, this->tolerances()->epsilon())) 2440 { 2441 // both variables' bounds are not affected by aggregation; choose the better aggregation coeff (aik/aij) 2442 if(spxAbs(aij) > spxAbs(aik)) 2443 flip_jk = false; 2444 else 2445 flip_jk = true; 2446 } 2447 else 2448 flip_jk = false; 2449 } 2450 else if(LE(new_lo_k, lower_k, this->tolerances()->epsilon()) 2451 && GE(new_up_k, upper_k, this->tolerances()->epsilon())) 2452 { 2453 flip_jk = true; 2454 } 2455 else 2456 { 2457 if(spxAbs(aij) > spxAbs(aik)) 2458 flip_jk = false; 2459 else 2460 flip_jk = true; 2461 } 2462 2463 if(flip_jk) 2464 { 2465 int _j = j; 2466 R _aij = aij; 2467 R _lower_j = lower_j; 2468 R _upper_j = upper_j; 2469 j = k; 2470 k = _j; 2471 aij = aik; 2472 aik = _aij; 2473 lower_j = lower_k; 2474 lower_k = _lower_j; 2475 upper_j = upper_k; 2476 upper_k = _upper_j; 2477 } 2478 2479 const SVectorBase<R>& col_j = lp.colVector(j); 2480 const SVectorBase<R>& col_k = lp.colVector(k); 2481 2482 // aggregation coefficients (x_j = aggr_coef * x_k + aggr_const) 2483 R aggr_coef = - (aik / aij); 2484 R aggr_const = rhs / aij; 2485 2486 SPxOut::debug(this, " removed, replacing x_{} with {} + {} * x_{}\n", j, aggr_const, aggr_coef, k); 2487 2488 // replace all occurrences of x_j 2489 for(int r = 0; r < col_j.size(); ++r) 2490 { 2491 int row_r = col_j.index(r); 2492 R arj = col_j.value(r); 2493 2494 // skip row i 2495 if(row_r == i) 2496 continue; 2497 2498 // adapt sides of row r 2499 R lhs_r = lp.lhs(row_r); 2500 R rhs_r = lp.rhs(row_r); 2501 2502 if(lhs_r > R(-infinity)) 2503 { 2504 lp.changeLhs(row_r, lhs_r - aggr_const * arj); 2505 this->m_chgLRhs++; 2506 } 2507 2508 if(rhs_r < R(infinity)) 2509 { 2510 lp.changeRhs(row_r, rhs_r - aggr_const * arj); 2511 this->m_chgLRhs++; 2512 } 2513 2514 R newcoef = aggr_coef * arj; 2515 int pos_rk = col_k.pos(row_r); 2516 2517 // check whether x_k is also present in row r and get its coefficient 2518 if(pos_rk >= 0) 2519 { 2520 R ark = col_k.value(pos_rk); 2521 newcoef += ark; 2522 this->m_remNzos++; 2523 } 2524 2525 // add new column k to row r or adapt the coefficient a_rk 2526 lp.changeElement(row_r, k, newcoef); 2527 } 2528 2529 // adapt objective function 2530 R obj_j = lp.obj(j); 2531 2532 if(isNotZero(obj_j, this->epsZero())) 2533 { 2534 this->addObjoffset(aggr_const * obj_j); 2535 R obj_k = lp.obj(k); 2536 lp.changeObj(k, obj_k + aggr_coef * obj_j); 2537 } 2538 2539 // adapt bounds of x_k 2540 R scale1 = maxAbs(rhs, aij * upper_j); 2541 R scale2 = maxAbs(rhs, aij * lower_j); 2542 2543 if(scale1 < 1.0) 2544 scale1 = 1.0; 2545 2546 if(scale2 < 1.0) 2547 scale2 = 1.0; 2548 2549 R z1 = (rhs / scale1) - (aij * upper_j / scale1); 2550 R z2 = (rhs / scale2) - (aij * lower_j / scale2); 2551 2552 // just some rounding 2553 if(isZero(z1, this->epsZero())) 2554 z1 = 0.0; 2555 2556 if(isZero(z2, this->epsZero())) 2557 z2 = 0.0; 2558 2559 // determine which side has to be used for the bounds comparison below 2560 if(aik * aij > 0.0) 2561 { 2562 new_lo_k = (upper_j >= R(infinity)) ? R(-infinity) : z1 * scale1 / aik; 2563 new_up_k = (lower_j <= R(-infinity)) ? R(infinity) : z2 * scale2 / aik; 2564 } 2565 else if(aik * aij < 0.0) 2566 { 2567 new_lo_k = (lower_j <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik; 2568 new_up_k = (upper_j >= R(infinity)) ? R(infinity) : z1 * scale1 / aik; 2569 } 2570 else 2571 throw SPxInternalCodeException("XMAISM12 This should never happen."); 2572 2573 // change bounds of x_k if the new ones are tighter 2574 R oldlower_k = lower_k; 2575 R oldupper_k = upper_k; 2576 2577 if(GT(new_lo_k, lower_k, this->epsZero())) 2578 { 2579 lp.changeLower(k, new_lo_k); 2580 this->m_chgBnds++; 2581 } 2582 2583 if(LT(new_up_k, upper_k, this->epsZero())) 2584 { 2585 lp.changeUpper(k, new_up_k); 2586 this->m_chgBnds++; 2587 } 2588 2589 std::shared_ptr<PostStep> ptr(new AggregationPS(lp, i, j, rhs, oldupper_k, oldlower_k, 2590 this->_tolerances)); 2591 m_hist.append(ptr); 2592 2593 removeRow(lp, i); 2594 removeCol(lp, j); 2595 2596 this->m_remRows++; 2597 this->m_remCols++; 2598 this->m_remNzos += 2; 2599 2600 ++m_stat[AGGREGATION]; 2601 2602 return this->OKAY; 2603 } 2604 2605 template <class R> 2606 typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyRows(SPxLPBase<R>& lp, bool& again) 2607 { 2608 2609 // This method simplifies the rows of the LP. 2610 // 2611 // The following operations are done: 2612 // 1. detect implied free variables 2613 // 2. detect implied free constraints 2614 // 3. detect infeasible constraints 2615 // 4. remove unconstrained constraints 2616 // 5. remove empty constraints 2617 // 6. remove row singletons and tighten the corresponding variable bounds if necessary 2618 // 7. remove doubleton equation, aka aggregation 2619 // 8. detect forcing rows and fix the corresponding variables 2620 2621 int remRows = 0; 2622 int remNzos = 0; 2623 int chgLRhs = 0; 2624 int chgBnds = 0; 2625 int keptBnds = 0; 2626 int keptLRhs = 0; 2627 2628 int oldRows = lp.nRows(); 2629 2630 bool redundantLower; 2631 bool redundantUpper; 2632 bool redundantLhs; 2633 bool redundantRhs; 2634 2635 for(int i = lp.nRows() - 1; i >= 0; --i) 2636 { 2637 const SVectorBase<R>& row = lp.rowVector(i); 2638 2639 2640 // compute bounds on constraint value 2641 R lhsBnd = 0.0; // minimal activity (finite summands) 2642 R rhsBnd = 0.0; // maximal activity (finite summands) 2643 int lhsCnt = 0; // number of R(-infinity) summands in minimal activity 2644 int rhsCnt = 0; // number of +R(infinity) summands in maximal activity 2645 2646 for(int k = 0; k < row.size(); ++k) 2647 { 2648 R aij = row.value(k); 2649 int j = row.index(k); 2650 2651 if(!isNotZero(aij, R(1.0 / R(infinity)))) 2652 { 2653 SPX_MSG_WARNING((*this->spxout), (*this->spxout) << "Warning: tiny nonzero coefficient " << aij << 2654 " in row " << i << "\n"); 2655 } 2656 2657 if(aij > 0.0) 2658 { 2659 if(lp.lower(j) <= R(-infinity)) 2660 ++lhsCnt; 2661 else 2662 lhsBnd += aij * lp.lower(j); 2663 2664 if(lp.upper(j) >= R(infinity)) 2665 ++rhsCnt; 2666 else 2667 rhsBnd += aij * lp.upper(j); 2668 } 2669 else if(aij < 0.0) 2670 { 2671 if(lp.lower(j) <= R(-infinity)) 2672 ++rhsCnt; 2673 else 2674 rhsBnd += aij * lp.lower(j); 2675 2676 if(lp.upper(j) >= R(infinity)) 2677 ++lhsCnt; 2678 else 2679 lhsBnd += aij * lp.upper(j); 2680 } 2681 } 2682 2683 #if SOPLEX_FREE_BOUNDS 2684 2685 // 1. detect implied free variables 2686 if(rhsCnt <= 1 || lhsCnt <= 1) 2687 { 2688 for(int k = 0; k < row.size(); ++k) 2689 { 2690 R aij = row.value(k); 2691 int j = row.index(k); 2692 2693 redundantLower = false; 2694 redundantUpper = false; 2695 2696 SOPLEX_ASSERT_WARN("WMAISM12", isNotZero(aij, R(1.0 / R(infinity)))); 2697 2698 if(aij > 0.0) 2699 { 2700 if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity) && rhsCnt <= 1 2701 && NErel(lp.lhs(i), rhsBnd, feastol()) 2702 // do not perform if strongly different orders of magnitude occur 2703 && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > this->tolerances()->epsilon()) 2704 { 2705 R lo = R(-infinity); 2706 R scale = maxAbs(lp.lhs(i), rhsBnd); 2707 2708 if(scale < 1.0) 2709 scale = 1.0; 2710 2711 R z = (lp.lhs(i) / scale) - (rhsBnd / scale); 2712 2713 if(isZero(z, this->epsZero())) 2714 z = 0.0; 2715 2716 assert(rhsCnt > 0 || lp.upper(j) < R(infinity)); 2717 2718 if(rhsCnt == 0) 2719 lo = lp.upper(j) + z * scale / aij; 2720 else if(lp.upper(j) >= R(infinity)) 2721 lo = z * scale / aij; 2722 2723 if(isZero(lo, this->epsZero())) 2724 lo = 0.0; 2725 2726 if(GErel(lo, lp.lower(j), feastol())) 2727 { 2728 SPxOut::debug(this, "IMAISM13 row {}: redundant lower bound on x{} -> lower={} ({})\n", i, j, lo, 2729 lp.lower(j)); 2730 2731 redundantLower = true; 2732 } 2733 2734 } 2735 2736 if(lp.rhs(i) < R(infinity) && lp.upper(j) < R(infinity) && lhsCnt <= 1 2737 && NErel(lp.rhs(i), lhsBnd, feastol()) 2738 // do not perform if strongly different orders of magnitude occur 2739 && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > this->tolerances()->epsilon()) 2740 { 2741 R up = R(infinity); 2742 R scale = maxAbs(lp.rhs(i), lhsBnd); 2743 2744 if(scale < 1.0) 2745 scale = 1.0; 2746 2747 R z = (lp.rhs(i) / scale) - (lhsBnd / scale); 2748 2749 if(isZero(z, this->epsZero())) 2750 z = 0.0; 2751 2752 assert(lhsCnt > 0 || lp.lower(j) > R(-infinity)); 2753 2754 if(lhsCnt == 0) 2755 up = lp.lower(j) + z * scale / aij; 2756 else if(lp.lower(j) <= R(-infinity)) 2757 up = z * scale / aij; 2758 2759 if(isZero(up, this->epsZero())) 2760 up = 0.0; 2761 2762 if(LErel(up, lp.upper(j), feastol())) 2763 { 2764 SPxOut::debug(this, "IMAISM14 row {}: redundant upper bound on x{} -> upper={} ({})\n", i, j, up, 2765 lp.upper(j)); 2766 2767 redundantUpper = true; 2768 } 2769 } 2770 2771 if(redundantLower) 2772 { 2773 // no upper bound on x_j OR redundant upper bound 2774 if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds)) 2775 { 2776 ++lhsCnt; 2777 lhsBnd -= aij * lp.lower(j); 2778 2779 lp.changeLower(j, R(-infinity)); 2780 ++chgBnds; 2781 } 2782 else 2783 ++keptBnds; 2784 } 2785 2786 if(redundantUpper) 2787 { 2788 // no lower bound on x_j OR redundant lower bound 2789 if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds)) 2790 { 2791 ++rhsCnt; 2792 rhsBnd -= aij * lp.upper(j); 2793 2794 lp.changeUpper(j, R(infinity)); 2795 ++chgBnds; 2796 } 2797 else 2798 ++keptBnds; 2799 } 2800 } 2801 else if(aij < 0.0) 2802 { 2803 if(lp.lhs(i) > R(-infinity) && lp.upper(j) < R(infinity) && rhsCnt <= 1 2804 && NErel(lp.lhs(i), rhsBnd, feastol()) 2805 // do not perform if strongly different orders of magnitude occur 2806 && spxAbs(lp.lhs(i) / maxAbs(rhsBnd, R(1.0))) > this->tolerances()->epsilon()) 2807 { 2808 R up = R(infinity); 2809 R scale = maxAbs(lp.lhs(i), rhsBnd); 2810 2811 if(scale < 1.0) 2812 scale = 1.0; 2813 2814 R z = (lp.lhs(i) / scale) - (rhsBnd / scale); 2815 2816 if(isZero(z, this->epsZero())) 2817 z = 0.0; 2818 2819 assert(rhsCnt > 0 || lp.lower(j) > R(-infinity)); 2820 2821 if(rhsCnt == 0) 2822 up = lp.lower(j) + z * scale / aij; 2823 else if(lp.lower(j) <= R(-infinity)) 2824 up = z * scale / aij; 2825 2826 if(isZero(up, this->epsZero())) 2827 up = 0.0; 2828 2829 if(LErel(up, lp.upper(j), feastol())) 2830 { 2831 SPxOut::debug(this, "IMAISM15 row {}: redundant upper bound on x{} -> upper={} ({})\n", i, j, up, 2832 lp.upper(j)); 2833 2834 redundantUpper = true; 2835 } 2836 } 2837 2838 if(lp.rhs(i) < R(infinity) && lp.lower(j) > R(-infinity) && lhsCnt <= 1 2839 && NErel(lp.rhs(i), lhsBnd, feastol()) 2840 // do not perform if strongly different orders of magnitude occur 2841 && spxAbs(lp.rhs(i) / maxAbs(lhsBnd, R(1.0))) > this->tolerances()->epsilon()) 2842 { 2843 R lo = R(-infinity); 2844 R scale = maxAbs(lp.rhs(i), lhsBnd); 2845 2846 if(scale < 1.0) 2847 scale = 1.0; 2848 2849 R z = (lp.rhs(i) / scale) - (lhsBnd / scale); 2850 2851 if(isZero(z, this->epsZero())) 2852 z = 0.0; 2853 2854 assert(lhsCnt > 0 || lp.upper(j) < R(infinity)); 2855 2856 if(lhsCnt == 0) 2857 lo = lp.upper(j) + z * scale / aij; 2858 else if(lp.upper(j) >= R(infinity)) 2859 lo = z * scale / aij; 2860 2861 if(isZero(lo, this->epsZero())) 2862 lo = 0.0; 2863 2864 if(GErel(lo, lp.lower(j), this->tolerances()->epsilon())) 2865 { 2866 SPxOut::debug(this, "IMAISM16 row {}: redundant lower bound on x{} -> lower={} ({})\n", i, j, lo, 2867 lp.lower(j)); 2868 2869 redundantLower = true; 2870 } 2871 } 2872 2873 if(redundantUpper) 2874 { 2875 // no lower bound on x_j OR redundant lower bound 2876 if((lp.lower(j) <= R(-infinity)) || redundantLower || (!m_keepbounds)) 2877 { 2878 ++lhsCnt; 2879 lhsBnd -= aij * lp.upper(j); 2880 2881 lp.changeUpper(j, R(infinity)); 2882 ++chgBnds; 2883 } 2884 else 2885 ++keptBnds; 2886 } 2887 2888 if(redundantLower) 2889 { 2890 // no upper bound on x_j OR redundant upper bound 2891 if((lp.upper(j) >= R(infinity)) || redundantUpper || (!m_keepbounds)) 2892 { 2893 ++rhsCnt; 2894 rhsBnd -= aij * lp.lower(j); 2895 2896 lp.changeLower(j, R(-infinity)); 2897 ++chgBnds; 2898 } 2899 else 2900 ++keptBnds; 2901 } 2902 } 2903 } 2904 } 2905 2906 #endif 2907 2908 #if SOPLEX_FREE_LHS_RHS 2909 2910 redundantLhs = false; 2911 redundantRhs = false; 2912 2913 // 2. detect implied free constraints 2914 if(lp.lhs(i) > R(-infinity) && lhsCnt == 0 && GErel(lhsBnd, lp.lhs(i), feastol())) 2915 { 2916 SPxOut::debug(this, "IMAISM17 row {}: redundant lhs -> lhsBnd={} lhs={}\n", i, lhsBnd, lp.lhs(i)); 2917 2918 redundantLhs = true; 2919 } 2920 2921 if(lp.rhs(i) < R(infinity) && rhsCnt == 0 && LErel(rhsBnd, lp.rhs(i), feastol())) 2922 { 2923 SPxOut::debug(this, "IMAISM18 row {}: redundant rhs -> rhsBnd={} rhs={}\n", i, rhsBnd, lp.rhs(i)); 2924 2925 redundantRhs = true; 2926 } 2927 2928 if(redundantLhs) 2929 { 2930 // no rhs for constraint i OR redundant rhs 2931 if((lp.rhs(i) >= R(infinity)) || redundantRhs || (!m_keepbounds)) 2932 { 2933 lp.changeLhs(i, R(-infinity)); 2934 ++chgLRhs; 2935 } 2936 else 2937 ++keptLRhs; 2938 } 2939 2940 if(redundantRhs) 2941 { 2942 // no lhs for constraint i OR redundant lhs 2943 if((lp.lhs(i) <= R(-infinity)) || redundantLhs || (!m_keepbounds)) 2944 { 2945 lp.changeRhs(i, R(infinity)); 2946 ++chgLRhs; 2947 } 2948 else 2949 ++keptLRhs; 2950 } 2951 2952 #endif 2953 2954 // 3. infeasible constraint 2955 if(LTrel(lp.rhs(i), lp.lhs(i), feastol()) || 2956 (LTrel(rhsBnd, lp.lhs(i), feastol()) && rhsCnt == 0) || 2957 (GTrel(lhsBnd, lp.rhs(i), feastol()) && lhsCnt == 0)) 2958 { 2959 SPxOut::debug(this, "IMAISM19 row {}: infeasible -> lhs={} rhs={} lhsBnd={} rhsBnd={}\n", i, 2960 lp.lhs(i), lp.rhs(i), lhsBnd, rhsBnd); 2961 return this->INFEASIBLE; 2962 } 2963 2964 #if SOPLEX_FREE_CONSTRAINT 2965 2966 // 4. unconstrained constraint 2967 if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity)) 2968 { 2969 SPxOut::debug(this, "IMAISM20 row {}: unconstrained -> removed\n", i); 2970 2971 std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances)); 2972 m_hist.append(ptr); 2973 2974 ++remRows; 2975 remNzos += row.size(); 2976 removeRow(lp, i); 2977 2978 ++m_stat[FREE_ROW]; 2979 2980 continue; 2981 } 2982 2983 #endif 2984 2985 #if SOPLEX_EMPTY_CONSTRAINT 2986 2987 // 5. empty constraint 2988 if(row.size() == 0) 2989 { 2990 SPxOut::debug(this, "IMAISM21 row {}: empty ->", i); 2991 2992 if(LT(lp.rhs(i), R(0.0), feastol()) || GT(lp.lhs(i), R(0.0), feastol())) 2993 { 2994 SPxOut::debug(this, " infeasible lhs={} rhs={}\n", lp.lhs(i), lp.rhs(i)); 2995 return this->INFEASIBLE; 2996 } 2997 2998 SPxOut::debug(this, " removed\n"); 2999 3000 std::shared_ptr<PostStep> ptr(new EmptyConstraintPS(lp, i, this->_tolerances)); 3001 m_hist.append(ptr); 3002 3003 ++remRows; 3004 removeRow(lp, i); 3005 3006 ++m_stat[EMPTY_ROW]; 3007 3008 continue; 3009 } 3010 3011 #endif 3012 3013 #if SOPLEX_ROW_SINGLETON 3014 3015 // 6. row singleton 3016 if(row.size() == 1) 3017 { 3018 removeRowSingleton(lp, row, i); 3019 continue; 3020 } 3021 3022 #endif 3023 3024 #if SOPLEX_AGGREGATE_VARS 3025 3026 // 7. row doubleton, aka. simple aggregation of two variables in an equation 3027 if(row.size() == 2 && EQrel(lp.lhs(i), lp.rhs(i), feastol())) 3028 { 3029 aggregateVars(lp, row, i); 3030 continue; 3031 } 3032 3033 #endif 3034 3035 #if SOPLEX_FORCE_CONSTRAINT 3036 3037 // 8. forcing constraint (postsolving) 3038 // fix variables to obtain the upper bound on constraint value 3039 if(rhsCnt == 0 && EQrel(rhsBnd, lp.lhs(i), feastol())) 3040 { 3041 SPxOut::debug(this, "IMAISM24 row {}: forcing constraint fix on lhs -> lhs={} rhsBnd={}\n", i, 3042 lp.lhs(i), rhsBnd); 3043 3044 DataArray<bool> fixedCol(row.size()); 3045 Array<R> lowers(row.size()); 3046 Array<R> uppers(row.size()); 3047 3048 for(int k = 0; k < row.size(); ++k) 3049 { 3050 R aij = row.value(k); 3051 int j = row.index(k); 3052 3053 fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), this->epsZero())); 3054 3055 lowers[k] = lp.lower(j); 3056 uppers[k] = lp.upper(j); 3057 3058 SOPLEX_ASSERT_WARN("WMAISM25", isNotZero(aij, R(1.0 / R(infinity)))); 3059 3060 if(aij > 0.0) 3061 lp.changeLower(j, lp.upper(j)); 3062 else 3063 lp.changeUpper(j, lp.lower(j)); 3064 } 3065 3066 std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, true, fixedCol, lowers, 3067 uppers, this->_tolerances)); 3068 m_hist.append(ptr); 3069 3070 ++remRows; 3071 remNzos += row.size(); 3072 removeRow(lp, i); 3073 3074 ++m_stat[FORCE_ROW]; 3075 3076 continue; 3077 } 3078 3079 // fix variables to obtain the lower bound on constraint value 3080 if(lhsCnt == 0 && EQrel(lhsBnd, lp.rhs(i), feastol())) 3081 { 3082 SPxOut::debug(this, "IMAISM26 row {}: forcing constraint fix on rhs -> rhs={} lhsBnd={}\n", i, 3083 lp.rhs(i), lhsBnd); 3084 3085 DataArray<bool> fixedCol(row.size()); 3086 Array<R> lowers(row.size()); 3087 Array<R> uppers(row.size()); 3088 3089 for(int k = 0; k < row.size(); ++k) 3090 { 3091 R aij = row.value(k); 3092 int j = row.index(k); 3093 3094 fixedCol[k] = !(EQrel(lp.upper(j), lp.lower(j), this->epsZero())); 3095 3096 lowers[k] = lp.lower(j); 3097 uppers[k] = lp.upper(j); 3098 3099 SOPLEX_ASSERT_WARN("WMAISM27", isNotZero(aij, R(1.0 / R(infinity)))); 3100 3101 if(aij > 0.0) 3102 lp.changeUpper(j, lp.lower(j)); 3103 else 3104 lp.changeLower(j, lp.upper(j)); 3105 } 3106 3107 std::shared_ptr<PostStep> ptr(new ForceConstraintPS(lp, i, false, fixedCol, lowers, 3108 uppers, this->_tolerances)); 3109 m_hist.append(ptr); 3110 3111 ++remRows; 3112 remNzos += row.size(); 3113 removeRow(lp, i); 3114 3115 ++m_stat[FORCE_ROW]; 3116 3117 continue; 3118 } 3119 3120 #endif 3121 } 3122 3123 assert(remRows > 0 || remNzos == 0); 3124 3125 if(remRows + chgLRhs + chgBnds > 0) 3126 { 3127 this->m_remRows += remRows; 3128 this->m_remNzos += remNzos; 3129 this->m_chgLRhs += chgLRhs; 3130 this->m_chgBnds += chgBnds; 3131 this->m_keptBnds += keptBnds; 3132 this->m_keptLRhs += keptLRhs; 3133 3134 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (rows) removed " 3135 << remRows << " rows, " 3136 << remNzos << " non-zeros, " 3137 << chgBnds << " col bounds, " 3138 << chgLRhs << " row bounds; kept " 3139 << keptBnds << " column bounds, " 3140 << keptLRhs << " row bounds" 3141 << std::endl;) 3142 3143 if(remRows > this->m_minReduction * oldRows) 3144 again = true; 3145 } 3146 3147 return this->OKAY; 3148 } 3149 3150 template <class R> 3151 typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyCols(SPxLPBase<R>& lp, bool& again) 3152 { 3153 3154 // This method simplifies the columns of the LP. 3155 // 3156 // The following operations are done: 3157 // 1. detect empty columns and fix corresponding variables 3158 // 2. detect variables that are unconstrained from below or above 3159 // and fix corresponding variables or remove involved constraints 3160 // 3. fix variables 3161 // 4. use column singleton variables with zero objective to adjust constraint bounds 3162 // 5. (not free) column singleton combined with doubleton equation are 3163 // used to make the column singleton variable free 3164 // 6. substitute (implied) free column singletons 3165 3166 int remRows = 0; 3167 int remCols = 0; 3168 int remNzos = 0; 3169 int chgBnds = 0; 3170 3171 int oldCols = lp.nCols(); 3172 int oldRows = lp.nRows(); 3173 3174 for(int j = lp.nCols() - 1; j >= 0; --j) 3175 { 3176 const SVectorBase<R>& col = lp.colVector(j); 3177 3178 // infeasible bounds 3179 if(GTrel(lp.lower(j), lp.upper(j), feastol())) 3180 { 3181 SPxOut::debug(this, "IMAISM29 col {}: infeasible bounds on x{} -> lower={} upper={}\n", j, j, 3182 lp.lower(j), lp.upper(j)); 3183 return this->INFEASIBLE; 3184 } 3185 3186 // 1. empty column 3187 if(col.size() == 0) 3188 { 3189 #if SOPLEX_COMPTY_COLUMN 3190 SPxOut::debug(this, "IMAISM30 col {}: empty -> maxObj={} lower={} upper={}\n", j, lp.maxObj(j), 3191 lp.lower(j), lp.upper(j)); 3192 3193 R val; 3194 3195 if(lp.maxObj(j) > R(0.0)) 3196 { 3197 if(lp.upper(j) >= R(infinity)) 3198 { 3199 SPxOut::debug(this, " unbounded\n"); 3200 return this->UNBOUNDED; 3201 } 3202 3203 val = lp.upper(j); 3204 } 3205 else if(lp.maxObj(j) < R(0.0)) 3206 { 3207 if(lp.lower(j) <= R(-infinity)) 3208 { 3209 SPxOut::debug(this, " unbounded\n"); 3210 return this->UNBOUNDED; 3211 } 3212 3213 val = lp.lower(j); 3214 } 3215 else 3216 { 3217 assert(isZero(lp.maxObj(j), this->epsZero())); 3218 3219 // any value within the bounds is ok 3220 if(lp.lower(j) > R(-infinity)) 3221 val = lp.lower(j); 3222 else if(lp.upper(j) < R(infinity)) 3223 val = lp.upper(j); 3224 else 3225 val = 0.0; 3226 } 3227 3228 SPxOut::debug(this, " removed\n"); 3229 3230 std::shared_ptr<PostStep> ptr1(new FixBoundsPS(lp, j, val, this->_tolerances)); 3231 std::shared_ptr<PostStep> ptr2(new FixVariablePS(lp, *this, j, val, this->_tolerances)); 3232 m_hist.append(ptr1); 3233 m_hist.append(ptr2); 3234 3235 ++remCols; 3236 removeCol(lp, j); 3237 3238 ++m_stat[EMPTY_COL]; 3239 3240 continue; 3241 #endif 3242 } 3243 3244 if(NErel(lp.lower(j), lp.upper(j), feastol())) 3245 { 3246 // will be set to false if any constraint implies a bound on the variable 3247 bool loFree = true; 3248 bool upFree = true; 3249 3250 // 1. fix and remove variables 3251 for(int k = 0; k < col.size(); ++k) 3252 { 3253 if(!loFree && !upFree) 3254 break; 3255 3256 int i = col.index(k); 3257 3258 // warn since this unhandled case may slip through unnoticed otherwise 3259 SOPLEX_ASSERT_WARN("WMAISM31", isNotZero(col.value(k), R(1.0 / R(infinity)))); 3260 3261 if(col.value(k) > 0.0) 3262 { 3263 if(lp.rhs(i) < R(infinity)) 3264 upFree = false; 3265 3266 if(lp.lhs(i) > R(-infinity)) 3267 loFree = false; 3268 } 3269 else if(col.value(k) < 0.0) 3270 { 3271 if(lp.rhs(i) < R(infinity)) 3272 loFree = false; 3273 3274 if(lp.lhs(i) > R(-infinity)) 3275 upFree = false; 3276 } 3277 } 3278 3279 // 2. detect variables that are unconstrained from below or above 3280 // max 3 x 3281 // s.t. 5 x >= 8 3282 if(GT(lp.maxObj(j), R(0.0), this->epsZero()) && upFree) 3283 { 3284 #if SOPLEX_FIX_VARIABLE 3285 SPxOut::debug(this, "IMAISM32 col {}: x{} unconstrained above ->", j, j); 3286 3287 if(lp.upper(j) >= R(infinity)) 3288 { 3289 SPxOut::debug(this, " unbounded\n"); 3290 3291 return this->UNBOUNDED; 3292 } 3293 3294 SPxOut::debug(this, " fixed at upper={}\n", lp.upper(j)); 3295 3296 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances)); 3297 m_hist.append(ptr); 3298 lp.changeLower(j, lp.upper(j)); 3299 } 3300 // max -3 x 3301 // s.t. 5 x <= 8 3302 else if(LT(lp.maxObj(j), R(0.0), this->epsZero()) && loFree) 3303 { 3304 SPxOut::debug(this, "IMAISM33 col {}: x{} unconstrained below ->", j, j); 3305 3306 if(lp.lower(j) <= R(-infinity)) 3307 { 3308 SPxOut::debug(this, " unbounded\n"); 3309 3310 return this->UNBOUNDED; 3311 } 3312 3313 SPxOut::debug(this, " fixed at lower={}\n", lp.lower(j)); 3314 3315 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances)); 3316 m_hist.append(ptr); 3317 lp.changeUpper(j, lp.lower(j)); 3318 #endif 3319 } 3320 else if(isZero(lp.maxObj(j), this->epsZero())) 3321 { 3322 #if SOPLEX_FREE_ZERO_OBJ_VARIABLE 3323 bool unconstrained_below = loFree && lp.lower(j) <= R(-infinity); 3324 bool unconstrained_above = upFree && lp.upper(j) >= R(infinity); 3325 3326 if(unconstrained_below || unconstrained_above) 3327 { 3328 SPxOut::debug(this, "IMAISM34 col {}: x{} unconstrained {} with zero objective ({})\n", j, j, 3329 (unconstrained_below ? "below" : "above"), lp.maxObj(j)); 3330 3331 DSVectorBase<R> col_idx_sorted(col); 3332 3333 // sort col elements by increasing idx 3334 IdxCompare compare; 3335 SPxQuicksort(col_idx_sorted.mem(), col_idx_sorted.size(), compare); 3336 3337 std::shared_ptr<PostStep> ptr(new FreeZeroObjVariablePS(lp, j, unconstrained_below, 3338 col_idx_sorted, this->_tolerances)); 3339 m_hist.append(ptr); 3340 3341 // we have to remove the rows with larger idx first, because otherwise the rows are reorder and indices 3342 // are out-of-date 3343 remRows += col.size(); 3344 3345 for(int k = col_idx_sorted.size() - 1; k >= 0; --k) 3346 removeRow(lp, col_idx_sorted.index(k)); 3347 3348 // remove column 3349 removeCol(lp, j); 3350 3351 // statistics 3352 for(int k = 0; k < col.size(); ++k) 3353 { 3354 int l = col.index(k); 3355 remNzos += lp.rowVector(l).size(); 3356 } 3357 3358 ++m_stat[FREE_ZOBJ_COL]; 3359 ++remCols; 3360 3361 continue; 3362 } 3363 3364 #endif 3365 } 3366 } 3367 3368 #if SOPLEX_FIX_VARIABLE 3369 3370 // 3. fix variable 3371 if(EQrel(lp.lower(j), lp.upper(j), feastol())) 3372 { 3373 SPxOut::debug(this, "IMAISM36 col {}: x{} fixed -> lower={} upper={}\n", j, j, lp.lower(j), 3374 lp.upper(j)); 3375 3376 fixColumn(lp, j); 3377 3378 ++remCols; 3379 remNzos += col.size(); 3380 removeCol(lp, j); 3381 3382 ++m_stat[FIX_COL]; 3383 3384 continue; 3385 } 3386 3387 #endif 3388 3389 // handle column singletons 3390 if(col.size() == 1) 3391 { 3392 R aij = col.value(0); 3393 int i = col.index(0); 3394 3395 // 4. column singleton with zero objective 3396 if(isZero(lp.maxObj(j), this->epsZero())) 3397 { 3398 #if SOPLEX_ZERO_OBJ_COL_SINGLETON 3399 SPxOut::debug(this, "IMAISM37 col {}: singleton in row {} with zero objective", j, i); 3400 3401 R lhs = R(-infinity); 3402 R rhs = +R(infinity); 3403 3404 if(GT(aij, R(0.0), this->epsZero())) 3405 { 3406 if(lp.lhs(i) > R(-infinity) && lp.upper(j) < R(infinity)) 3407 lhs = lp.lhs(i) - aij * lp.upper(j); 3408 3409 if(lp.rhs(i) < R(infinity) && lp.lower(j) > R(-infinity)) 3410 rhs = lp.rhs(i) - aij * lp.lower(j); 3411 } 3412 else if(LT(aij, R(0.0), this->epsZero())) 3413 { 3414 if(lp.lhs(i) > R(-infinity) && lp.lower(j) > R(-infinity)) 3415 lhs = lp.lhs(i) - aij * lp.lower(j); 3416 3417 if(lp.rhs(i) < R(infinity) && lp.upper(j) < R(infinity)) 3418 rhs = lp.rhs(i) - aij * lp.upper(j); 3419 } 3420 else 3421 { 3422 lhs = lp.lhs(i); 3423 rhs = lp.rhs(i); 3424 } 3425 3426 if(isZero(lhs, this->epsZero())) 3427 lhs = 0.0; 3428 3429 if(isZero(rhs, this->epsZero())) 3430 rhs = 0.0; 3431 3432 SPxOut::debug(this, " removed -> lhs={} ({}) rhs={} ({})\n", lhs, lp.lhs(i), rhs, lp.rhs(i)); 3433 3434 std::shared_ptr<PostStep> ptr(new ZeroObjColSingletonPS(lp, *this, j, i, this->_tolerances)); 3435 m_hist.append(ptr); 3436 3437 lp.changeRange(i, lhs, rhs); 3438 3439 ++remCols; 3440 ++remNzos; 3441 removeCol(lp, j); 3442 3443 ++m_stat[ZOBJ_SINGLETON_COL]; 3444 3445 if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity)) 3446 { 3447 std::shared_ptr<PostStep> ptr2(new FreeConstraintPS(lp, i, this->_tolerances)); 3448 m_hist.append(ptr2); 3449 3450 ++remRows; 3451 removeRow(lp, i); 3452 3453 ++m_stat[FREE_ROW]; 3454 } 3455 3456 continue; 3457 #endif 3458 } 3459 3460 // 5. not free column singleton combined with doubleton equation 3461 else if(EQrel(lp.lhs(i), lp.rhs(i), feastol()) && 3462 lp.rowVector(i).size() == 2 && 3463 (lp.lower(j) > R(-infinity) || lp.upper(j) < R(infinity))) 3464 { 3465 #if SOPLEX_DOUBLETON_EQUATION 3466 SPxOut::debug(this, "IMAISM38 col {}: singleton in row {} with doubleton equation ->", j, i); 3467 3468 R lhs = lp.lhs(i); 3469 3470 const SVectorBase<R>& row = lp.rowVector(i); 3471 3472 R aik; 3473 int k; 3474 3475 if(row.index(0) == j) 3476 { 3477 aik = row.value(1); 3478 k = row.index(1); 3479 } 3480 else if(row.index(1) == j) 3481 { 3482 aik = row.value(0); 3483 k = row.index(0); 3484 } 3485 else 3486 throw SPxInternalCodeException("XMAISM11 This should never happen."); 3487 3488 SOPLEX_ASSERT_WARN("WMAISM39", isNotZero(aik, R(1.0 / R(infinity)))); 3489 3490 R lo, up; 3491 R oldLower = lp.lower(k); 3492 R oldUpper = lp.upper(k); 3493 3494 R scale1 = maxAbs(lhs, aij * lp.upper(j)); 3495 R scale2 = maxAbs(lhs, aij * lp.lower(j)); 3496 3497 if(scale1 < 1.0) 3498 scale1 = 1.0; 3499 3500 if(scale2 < 1.0) 3501 scale2 = 1.0; 3502 3503 R z1 = (lhs / scale1) - (aij * lp.upper(j) / scale1); 3504 R z2 = (lhs / scale2) - (aij * lp.lower(j) / scale2); 3505 3506 if(isZero(z1, this->epsZero())) 3507 z1 = 0.0; 3508 3509 if(isZero(z2, this->epsZero())) 3510 z2 = 0.0; 3511 3512 if(aij * aik > 0.0) 3513 { 3514 lo = (lp.upper(j) >= R(infinity)) ? R(-infinity) : z1 * scale1 / aik; 3515 up = (lp.lower(j) <= R(-infinity)) ? R(infinity) : z2 * scale2 / aik; 3516 } 3517 else if(aij * aik < 0.0) 3518 { 3519 lo = (lp.lower(j) <= R(-infinity)) ? R(-infinity) : z2 * scale2 / aik; 3520 up = (lp.upper(j) >= R(infinity)) ? R(infinity) : z1 * scale1 / aik; 3521 } 3522 else 3523 throw SPxInternalCodeException("XMAISM12 This should never happen."); 3524 3525 if(GTrel(lo, lp.lower(k), this->epsZero())) 3526 lp.changeLower(k, lo); 3527 3528 if(LTrel(up, lp.upper(k), this->epsZero())) 3529 lp.changeUpper(k, up); 3530 3531 SPxOut::debug(this, " made free, bounds on x{}: lower={} ({}) upper={} ({})\n", k, lp.lower(k), 3532 oldLower, lp.upper(k), oldUpper); 3533 3534 // infeasible bounds 3535 if(GTrel(lp.lower(k), lp.upper(k), feastol())) 3536 { 3537 SPxOut::debug(this, "new bounds are infeasible\n"); 3538 return this->INFEASIBLE; 3539 } 3540 3541 std::shared_ptr<PostStep> ptr(new DoubletonEquationPS(lp, j, k, i, oldLower, oldUpper, 3542 this->_tolerances)); 3543 m_hist.append(ptr); 3544 3545 if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity)) 3546 chgBnds += 2; 3547 else 3548 ++chgBnds; 3549 3550 lp.changeBounds(j, R(-infinity), R(infinity)); 3551 3552 ++m_stat[DOUBLETON_ROW]; 3553 #endif 3554 } 3555 3556 // 6. (implied) free column singleton 3557 if(lp.lower(j) <= R(-infinity) && lp.upper(j) >= R(infinity)) 3558 { 3559 #if SOPLEX_FREE_COL_SINGLETON 3560 R slackVal = lp.lhs(i); 3561 3562 // constraint i is an inequality constraint -> transform into equation type 3563 if(NErel(lp.lhs(i), lp.rhs(i), feastol())) 3564 { 3565 SPxOut::debug(this, "IMAISM40 col {}: free singleton in inequality constraint\n", j); 3566 3567 // do nothing if constraint i is unconstrained 3568 if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity)) 3569 continue; 3570 3571 // introduce slack variable to obtain equality constraint 3572 R sMaxObj = lp.maxObj(j) / aij; // after substituting variable j in objective 3573 R sLo = lp.lhs(i); 3574 R sUp = lp.rhs(i); 3575 3576 if(GT(sMaxObj, R(0.0), this->epsZero())) 3577 { 3578 if(sUp >= R(infinity)) 3579 { 3580 SPxOut::debug(this, " -> problem unbounded\n"); 3581 return this->UNBOUNDED; 3582 } 3583 3584 slackVal = sUp; 3585 } 3586 else if(LT(sMaxObj, R(0.0), this->epsZero())) 3587 { 3588 if(sLo <= R(-infinity)) 3589 { 3590 SPxOut::debug(this, " -> problem unbounded\n"); 3591 return this->UNBOUNDED; 3592 } 3593 3594 slackVal = sLo; 3595 } 3596 else 3597 { 3598 assert(isZero(sMaxObj, this->epsZero())); 3599 3600 // any value within the bounds is ok 3601 if(sLo > R(-infinity)) 3602 slackVal = sLo; 3603 else if(sUp < R(infinity)) 3604 slackVal = sUp; 3605 else 3606 throw SPxInternalCodeException("XMAISM13 This should never happen."); 3607 } 3608 } 3609 3610 std::shared_ptr<PostStep> ptr(new FreeColSingletonPS(lp, *this, j, i, slackVal, this->_tolerances)); 3611 m_hist.append(ptr); 3612 3613 SPxOut::debug(this, "IMAISM41 col {}: free singleton removed\n", j); 3614 3615 const SVectorBase<R>& row = lp.rowVector(i); 3616 3617 for(int h = 0; h < row.size(); ++h) 3618 { 3619 int k = row.index(h); 3620 3621 if(k != j) 3622 { 3623 R new_obj = lp.obj(k) - (lp.obj(j) * row.value(h) / aij); 3624 lp.changeObj(k, new_obj); 3625 } 3626 } 3627 3628 ++remRows; 3629 ++remCols; 3630 remNzos += row.size(); 3631 removeRow(lp, i); 3632 removeCol(lp, j); 3633 3634 ++m_stat[FREE_SINGLETON_COL]; 3635 3636 continue; 3637 #endif 3638 } 3639 } 3640 } 3641 3642 if(remCols + remRows > 0) 3643 { 3644 this->m_remRows += remRows; 3645 this->m_remCols += remCols; 3646 this->m_remNzos += remNzos; 3647 this->m_chgBnds += chgBnds; 3648 3649 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (columns) removed " 3650 << remRows << " rows, " 3651 << remCols << " cols, " 3652 << remNzos << " non-zeros, " 3653 << chgBnds << " col bounds" 3654 << std::endl;) 3655 3656 if(remCols + remRows > this->m_minReduction * (oldCols + oldRows)) 3657 again = true; 3658 } 3659 3660 return this->OKAY; 3661 } 3662 3663 template <class R> 3664 typename SPxSimplifier<R>::Result SPxMainSM<R>::simplifyDual(SPxLPBase<R>& lp, bool& again) 3665 { 3666 3667 // This method simplifies LP using the following dual structures: 3668 // 3669 // 1. dominated columns 3670 // 2. weakly dominated columns 3671 // 3672 // For constructing the dual variables, it is assumed that the objective sense is max 3673 3674 int remRows = 0; 3675 int remCols = 0; 3676 int remNzos = 0; 3677 3678 int oldRows = lp.nRows(); 3679 int oldCols = lp.nCols(); 3680 3681 DataArray<bool> colSingleton(lp.nCols()); 3682 VectorBase<R> dualVarLo(lp.nRows()); 3683 VectorBase<R> dualVarUp(lp.nRows()); 3684 VectorBase<R> dualConsLo(lp.nCols()); 3685 VectorBase<R> dualConsUp(lp.nCols()); 3686 3687 // init 3688 for(int i = lp.nRows() - 1; i >= 0; --i) 3689 { 3690 // check for unconstrained constraints 3691 if(lp.lhs(i) <= R(-infinity) && lp.rhs(i) >= R(infinity)) 3692 { 3693 SPxOut::debug(this, "IMAISM43 row {}: unconstrained\n", i); 3694 3695 std::shared_ptr<PostStep> ptr(new FreeConstraintPS(lp, i, this->_tolerances)); 3696 m_hist.append(ptr); 3697 3698 ++remRows; 3699 remNzos += lp.rowVector(i).size(); 3700 removeRow(lp, i); 3701 3702 ++m_stat[FREE_ROW]; 3703 3704 continue; 3705 } 3706 3707 // corresponds to maximization sense 3708 dualVarLo[i] = (lp.lhs(i) <= R(-infinity)) ? 0.0 : R(-infinity); 3709 dualVarUp[i] = (lp.rhs(i) >= R(infinity)) ? 0.0 : R(infinity); 3710 } 3711 3712 // compute bounds on the dual variables using column singletons 3713 for(int j = 0; j < lp.nCols(); ++j) 3714 { 3715 if(lp.colVector(j).size() == 1) 3716 { 3717 int i = lp.colVector(j).index(0); 3718 R aij = lp.colVector(j).value(0); 3719 3720 SOPLEX_ASSERT_WARN("WMAISM44", isNotZero(aij, R(1.0 / R(infinity)))); 3721 3722 R bound = lp.maxObj(j) / aij; 3723 3724 if(aij > 0) 3725 { 3726 if(lp.lower(j) <= R(-infinity) && bound < dualVarUp[i]) 3727 dualVarUp[i] = bound; 3728 3729 if(lp.upper(j) >= R(infinity) && bound > dualVarLo[i]) 3730 dualVarLo[i] = bound; 3731 } 3732 else if(aij < 0) 3733 { 3734 if(lp.lower(j) <= R(-infinity) && bound > dualVarLo[i]) 3735 dualVarLo[i] = bound; 3736 3737 if(lp.upper(j) >= R(infinity) && bound < dualVarUp[i]) 3738 dualVarUp[i] = bound; 3739 } 3740 } 3741 3742 } 3743 3744 // compute bounds on the dual constraints 3745 for(int j = 0; j < lp.nCols(); ++j) 3746 { 3747 dualConsLo[j] = dualConsUp[j] = 0.0; 3748 3749 const SVectorBase<R>& col = lp.colVector(j); 3750 3751 for(int k = 0; k < col.size(); ++k) 3752 { 3753 if(dualConsLo[j] <= R(-infinity) && dualConsUp[j] >= R(infinity)) 3754 break; 3755 3756 R aij = col.value(k); 3757 int i = col.index(k); 3758 3759 SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity)))); 3760 3761 if(aij > 0) 3762 { 3763 if(dualVarLo[i] <= R(-infinity)) 3764 dualConsLo[j] = R(-infinity); 3765 else 3766 dualConsLo[j] += aij * dualVarLo[i]; 3767 3768 if(dualVarUp[i] >= R(infinity)) 3769 dualConsUp[j] = R(infinity); 3770 else 3771 dualConsUp[j] += aij * dualVarUp[i]; 3772 } 3773 else if(aij < 0) 3774 { 3775 if(dualVarLo[i] <= R(-infinity)) 3776 dualConsUp[j] = R(infinity); 3777 else 3778 dualConsUp[j] += aij * dualVarLo[i]; 3779 3780 if(dualVarUp[i] >= R(infinity)) 3781 dualConsLo[j] = R(-infinity); 3782 else 3783 dualConsLo[j] += aij * dualVarUp[i]; 3784 } 3785 } 3786 } 3787 3788 for(int j = lp.nCols() - 1; j >= 0; --j) 3789 { 3790 if(lp.colVector(j).size() <= 1) 3791 continue; 3792 3793 // dual infeasibility checks 3794 if(LTrel(dualConsUp[j], dualConsLo[j], opttol())) 3795 { 3796 SPxOut::debug(this, "IMAISM46 col {}: dual infeasible -> dual lhs bound={} dual rhs bound={}\n", j, 3797 dualConsLo[j], dualConsUp[j]); 3798 return this->DUAL_INFEASIBLE; 3799 } 3800 3801 R obj = lp.maxObj(j); 3802 3803 // 1. dominated column 3804 // Is the problem really unbounded in the cases below ??? Or is only dual infeasibility be shown 3805 if(GTrel(obj, dualConsUp[j], opttol())) 3806 { 3807 #if SOPLEX_DOMINATED_COLUMN 3808 SPxOut::debug(this, "IMAISM47 col{}: dominated -> maxObj={} dual rhs bound={}\n", j, obj, 3809 dualConsUp[j]); 3810 3811 if(lp.upper(j) >= R(infinity)) 3812 { 3813 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;) 3814 return this->UNBOUNDED; 3815 } 3816 3817 SPxOut::debug(this, " fixed at upper={}\n", lp.upper(j)); 3818 3819 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances)); 3820 m_hist.append(ptr); 3821 lp.changeLower(j, lp.upper(j)); 3822 3823 ++m_stat[DOMINATED_COL]; 3824 #endif 3825 } 3826 else if(LTrel(obj, dualConsLo[j], opttol())) 3827 { 3828 #if SOPLEX_DOMINATED_COLUMN 3829 SPxOut::debug(this, "IMAISM48 col {}: dominated -> maxObj={} dual lhs bound={}\n", j, obj, 3830 dualConsLo[j]); 3831 3832 if(lp.lower(j) <= R(-infinity)) 3833 { 3834 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << " unbounded" << std::endl;) 3835 return this->UNBOUNDED; 3836 } 3837 3838 SPxOut::debug(this, " fixed at lower={}\n", lp.lower(j)); 3839 3840 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances)); 3841 m_hist.append(ptr); 3842 lp.changeUpper(j, lp.lower(j)); 3843 3844 ++m_stat[DOMINATED_COL]; 3845 #endif 3846 } 3847 3848 // 2. weakly dominated column (no postsolving) 3849 else if(lp.upper(j) < R(infinity) && EQrel(obj, dualConsUp[j], opttol())) 3850 { 3851 #if SOPLEX_WEAKLY_DOMINATED_COLUMN 3852 SPxOut::debug(this, "IMAISM49 col {}: weakly dominated -> maxObj={} dual rhs bound={}\n", j, obj, 3853 dualConsUp[j]); 3854 3855 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.upper(j), this->_tolerances)); 3856 m_hist.append(ptr); 3857 lp.changeLower(j, lp.upper(j)); 3858 3859 ++m_stat[WEAKLY_DOMINATED_COL]; 3860 #endif 3861 } 3862 else if(lp.lower(j) > R(-infinity) && EQrel(obj, dualConsLo[j], opttol())) 3863 { 3864 #if SOPLEX_WEAKLY_DOMINATED_COLUMN 3865 SPxOut::debug(this, "IMAISM50 col {}: weakly dominated -> maxObj={} dual lhs bound={}\n", j, obj, 3866 dualConsLo[j]); 3867 3868 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j, lp.lower(j), this->_tolerances)); 3869 m_hist.append(ptr); 3870 lp.changeUpper(j, lp.lower(j)); 3871 3872 ++m_stat[WEAKLY_DOMINATED_COL]; 3873 #endif 3874 } 3875 3876 // fix column 3877 if(EQrel(lp.lower(j), lp.upper(j), feastol())) 3878 { 3879 #if SOPLEX_FIX_VARIABLE 3880 fixColumn(lp, j); 3881 3882 ++remCols; 3883 remNzos += lp.colVector(j).size(); 3884 removeCol(lp, j); 3885 3886 ++m_stat[FIX_COL]; 3887 #endif 3888 } 3889 } 3890 3891 3892 assert(remRows > 0 || remCols > 0 || remNzos == 0); 3893 3894 if(remCols + remRows > 0) 3895 { 3896 this->m_remRows += remRows; 3897 this->m_remCols += remCols; 3898 this->m_remNzos += remNzos; 3899 3900 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (dual) removed " 3901 << remRows << " rows, " 3902 << remCols << " cols, " 3903 << remNzos << " non-zeros" 3904 << std::endl;) 3905 3906 if(remCols + remRows > this->m_minReduction * (oldCols + oldRows)) 3907 again = true; 3908 } 3909 3910 return this->OKAY; 3911 } 3912 3913 3914 3915 template <class R> 3916 typename SPxSimplifier<R>::Result SPxMainSM<R>::multiaggregation(SPxLPBase<R>& lp, bool& again) 3917 { 3918 // this simplifier eliminates rows and columns by performing multi aggregations as identified by the constraint 3919 // activities. 3920 int remRows = 0; 3921 int remCols = 0; 3922 int remNzos = 0; 3923 3924 int oldRows = lp.nRows(); 3925 int oldCols = lp.nCols(); 3926 3927 VectorBase<R> upLocks(lp.nCols()); 3928 VectorBase<R> downLocks(lp.nCols()); 3929 3930 for(int j = lp.nCols() - 1; j >= 0; --j) 3931 { 3932 // setting the locks on the variables 3933 upLocks[j] = 0; 3934 downLocks[j] = 0; 3935 3936 if(lp.colVector(j).size() <= 1) 3937 continue; 3938 3939 const SVectorBase<R>& col = lp.colVector(j); 3940 3941 for(int k = 0; k < col.size(); ++k) 3942 { 3943 R aij = col.value(k); 3944 3945 SOPLEX_ASSERT_WARN("WMAISM45", isNotZero(aij, R(1.0 / R(infinity)))); 3946 3947 if(lp.lhs(col.index(k)) > R(-infinity) && lp.rhs(col.index(k)) < R(infinity)) 3948 { 3949 upLocks[j]++; 3950 downLocks[j]++; 3951 } 3952 else if(lp.lhs(col.index(k)) > R(-infinity)) 3953 { 3954 if(aij > 0) 3955 downLocks[j]++; 3956 else if(aij < 0) 3957 upLocks[j]++; 3958 } 3959 else if(lp.rhs(col.index(k)) < R(infinity)) 3960 { 3961 if(aij > 0) 3962 upLocks[j]++; 3963 else if(aij < 0) 3964 downLocks[j]++; 3965 } 3966 } 3967 3968 // multi-aggregate column 3969 if(upLocks[j] == 1 || downLocks[j] == 1) 3970 { 3971 R lower = lp.lower(j); 3972 R upper = lp.upper(j); 3973 int maxOtherLocks; 3974 int bestpos = -1; 3975 bool bestislhs = true; 3976 3977 for(int k = 0; k < col.size(); ++k) 3978 { 3979 int rowNumber; 3980 R lhs; 3981 R rhs; 3982 bool lhsExists; 3983 bool rhsExists; 3984 bool aggLhs; 3985 bool aggRhs; 3986 3987 R val = col.value(k); 3988 3989 rowNumber = col.index(k); 3990 lhs = lp.lhs(rowNumber); 3991 rhs = lp.rhs(rowNumber); 3992 3993 if(EQ(lhs, rhs, feastol())) 3994 continue; 3995 3996 lhsExists = lhs > R(-infinity); 3997 rhsExists = rhs < R(infinity); 3998 3999 if(lp.rowVector(rowNumber).size() <= 2) 4000 maxOtherLocks = INT_MAX; 4001 else if(lp.rowVector(rowNumber).size() == 3) 4002 maxOtherLocks = 3; 4003 else if(lp.rowVector(rowNumber).size() == 4) 4004 maxOtherLocks = 2; 4005 else 4006 maxOtherLocks = 1; 4007 4008 aggLhs = lhsExists 4009 && ((col.value(k) > 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks) 4010 || (col.value(k) < 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks)); 4011 aggRhs = rhsExists 4012 && ((col.value(k) > 0.0 && lp.maxObj(j) >= 0.0 && upLocks[j] == 1 && downLocks[j] <= maxOtherLocks) 4013 || (col.value(k) < 0.0 && lp.maxObj(j) <= 0.0 && downLocks[j] == 1 && upLocks[j] <= maxOtherLocks)); 4014 4015 if(aggLhs || aggRhs) 4016 { 4017 R minRes = 0; // this is the minimum value that the aggregation can attain 4018 R maxRes = 0; // this is the maximum value that the aggregation can attain 4019 4020 // computing the minimum and maximum residuals if variable j is set to zero. 4021 computeMinMaxResidualActivity(lp, rowNumber, j, minRes, maxRes); 4022 4023 // we will try to aggregate to the lhs 4024 if(aggLhs) 4025 { 4026 R minVal; 4027 R maxVal; 4028 4029 // computing the values of the upper and lower bounds for the aggregated variables 4030 computeMinMaxValues(lp, lhs, val, minRes, maxRes, minVal, maxVal); 4031 4032 assert(LE(minVal, maxVal, this->tolerances()->epsilon())); 4033 4034 // if the bounds of the aggregation and the original variable are equivalent, then we can reduce 4035 if((minVal > R(-infinity) && GT(minVal, lower, feastol())) 4036 && (maxVal < R(infinity) && LT(maxVal, upper, feastol()))) 4037 { 4038 bestpos = col.index(k); 4039 bestislhs = true; 4040 break; 4041 } 4042 } 4043 4044 // we will try to aggregate to the rhs 4045 if(aggRhs) 4046 { 4047 R minVal; 4048 R maxVal; 4049 4050 // computing the values of the upper and lower bounds for the aggregated variables 4051 computeMinMaxValues(lp, rhs, val, minRes, maxRes, minVal, maxVal); 4052 4053 assert(LE(minVal, maxVal, this->tolerances()->epsilon())); 4054 4055 if((minVal > R(-infinity) && GT(minVal, lower, feastol())) 4056 && (maxVal < R(infinity) && LT(maxVal, upper, feastol()))) 4057 { 4058 bestpos = col.index(k); 4059 bestislhs = false; 4060 break; 4061 } 4062 } 4063 } 4064 } 4065 4066 // it is only possible to aggregate if a best position has been found 4067 if(bestpos >= 0) 4068 { 4069 const SVectorBase<R>& bestRow = lp.rowVector(bestpos); 4070 // aggregating the variable and applying the fixings to the all other constraints 4071 R aggConstant = (bestislhs ? lp.lhs(bestpos) : lp.rhs( 4072 bestpos)); // this is the lhs or rhs of the aggregated row 4073 R aggAij = 4074 bestRow[j]; // this is the coefficient of the deleted col 4075 4076 SPxOut::debug(this, 4077 "IMAISM51 col {}: Aggregating row: {} Aggregation Constant={} Coefficient of aggregated col={}\n", 4078 j, bestpos, aggConstant, aggAij); 4079 4080 std::shared_ptr<PostStep> ptr(new MultiAggregationPS(lp, *this, bestpos, j, 4081 aggConstant, this->_tolerances)); 4082 m_hist.append(ptr); 4083 4084 for(int k = 0; k < col.size(); ++k) 4085 { 4086 if(col.index(k) != bestpos) 4087 { 4088 int rowNumber = col.index(k); 4089 VectorBase<R> updateRow(lp.nCols()); 4090 R updateRhs = lp.rhs(col.index(k)); 4091 R updateLhs = lp.lhs(col.index(k)); 4092 4093 updateRow = lp.rowVector(col.index(k)); 4094 4095 // updating the row with the best row 4096 for(int l = 0; l < bestRow.size(); l++) 4097 { 4098 if(bestRow.index(l) != j) 4099 { 4100 if(lp.rowVector(rowNumber).pos(bestRow.index(l)) >= 0) 4101 lp.changeElement(rowNumber, bestRow.index(l), updateRow[bestRow.index(l)] 4102 - updateRow[j]*bestRow.value(l) / aggAij); 4103 else 4104 lp.changeElement(rowNumber, bestRow.index(l), -1.0 * updateRow[j]*bestRow.value(l) / aggAij); 4105 } 4106 } 4107 4108 // NOTE: I don't know whether we should change the LHS and RHS if they are currently at R(infinity) 4109 if(lp.rhs(rowNumber) < R(infinity)) 4110 lp.changeRhs(rowNumber, updateRhs - updateRow[j]*aggConstant / aggAij); 4111 4112 if(lp.lhs(rowNumber) > R(-infinity)) 4113 lp.changeLhs(rowNumber, updateLhs - updateRow[j]*aggConstant / aggAij); 4114 4115 assert(LE(lp.lhs(rowNumber), lp.rhs(rowNumber), this->tolerances()->epsilon())); 4116 } 4117 } 4118 4119 for(int l = 0; l < bestRow.size(); l++) 4120 { 4121 if(bestRow.index(l) != j) 4122 lp.changeMaxObj(bestRow.index(l), 4123 lp.maxObj(bestRow.index(l)) - lp.maxObj(j)*bestRow.value(l) / aggAij); 4124 } 4125 4126 ++remCols; 4127 remNzos += lp.colVector(j).size(); 4128 removeCol(lp, j); 4129 ++remRows; 4130 remNzos += lp.rowVector(bestpos).size(); 4131 removeRow(lp, bestpos); 4132 4133 ++m_stat[MULTI_AGG]; 4134 } 4135 } 4136 } 4137 4138 4139 assert(remRows > 0 || remCols > 0 || remNzos == 0); 4140 4141 if(remCols + remRows > 0) 4142 { 4143 this->m_remRows += remRows; 4144 this->m_remCols += remCols; 4145 this->m_remNzos += remNzos; 4146 4147 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (multi-aggregation) removed " 4148 << remRows << " rows, " 4149 << remCols << " cols, " 4150 << remNzos << " non-zeros" 4151 << std::endl;) 4152 4153 if(remCols + remRows > this->m_minReduction * (oldCols + oldRows)) 4154 again = true; 4155 } 4156 4157 return this->OKAY; 4158 } 4159 4160 4161 4162 template <class R> 4163 typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateRows(SPxLPBase<R>& lp, bool& again) 4164 { 4165 4166 // This method simplifies the LP by removing duplicate rows 4167 // Duplicates are detected using the algorithm of Bixby and Wagner [1987] 4168 4169 // Possible extension: use generalized definition of duplicate rows according to Andersen and Andersen 4170 // However: the resulting sparsification is often very small since the involved rows are usually very sparse 4171 4172 int remRows = 0; 4173 int remNzos = 0; 4174 4175 int oldRows = lp.nRows(); 4176 4177 // remove empty rows and columns 4178 typename SPxSimplifier<R>::Result ret = removeEmpty(lp); 4179 4180 if(ret != this->OKAY) 4181 return ret; 4182 4183 #if SOPLEX_ROW_SINGLETON 4184 int rs_remRows = 0; 4185 4186 for(int i = 0; i < lp.nRows(); ++i) 4187 { 4188 const SVectorBase<R>& row = lp.rowVector(i); 4189 4190 if(row.size() == 1) 4191 { 4192 removeRowSingleton(lp, row, i); 4193 rs_remRows++; 4194 } 4195 } 4196 4197 if(rs_remRows > 0) 4198 { 4199 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << 4200 "Simplifier duplicate rows (row singleton stage) removed " 4201 << rs_remRows << " rows, " 4202 << rs_remRows << " non-zeros" 4203 << std::endl;) 4204 } 4205 4206 #endif 4207 4208 if(lp.nRows() < 2) 4209 return this->OKAY; 4210 4211 DataArray<int> pClass(lp.nRows()); // class of parallel rows 4212 DataArray<int> classSize(lp.nRows()); // size of each class 4213 Array<R> scale(lp.nRows()); // scaling factor for each row 4214 int* idxMem = 0; 4215 4216 try 4217 { 4218 spx_alloc(idxMem, lp.nRows()); 4219 } 4220 catch(const SPxMemoryException& x) 4221 { 4222 spx_free(idxMem); 4223 throw x; 4224 } 4225 4226 IdxSet idxSet(lp.nRows(), idxMem); // set of feasible indices for new pClass 4227 4228 // init 4229 pClass[0] = 0; 4230 scale[0] = 0.0; 4231 classSize[0] = lp.nRows(); 4232 4233 for(int i = 1; i < lp.nRows(); ++i) 4234 { 4235 pClass[i] = 0; 4236 scale[i] = 0.0; 4237 classSize[i] = 0; 4238 idxSet.addIdx(i); 4239 } 4240 4241 R oldVal = 0.0; 4242 4243 // main loop 4244 for(int j = 0; j < lp.nCols(); ++j) 4245 { 4246 const SVectorBase<R>& col = lp.colVector(j); 4247 4248 for(int k = 0; k < col.size(); ++k) 4249 { 4250 R aij = col.value(k); 4251 int i = col.index(k); 4252 4253 if(scale[i] == 0.0) 4254 scale[i] = aij; 4255 4256 m_classSetRows[pClass[i]].add(i, aij / scale[i]); 4257 4258 if(--classSize[pClass[i]] == 0) 4259 idxSet.addIdx(pClass[i]); 4260 } 4261 4262 // update each parallel class with non-zero column entry 4263 for(int m = 0; m < col.size(); ++m) 4264 { 4265 int k = pClass[col.index(m)]; 4266 4267 if(m_classSetRows[k].size() > 0) 4268 { 4269 // sort classSet[k] w.r.t. scaled column values 4270 ElementCompare compare(this->tolerances()->epsilon()); 4271 4272 if(m_classSetRows[k].size() > 1) 4273 SPxQuicksort(m_classSetRows[k].mem(), m_classSetRows[k].size(), compare); 4274 4275 // use new index first 4276 int classIdx = idxSet.index(0); 4277 idxSet.remove(0); 4278 4279 for(int l = 0; l < m_classSetRows[k].size(); ++l) 4280 { 4281 if(l != 0 && NErel(m_classSetRows[k].value(l), oldVal, this->epsZero())) 4282 { 4283 classIdx = idxSet.index(0); 4284 idxSet.remove(0); 4285 } 4286 4287 pClass[m_classSetRows[k].index(l)] = classIdx; 4288 ++classSize[classIdx]; 4289 4290 oldVal = m_classSetRows[k].value(l); 4291 } 4292 4293 m_classSetRows[k].clear(); 4294 } 4295 } 4296 } 4297 4298 spx_free(idxMem); 4299 4300 DataArray<bool> remRow(lp.nRows()); 4301 4302 for(int k = 0; k < lp.nRows(); ++k) 4303 m_dupRows[k].clear(); 4304 4305 for(int k = 0; k < lp.nRows(); ++k) 4306 { 4307 remRow[k] = false; 4308 m_dupRows[pClass[k]].add(k, 0.0); 4309 } 4310 4311 const int nRowsOld_tmp = lp.nRows(); 4312 int* perm_tmp = 0; 4313 spx_alloc(perm_tmp, nRowsOld_tmp); 4314 4315 for(int j = 0; j < nRowsOld_tmp; ++j) 4316 { 4317 perm_tmp[j] = 0; 4318 } 4319 4320 int idxFirstDupRows = -1; 4321 int idxLastDupRows = -1; 4322 int numDelRows = 0; 4323 4324 for(int k = 0; k < lp.nRows(); ++k) 4325 { 4326 if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1)) 4327 { 4328 idxLastDupRows = k; 4329 4330 if(idxFirstDupRows < 0) 4331 { 4332 idxFirstDupRows = k; 4333 } 4334 4335 for(int l = 1; l < m_dupRows[k].size(); ++l) 4336 { 4337 int i = m_dupRows[k].index(l); 4338 perm_tmp[i] = -1; 4339 } 4340 4341 numDelRows += (m_dupRows[k].size() - 1); 4342 } 4343 } 4344 4345 { 4346 int k_tmp, j_tmp = -1; 4347 4348 for(k_tmp = j_tmp = 0; k_tmp < nRowsOld_tmp; ++k_tmp) 4349 { 4350 if(perm_tmp[k_tmp] >= 0) 4351 perm_tmp[k_tmp] = j_tmp++; 4352 } 4353 } 4354 4355 // store rhs and lhs changes for combined update 4356 bool doChangeRanges = false; 4357 VectorBase<R> newLhsVec(lp.lhs()); 4358 VectorBase<R> newRhsVec(lp.rhs()); 4359 4360 for(int k = 0; k < lp.nRows(); ++k) 4361 { 4362 if(m_dupRows[k].size() > 1 && !(lp.rowVector(m_dupRows[k].index(0)).size() == 1)) 4363 { 4364 SPxOut::debug(this, "IMAISM53 {} duplicate rows found\n", m_dupRows[k].size()); 4365 4366 m_stat[DUPLICATE_ROW] += m_dupRows[k].size() - 1; 4367 4368 // index of one non-column singleton row in dupRows[k] 4369 int rowIdx = -1; 4370 int maxLhsIdx = -1; 4371 int minRhsIdx = -1; 4372 R maxLhs = R(-infinity); 4373 R minRhs = +R(infinity); 4374 4375 DataArray<bool> isLhsEqualRhs(m_dupRows[k].size()); 4376 4377 // determine strictest bounds on constraint 4378 for(int l = 0; l < m_dupRows[k].size(); ++l) 4379 { 4380 int i = m_dupRows[k].index(l); 4381 isLhsEqualRhs[l] = (lp.lhs(i) == lp.rhs(i)); 4382 4383 SOPLEX_ASSERT_WARN("WMAISM54", isNotZero(scale[i], R(1.0 / R(infinity)))); 4384 4385 if(rowIdx == -1) 4386 { 4387 rowIdx = i; 4388 maxLhs = lp.lhs(rowIdx); 4389 minRhs = lp.rhs(rowIdx); 4390 } 4391 else 4392 { 4393 R scaledLhs, scaledRhs; 4394 R factor = scale[rowIdx] / scale[i]; 4395 4396 if(factor > 0) 4397 { 4398 scaledLhs = (lp.lhs(i) <= R(-infinity)) ? R(-infinity) : lp.lhs(i) * factor; 4399 scaledRhs = (lp.rhs(i) >= R(infinity)) ? R(infinity) : lp.rhs(i) * factor; 4400 } 4401 else 4402 { 4403 scaledLhs = (lp.rhs(i) >= R(infinity)) ? R(-infinity) : lp.rhs(i) * factor; 4404 scaledRhs = (lp.lhs(i) <= R(-infinity)) ? R(infinity) : lp.lhs(i) * factor; 4405 } 4406 4407 if(scaledLhs > maxLhs) 4408 { 4409 maxLhs = scaledLhs; 4410 maxLhsIdx = i; 4411 } 4412 4413 if(scaledRhs < minRhs) 4414 { 4415 minRhs = scaledRhs; 4416 minRhsIdx = i; 4417 } 4418 4419 remRow[i] = true; 4420 } 4421 } 4422 4423 if(rowIdx != -1) 4424 { 4425 R newLhs = (maxLhs > lp.lhs(rowIdx)) ? maxLhs : lp.lhs(rowIdx); 4426 R newRhs = (minRhs < lp.rhs(rowIdx)) ? minRhs : lp.rhs(rowIdx); 4427 4428 if(k == idxLastDupRows) 4429 { 4430 DataArray<int> da_perm(nRowsOld_tmp); 4431 4432 for(int j = 0; j < nRowsOld_tmp; ++j) 4433 { 4434 da_perm[j] = perm_tmp[j]; 4435 } 4436 4437 std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx, 4438 m_dupRows[k], scale, da_perm, isLhsEqualRhs, true, 4439 EQrel(newLhs, newRhs, this->tolerances()->epsilon()), this->_tolerances, k == idxFirstDupRows)); 4440 m_hist.append(ptr); 4441 } 4442 else 4443 { 4444 DataArray<int> da_perm_empty(0); 4445 std::shared_ptr<PostStep> ptr(new DuplicateRowsPS(lp, rowIdx, maxLhsIdx, minRhsIdx, 4446 m_dupRows[k], scale, da_perm_empty, isLhsEqualRhs, false, EQrel(newLhs, newRhs, 4447 this->tolerances()->epsilon()), 4448 this->_tolerances, k == idxFirstDupRows)); 4449 m_hist.append(ptr); 4450 } 4451 4452 if(maxLhs > lp.lhs(rowIdx) || minRhs < lp.rhs(rowIdx)) 4453 { 4454 // modify lhs and rhs of constraint rowIdx 4455 doChangeRanges = true; 4456 4457 if(LTrel(newRhs, newLhs, feastol())) 4458 { 4459 SPxOut::debug(this, "IMAISM55 duplicate rows yield infeasible bounds: lhs={} rhs={}\n", newLhs, 4460 newRhs); 4461 spx_free(perm_tmp); 4462 return this->INFEASIBLE; 4463 } 4464 4465 // if we accept the infeasibility we should clean up the values to avoid problems later 4466 if(newRhs < newLhs) 4467 newRhs = newLhs; 4468 4469 newLhsVec[rowIdx] = newLhs; 4470 newRhsVec[rowIdx] = newRhs; 4471 } 4472 } 4473 } 4474 } 4475 4476 // change ranges for all modified constraints by one single call (more efficient) 4477 if(doChangeRanges) 4478 { 4479 lp.changeRange(newLhsVec, newRhsVec); 4480 } 4481 4482 // remove all rows by one single method call (more efficient) 4483 const int nRowsOld = lp.nRows(); 4484 int* perm = 0; 4485 spx_alloc(perm, nRowsOld); 4486 4487 for(int i = 0; i < nRowsOld; ++i) 4488 { 4489 if(remRow[i]) 4490 { 4491 perm[i] = -1; 4492 ++remRows; 4493 remNzos += lp.rowVector(i).size(); 4494 } 4495 else 4496 perm[i] = 0; 4497 } 4498 4499 lp.removeRows(perm); 4500 4501 for(int i = 0; i < nRowsOld; ++i) 4502 { 4503 // assert that the pre-computed permutation was correct 4504 assert(perm[i] == perm_tmp[i]); 4505 4506 // update the global index mapping 4507 if(perm[i] >= 0) 4508 m_rIdx[perm[i]] = m_rIdx[i]; 4509 } 4510 4511 spx_free(perm); 4512 spx_free(perm_tmp); 4513 4514 if(remRows + remNzos > 0) 4515 { 4516 this->m_remRows += remRows; 4517 this->m_remNzos += remNzos; 4518 4519 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate rows) removed " 4520 << remRows << " rows, " 4521 << remNzos << " non-zeros" 4522 << std::endl;) 4523 4524 if(remRows > this->m_minReduction * oldRows) 4525 again = true; 4526 4527 } 4528 4529 return this->OKAY; 4530 } 4531 4532 template <class R> 4533 typename SPxSimplifier<R>::Result SPxMainSM<R>::duplicateCols(SPxLPBase<R>& lp, bool& again) 4534 { 4535 4536 // This method simplifies the LP by removing duplicate columns 4537 // Duplicates are detected using the algorithm of Bixby and Wagner [1987] 4538 4539 int remCols = 0; 4540 int remNzos = 0; 4541 4542 // remove empty rows and columns 4543 typename SPxSimplifier<R>::Result ret = removeEmpty(lp); 4544 4545 if(ret != this->OKAY) 4546 return ret; 4547 4548 if(lp.nCols() < 2) 4549 return this->OKAY; 4550 4551 DataArray<int> pClass(lp.nCols()); // class of parallel columns 4552 DataArray<int> classSize(lp.nCols()); // size of each class 4553 Array<R> scale(lp.nCols()); // scaling factor for each column 4554 int* idxMem = 0; 4555 4556 try 4557 { 4558 spx_alloc(idxMem, lp.nCols()); 4559 } 4560 catch(const SPxMemoryException& x) 4561 { 4562 spx_free(idxMem); 4563 throw x; 4564 } 4565 4566 IdxSet idxSet(lp.nCols(), idxMem); // set of feasible indices for new pClass 4567 4568 // init 4569 pClass[0] = 0; 4570 scale[0] = 0.0; 4571 classSize[0] = lp.nCols(); 4572 4573 for(int j = 1; j < lp.nCols(); ++j) 4574 { 4575 pClass[j] = 0; 4576 scale[j] = 0.0; 4577 classSize[j] = 0; 4578 idxSet.addIdx(j); 4579 } 4580 4581 R oldVal = 0.0; 4582 4583 // main loop 4584 for(int i = 0; i < lp.nRows(); ++i) 4585 { 4586 const SVectorBase<R>& row = lp.rowVector(i); 4587 4588 for(int k = 0; k < row.size(); ++k) 4589 { 4590 R aij = row.value(k); 4591 int j = row.index(k); 4592 4593 if(scale[j] == 0.0) 4594 scale[j] = aij; 4595 4596 m_classSetCols[pClass[j]].add(j, aij / scale[j]); 4597 4598 if(--classSize[pClass[j]] == 0) 4599 idxSet.addIdx(pClass[j]); 4600 } 4601 4602 // update each parallel class with non-zero row entry 4603 for(int m = 0; m < row.size(); ++m) 4604 { 4605 int k = pClass[row.index(m)]; 4606 4607 if(m_classSetCols[k].size() > 0) 4608 { 4609 // sort classSet[k] w.r.t. scaled row values 4610 ElementCompare compare(this->tolerances()->epsilon()); 4611 4612 if(m_classSetCols[k].size() > 1) 4613 SPxQuicksort(m_classSetCols[k].mem(), m_classSetCols[k].size(), compare); 4614 4615 // use new index first 4616 int classIdx = idxSet.index(0); 4617 idxSet.remove(0); 4618 4619 for(int l = 0; l < m_classSetCols[k].size(); ++l) 4620 { 4621 if(l != 0 && NErel(m_classSetCols[k].value(l), oldVal, this->epsZero())) 4622 { 4623 // start new parallel class 4624 classIdx = idxSet.index(0); 4625 idxSet.remove(0); 4626 } 4627 4628 pClass[m_classSetCols[k].index(l)] = classIdx; 4629 ++classSize[classIdx]; 4630 4631 oldVal = m_classSetCols[k].value(l); 4632 } 4633 4634 m_classSetCols[k].clear(); 4635 } 4636 } 4637 } 4638 4639 spx_free(idxMem); 4640 4641 DataArray<bool> remCol(lp.nCols()); 4642 DataArray<bool> fixAndRemCol(lp.nCols()); 4643 4644 for(int k = 0; k < lp.nCols(); ++k) 4645 m_dupCols[k].clear(); 4646 4647 for(int k = 0; k < lp.nCols(); ++k) 4648 { 4649 remCol[k] = false; 4650 fixAndRemCol[k] = false; 4651 m_dupCols[pClass[k]].add(k, 0.0); 4652 } 4653 4654 bool hasDuplicateCol = false; 4655 DataArray<int> m_perm_empty(0); 4656 4657 for(int k = 0; k < lp.nCols(); ++k) 4658 { 4659 if(m_dupCols[k].size() > 1 && !(lp.colVector(m_dupCols[k].index(0)).size() == 1)) 4660 { 4661 SPxOut::debug(this, "IMAISM58 {} duplicate columns found\n", m_dupCols[k].size()); 4662 4663 if(!hasDuplicateCol) 4664 { 4665 std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, m_perm_empty, this->_tolerances, 4666 true)); 4667 m_hist.append(ptr); 4668 hasDuplicateCol = true; 4669 } 4670 4671 for(int l = 0; l < m_dupCols[k].size(); ++l) 4672 { 4673 for(int m = 0; m < m_dupCols[k].size(); ++m) 4674 { 4675 int j1 = m_dupCols[k].index(l); 4676 int j2 = m_dupCols[k].index(m); 4677 4678 if(l != m && !remCol[j1] && !remCol[j2]) 4679 { 4680 R cj1 = lp.maxObj(j1); 4681 R cj2 = lp.maxObj(j2); 4682 4683 // A.j1 = factor * A.j2 4684 R factor = scale[j1] / scale[j2]; 4685 R objDif = cj1 - cj2 * scale[j1] / scale[j2]; 4686 4687 SOPLEX_ASSERT_WARN("WMAISM59", isNotZero(factor, this->epsZero())); 4688 4689 if(isZero(objDif, this->epsZero())) 4690 { 4691 // case 1: objectives also duplicate 4692 4693 // if 0 is not within the column bounds, we are not able to postsolve if the aggregated column has 4694 // status ZERO, hence we skip this case 4695 if(LErel(lp.lower(j1), R(0.0), this->tolerances()->epsilon()) 4696 && GErel(lp.upper(j1), R(0.0), this->tolerances()->epsilon()) 4697 && LErel(lp.lower(j2), R(0.0), this->tolerances()->epsilon()) 4698 && GErel(lp.upper(j2), R(0.0), this->tolerances()->epsilon())) 4699 { 4700 std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, j1, j2, factor, m_perm_empty, 4701 this->_tolerances)); 4702 // variable substitution xj2' := xj2 + factor * xj1 <=> xj2 = -factor * xj1 + xj2' 4703 m_hist.append(ptr); 4704 4705 // update bounds of remaining column j2 (new column j2') 4706 if(factor > 0) 4707 { 4708 if(lp.lower(j2) <= R(-infinity) || lp.lower(j1) <= R(-infinity)) 4709 lp.changeLower(j2, R(-infinity)); 4710 else 4711 lp.changeLower(j2, lp.lower(j2) + factor * lp.lower(j1)); 4712 4713 if(lp.upper(j2) >= R(infinity) || lp.upper(j1) >= R(infinity)) 4714 lp.changeUpper(j2, R(infinity)); 4715 else 4716 lp.changeUpper(j2, lp.upper(j2) + factor * lp.upper(j1)); 4717 } 4718 else if(factor < 0) 4719 { 4720 if(lp.lower(j2) <= R(-infinity) || lp.upper(j1) >= R(infinity)) 4721 lp.changeLower(j2, R(-infinity)); 4722 else 4723 lp.changeLower(j2, lp.lower(j2) + factor * lp.upper(j1)); 4724 4725 if(lp.upper(j2) >= R(infinity) || lp.lower(j1) <= R(-infinity)) 4726 lp.changeUpper(j2, R(infinity)); 4727 else 4728 lp.changeUpper(j2, lp.upper(j2) + factor * lp.lower(j1)); 4729 } 4730 4731 SPxOut::debug(this, "IMAISM60 two duplicate columns {}, {} replaced by one\n", j1, j2); 4732 4733 remCol[j1] = true; 4734 4735 ++m_stat[SUB_DUPLICATE_COL]; 4736 } 4737 else 4738 { 4739 SPxOut::debug(this, 4740 "IMAISM80 not removing two duplicate columns {}, {} because zero not contained in their bounds\n", 4741 j1, j2); 4742 } 4743 } 4744 else 4745 { 4746 // case 2: objectives not duplicate 4747 // considered for maximization sense 4748 if(lp.lower(j2) <= R(-infinity)) 4749 { 4750 if(factor > 0 && objDif > 0) 4751 { 4752 if(lp.upper(j1) >= R(infinity)) 4753 { 4754 SPxOut::debug(this, "IMAISM75 LP unbounded\n"); 4755 return this->UNBOUNDED; 4756 } 4757 4758 // fix j1 at upper bound 4759 SPxOut::debug(this, "IMAISM61 two duplicate columns {}, {} first one fixed at upper bound={}\n", j1, 4760 j2, lp.upper(j1)); 4761 4762 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1), this->_tolerances)); 4763 m_hist.append(ptr); 4764 lp.changeLower(j1, lp.upper(j1)); 4765 } 4766 else if(factor < 0 && objDif < 0) 4767 { 4768 if(lp.lower(j1) <= R(-infinity)) 4769 { 4770 SPxOut::debug(this, "IMAISM76 LP unbounded\n"); 4771 return this->UNBOUNDED; 4772 } 4773 4774 // fix j1 at lower bound 4775 SPxOut::debug(this, "IMAISM62 two duplicate columns {}, {} first one fixed at lower bound={}\n", j1, 4776 j2, lp.lower(j1)); 4777 4778 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1), this->_tolerances)); 4779 m_hist.append(ptr); 4780 lp.changeUpper(j1, lp.lower(j1)); 4781 } 4782 } 4783 else if(lp.upper(j2) >= R(infinity)) 4784 { 4785 // fix j1 at upper bound 4786 if(factor < 0 && objDif > 0) 4787 { 4788 if(lp.upper(j1) >= R(infinity)) 4789 { 4790 SPxOut::debug(this, "IMAISM77 LP unbounded\n"); 4791 return this->UNBOUNDED; 4792 } 4793 4794 // fix j1 at upper bound 4795 SPxOut::debug(this, "IMAISM63 two duplicate columns {}, {}, first one fixed at upper bound={}\n", 4796 j1, j2, lp.upper(j1)); 4797 4798 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.upper(j1), this->_tolerances)); 4799 m_hist.append(ptr); 4800 lp.changeLower(j1, lp.upper(j1)); 4801 } 4802 4803 // fix j1 at lower bound 4804 else if(factor > 0 && objDif < 0) 4805 { 4806 if(lp.lower(j1) <= R(-infinity)) 4807 { 4808 SPxOut::debug(this, "IMAISM78 LP unbounded\n"); 4809 return this->UNBOUNDED; 4810 } 4811 4812 // fix j1 at lower bound 4813 SPxOut::debug(this, "IMAISM64 two duplicate columns {}, {} first one fixed at lower bound={}\n", j1, 4814 j2, lp.lower(j1)); 4815 4816 std::shared_ptr<PostStep> ptr(new FixBoundsPS(lp, j1, lp.lower(j1), this->_tolerances)); 4817 m_hist.append(ptr); 4818 lp.changeUpper(j1, lp.lower(j1)); 4819 } 4820 } 4821 4822 if(EQrel(lp.lower(j1), lp.upper(j1), feastol())) 4823 { 4824 remCol[j1] = true; 4825 fixAndRemCol[j1] = true; 4826 4827 ++m_stat[FIX_DUPLICATE_COL]; 4828 } 4829 } 4830 } 4831 } 4832 } 4833 } 4834 } 4835 4836 for(int j = 0; j < lp.nCols(); ++j) 4837 { 4838 if(fixAndRemCol[j]) 4839 { 4840 assert(remCol[j]); 4841 4842 // correctIdx == false, because the index mapping will be handled by the postsolving in DuplicateColsPS 4843 fixColumn(lp, j, false); 4844 } 4845 } 4846 4847 // remove all columns by one single method call (more efficient) 4848 const int nColsOld = lp.nCols(); 4849 int* perm = 0; 4850 spx_alloc(perm, nColsOld); 4851 4852 for(int j = 0; j < nColsOld; ++j) 4853 { 4854 if(remCol[j]) 4855 { 4856 perm[j] = -1; 4857 ++remCols; 4858 remNzos += lp.colVector(j).size(); 4859 } 4860 else 4861 perm[j] = 0; 4862 } 4863 4864 lp.removeCols(perm); 4865 4866 for(int j = 0; j < nColsOld; ++j) 4867 { 4868 if(perm[j] >= 0) 4869 m_cIdx[perm[j]] = m_cIdx[j]; 4870 } 4871 4872 DataArray<int> da_perm(nColsOld); 4873 4874 for(int j = 0; j < nColsOld; ++j) 4875 { 4876 da_perm[j] = perm[j]; 4877 } 4878 4879 if(hasDuplicateCol) 4880 { 4881 std::shared_ptr<PostStep> ptr(new DuplicateColsPS(lp, 0, 0, 1.0, da_perm, this->_tolerances, false, 4882 true)); 4883 m_hist.append(ptr); 4884 } 4885 4886 spx_free(perm); 4887 4888 assert(remCols > 0 || remNzos == 0); 4889 4890 if(remCols > 0) 4891 { 4892 this->m_remCols += remCols; 4893 this->m_remNzos += remNzos; 4894 4895 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier (duplicate columns) removed " 4896 << remCols << " cols, " 4897 << remNzos << " non-zeros" 4898 << std::endl;) 4899 4900 if(remCols > this->m_minReduction * nColsOld) 4901 again = true; 4902 } 4903 4904 return this->OKAY; 4905 } 4906 4907 template <class R> 4908 void SPxMainSM<R>::fixColumn(SPxLPBase<R>& lp, int j, bool correctIdx) 4909 { 4910 4911 assert(EQrel(lp.lower(j), lp.upper(j), feastol())); 4912 4913 R lo = lp.lower(j); 4914 R up = lp.upper(j); 4915 const SVectorBase<R>& col = lp.colVector(j); 4916 R mid = lo; 4917 4918 // use the center value between slightly different bounds to improve numerics 4919 if(NE(lo, up, this->tolerances()->epsilon())) 4920 mid = (up + lo) / 2.0; 4921 4922 assert(lo < R(infinity) && lo > R(-infinity)); 4923 assert(up < R(infinity) && up > R(-infinity)); 4924 4925 SPxOut::debug(this, "IMAISM66 fix variable x{}: lower={} upper={} to new value: {}\n", j, lo, up, 4926 mid); 4927 4928 if(isNotZero(lo, this->epsZero())) 4929 { 4930 for(int k = 0; k < col.size(); ++k) 4931 { 4932 int i = col.index(k); 4933 4934 if(lp.rhs(i) < R(infinity)) 4935 { 4936 R y = mid * col.value(k); 4937 R scale = maxAbs(lp.rhs(i), y); 4938 4939 if(scale < 1.0) 4940 scale = 1.0; 4941 4942 R rhs = (lp.rhs(i) / scale) - (y / scale); 4943 4944 if(isZero(rhs, this->epsZero())) 4945 rhs = 0.0; 4946 else 4947 rhs *= scale; 4948 4949 SPxOut::debug(this, "IMAISM67 row {}: rhs={} ({}) aij={}\n", i, rhs, lp.rhs(i), col.value(k)); 4950 4951 lp.changeRhs(i, rhs); 4952 } 4953 4954 if(lp.lhs(i) > R(-infinity)) 4955 { 4956 R y = mid * col.value(k); 4957 R scale = maxAbs(lp.lhs(i), y); 4958 4959 if(scale < 1.0) 4960 scale = 1.0; 4961 4962 R lhs = (lp.lhs(i) / scale) - (y / scale); 4963 4964 if(isZero(lhs, this->epsZero())) 4965 lhs = 0.0; 4966 else 4967 lhs *= scale; 4968 4969 SPxOut::debug(this, "IMAISM68 row {}: lhs={} ({}) aij={}\n", i, lhs, lp.lhs(i), col.value(k)); 4970 4971 lp.changeLhs(i, lhs); 4972 } 4973 4974 assert(lp.lhs(i) <= lp.rhs(i) + feastol()); 4975 } 4976 } 4977 4978 std::shared_ptr<PostStep> ptr(new FixVariablePS(lp, *this, j, lp.lower(j), this->_tolerances, 4979 correctIdx)); 4980 m_hist.append(ptr); 4981 } 4982 4983 template <class R> 4984 typename SPxSimplifier<R>::Result SPxMainSM<R>::simplify(SPxLPBase<R>& lp, Real remainingTime, 4985 bool keepbounds, uint32_t seed) 4986 { 4987 // transfer message handler 4988 this->spxout = lp.spxout; 4989 assert(this->spxout != 0); 4990 4991 m_thesense = lp.spxSense(); 4992 this->m_timeUsed->reset(); 4993 this->m_timeUsed->start(); 4994 4995 this->m_objoffset = 0.0; 4996 m_cutoffbound = R(-infinity); 4997 m_pseudoobj = R(-infinity); 4998 4999 this->m_remRows = 0; 5000 this->m_remCols = 0; 5001 this->m_remNzos = 0; 5002 this->m_chgBnds = 0; 5003 this->m_chgLRhs = 0; 5004 this->m_keptBnds = 0; 5005 this->m_keptLRhs = 0; 5006 5007 m_result = this->OKAY; 5008 bool again = true; 5009 int nrounds = 0; 5010 5011 if(m_hist.size() > 0) 5012 { 5013 m_hist.clear(); 5014 } 5015 5016 m_hist.reSize(0); 5017 m_postsolved = false; 5018 5019 SPX_MSG_INFO2((*this->spxout), 5020 int numRangedRows = 0; 5021 int numBoxedCols = 0; 5022 int numEqualities = 0; 5023 5024 for(int i = 0; i < lp.nRows(); ++i) 5025 { 5026 if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity)) 5027 { 5028 if(EQ(lp.lhs(i), lp.rhs(i), this->tolerances()->epsilon())) 5029 ++numEqualities; 5030 else 5031 ++numRangedRows; 5032 } 5033 } 5034 for(int j = 0; j < lp.nCols(); ++j) 5035 if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity)) 5036 ++numBoxedCols; 5037 5038 (*this->spxout) << "LP has " 5039 << numEqualities << " equations, " 5040 << numRangedRows << " ranged rows, " 5041 << numBoxedCols << " boxed columns" 5042 << std::endl; 5043 ) 5044 5045 m_stat.reSize(17); 5046 5047 for(int k = 0; k < m_stat.size(); ++k) 5048 m_stat[k] = 0; 5049 5050 m_addedcols = 0; 5051 handleRowObjectives(lp); 5052 5053 m_prim.reDim(lp.nCols()); 5054 m_slack.reDim(lp.nRows()); 5055 m_dual.reDim(lp.nRows()); 5056 m_redCost.reDim(lp.nCols()); 5057 m_cBasisStat.reSize(lp.nCols()); 5058 m_rBasisStat.reSize(lp.nRows()); 5059 m_cIdx.reSize(lp.nCols()); 5060 m_rIdx.reSize(lp.nRows()); 5061 5062 m_classSetRows.reSize(lp.nRows()); 5063 m_classSetCols.reSize(lp.nCols()); 5064 m_dupRows.reSize(lp.nRows()); 5065 m_dupCols.reSize(lp.nCols()); 5066 5067 m_keepbounds = keepbounds; 5068 5069 for(int i = 0; i < lp.nRows(); ++i) 5070 m_rIdx[i] = i; 5071 5072 for(int j = 0; j < lp.nCols(); ++j) 5073 m_cIdx[j] = j; 5074 5075 // round extreme values (set all values smaller than this->eps to zero and all values bigger than R(infinity)/5 to R(infinity)) 5076 #if SOPLEX_EXTREMES 5077 handleExtremes(lp); 5078 #endif 5079 5080 // main presolving loop 5081 while(again && m_result == this->OKAY) 5082 { 5083 nrounds++; 5084 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "Round " << nrounds << ":" << std::endl;) 5085 again = false; 5086 5087 #if SOPLEX_ROWS_SPXMAINSM 5088 5089 if(m_result == this->OKAY) 5090 m_result = simplifyRows(lp, again); 5091 5092 #endif 5093 5094 #if SOPLEX_COLS_SPXMAINSM 5095 5096 if(m_result == this->OKAY) 5097 m_result = simplifyCols(lp, again); 5098 5099 #endif 5100 5101 if(again) 5102 continue; 5103 5104 #if SOPLEX_DUAL_SPXMAINSM 5105 5106 if(m_result == this->OKAY) 5107 m_result = simplifyDual(lp, again); 5108 5109 #endif 5110 5111 if(again) 5112 continue; 5113 5114 #if SOPLEX_DUPLICATE_ROWS 5115 5116 if(m_result == this->OKAY) 5117 m_result = duplicateRows(lp, again); 5118 5119 #endif 5120 5121 #if SOPLEX_DUPLICATE_COLS 5122 5123 if(m_result == this->OKAY) 5124 m_result = duplicateCols(lp, again); 5125 5126 #endif 5127 5128 if(!again) 5129 { 5130 #if SOPLEX_TRIVIAL_HEURISTICS 5131 trivialHeuristic(lp); 5132 #endif 5133 5134 #if SOPLEX_PSEUDOOBJ 5135 propagatePseudoobj(lp); 5136 #endif 5137 5138 #if SOPLEX_MULTI_AGGREGATE 5139 5140 if(m_result == this->OKAY) 5141 m_result = multiaggregation(lp, again); 5142 5143 #endif 5144 } 5145 5146 } 5147 5148 // preprocessing detected infeasibility or unboundedness 5149 if(m_result != this->OKAY) 5150 { 5151 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier result: " << static_cast<int> 5152 (m_result) << std::endl;) 5153 return m_result; 5154 } 5155 5156 this->m_remCols -= m_addedcols; 5157 this->m_remNzos -= m_addedcols; 5158 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed " 5159 << this->m_remRows << " rows, " 5160 << this->m_remCols << " columns, " 5161 << this->m_remNzos << " nonzeros, " 5162 << this->m_chgBnds << " col bounds, " 5163 << this->m_chgLRhs << " row bounds" 5164 << std::endl;) 5165 5166 if(keepbounds) 5167 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "Simplifier kept " 5168 << this->m_keptBnds << " column bounds, " 5169 << this->m_keptLRhs << " row bounds" 5170 << std::endl;) 5171 5172 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Reduced LP has " 5173 << lp.nRows() << " rows " 5174 << lp.nCols() << " columns " 5175 << lp.nNzos() << " nonzeros" 5176 << std::endl;) 5177 5178 SPX_MSG_INFO2((*this->spxout), 5179 int numRangedRows = 0; 5180 int numBoxedCols = 0; 5181 int numEqualities = 0; 5182 5183 for(int i = 0; i < lp.nRows(); ++i) 5184 { 5185 if(lp.lhs(i) > R(-infinity) && lp.rhs(i) < R(infinity)) 5186 { 5187 if(EQ(lp.lhs(i), lp.rhs(i), this->tolerances()->epsilon())) 5188 ++numEqualities; 5189 else 5190 ++numRangedRows; 5191 } 5192 } 5193 for(int j = 0; j < lp.nCols(); ++j) 5194 if(lp.lower(j) > R(-infinity) && lp.upper(j) < R(infinity)) 5195 ++numBoxedCols; 5196 5197 (*this->spxout) << "Reduced LP has " 5198 << numEqualities << " equations, " 5199 << numRangedRows << " ranged rows, " 5200 << numBoxedCols << " boxed columns" 5201 << std::endl; 5202 ) 5203 5204 if(lp.nCols() == 0 && lp.nRows() == 0) 5205 { 5206 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Simplifier removed all rows and columns" << 5207 std::endl;) 5208 m_result = this->VANISHED; 5209 } 5210 5211 SPX_MSG_INFO2((*this->spxout), (*this->spxout) << "\nSimplifier performed:\n" 5212 << m_stat[EMPTY_ROW] << " empty rows\n" 5213 << m_stat[FREE_ROW] << " free rows\n" 5214 << m_stat[SINGLETON_ROW] << " singleton rows\n" 5215 << m_stat[FORCE_ROW] << " forcing rows\n" 5216 << m_stat[EMPTY_COL] << " empty columns\n" 5217 << m_stat[FIX_COL] << " fixed columns\n" 5218 << m_stat[FREE_ZOBJ_COL] << " free columns with zero objective\n" 5219 << m_stat[ZOBJ_SINGLETON_COL] << " singleton columns with zero objective\n" 5220 << m_stat[DOUBLETON_ROW] << " singleton columns combined with a doubleton equation\n" 5221 << m_stat[FREE_SINGLETON_COL] << " free singleton columns\n" 5222 << m_stat[DOMINATED_COL] << " dominated columns\n" 5223 << m_stat[WEAKLY_DOMINATED_COL] << " weakly dominated columns\n" 5224 << m_stat[DUPLICATE_ROW] << " duplicate rows\n" 5225 << m_stat[FIX_DUPLICATE_COL] << " duplicate columns (fixed)\n" 5226 << m_stat[SUB_DUPLICATE_COL] << " duplicate columns (substituted)\n" 5227 << m_stat[AGGREGATION] << " variable aggregations\n" 5228 << m_stat[MULTI_AGG] << " multi aggregations\n" 5229 << std::endl;); 5230 5231 this->m_timeUsed->stop(); 5232 5233 return m_result; 5234 } 5235 5236 template <class R> 5237 void SPxMainSM<R>::unsimplify(const VectorBase<R>& x, const VectorBase<R>& y, 5238 const VectorBase<R>& s, const VectorBase<R>& r, 5239 const typename SPxSolverBase<R>::VarStatus rows[], 5240 const typename SPxSolverBase<R>::VarStatus cols[], bool isOptimal) 5241 { 5242 SPX_MSG_INFO1((*this->spxout), 5243 (*this->spxout) << " --- unsimplifying solution and basis" << std::endl;) 5244 assert(x.dim() <= m_prim.dim()); 5245 assert(y.dim() <= m_dual.dim()); 5246 assert(x.dim() == r.dim()); 5247 assert(y.dim() == s.dim()); 5248 5249 // assign values of variables in reduced LP 5250 // NOTE: for maximization problems, we have to switch signs of dual and reduced cost values, 5251 // since simplifier assumes minimization problem 5252 for(int j = 0; j < x.dim(); ++j) 5253 { 5254 m_prim[j] = isZero(x[j], this->epsZero()) ? 0.0 : x[j]; 5255 m_redCost[j] = isZero(r[j], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -r[j] : 5256 r[j]); 5257 m_cBasisStat[j] = cols[j]; 5258 } 5259 5260 for(int i = 0; i < y.dim(); ++i) 5261 { 5262 m_dual[i] = isZero(y[i], this->epsZero()) ? 0.0 : (m_thesense == SPxLPBase<R>::MAXIMIZE ? -y[i] : 5263 y[i]); 5264 m_slack[i] = isZero(s[i], this->epsZero()) ? 0.0 : s[i]; 5265 m_rBasisStat[i] = rows[i]; 5266 } 5267 5268 // undo preprocessing 5269 for(int k = m_hist.size() - 1; k >= 0; --k) 5270 { 5271 SPxOut::debug(this, "unsimplifying {} \n", m_hist[k]->getName()); 5272 5273 try 5274 { 5275 m_hist[k]->execute(m_prim, m_dual, m_slack, m_redCost, m_cBasisStat, m_rBasisStat, isOptimal); 5276 } 5277 catch(const SPxException& ex) 5278 { 5279 SPX_MSG_INFO1((*this->spxout), (*this->spxout) << "Exception thrown while unsimplifying " << 5280 m_hist[k]->getName() << ":\n" << ex.what() << "\n"); 5281 throw SPxInternalCodeException("XMAISM00 Exception thrown during unsimply()."); 5282 } 5283 5284 m_hist.reSize(k); 5285 } 5286 5287 // for maximization problems, we have to switch signs of dual and reduced cost values back 5288 if(m_thesense == SPxLPBase<R>::MAXIMIZE) 5289 { 5290 for(int j = 0; j < m_redCost.dim(); ++j) 5291 m_redCost[j] = -m_redCost[j]; 5292 5293 for(int i = 0; i < m_dual.dim(); ++i) 5294 m_dual[i] = -m_dual[i]; 5295 } 5296 5297 if(m_addedcols > 0) 5298 { 5299 assert(m_prim.dim() >= m_addedcols); 5300 m_prim.reDim(m_prim.dim() - m_addedcols); 5301 m_redCost.reDim(m_redCost.dim() - m_addedcols); 5302 m_cBasisStat.reSize(m_cBasisStat.size() - m_addedcols); 5303 m_cIdx.reSize(m_cIdx.size() - m_addedcols); 5304 } 5305 5306 #ifdef SOPLEX_CHECK_BASIS_DIM 5307 int numBasis = 0; 5308 5309 for(int rs = 0; rs < m_rBasisStat.size(); ++rs) 5310 { 5311 if(m_rBasisStat[rs] == SPxSolverBase<R>::BASIC) 5312 numBasis ++; 5313 } 5314 5315 for(int cs = 0; cs < m_cBasisStat.size(); ++cs) 5316 { 5317 if(m_cBasisStat[cs] == SPxSolverBase<R>::BASIC) 5318 numBasis ++; 5319 } 5320 5321 if(numBasis != m_rBasisStat.size()) 5322 { 5323 throw SPxInternalCodeException("XMAISM26 Dimension doesn't match after this step."); 5324 } 5325 5326 #endif 5327 5328 m_hist.clear(); 5329 m_postsolved = true; 5330 } 5331 5332 // Pretty-printing of solver status. 5333 template <class R> 5334 std::ostream& operator<<(std::ostream& os, const typename SPxSimplifier<R>::Result& status) 5335 { 5336 switch(status) 5337 { 5338 case SPxSimplifier<R>::OKAY: 5339 os << "SUCCESS"; 5340 break; 5341 5342 case SPxSimplifier<R>::INFEASIBLE: 5343 os << "INFEASIBLE"; 5344 break; 5345 5346 case SPxSimplifier<R>::DUAL_INFEASIBLE: 5347 os << "DUAL_INFEASIBLE"; 5348 break; 5349 5350 case SPxSimplifier<R>::UNBOUNDED: 5351 os << "UNBOUNDED"; 5352 break; 5353 5354 case SPxSimplifier<R>::VANISHED: 5355 os << "VANISHED"; 5356 break; 5357 5358 default: 5359 os << "UNKNOWN"; 5360 break; 5361 } 5362 5363 return os; 5364 } 5365 5366 } //namespace soplex 5367