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/spxdefines.h" 28 #include "soplex/spxbasis.h" 29 #include "soplex/spxsolver.h" 30 #include "soplex/spxout.h" 31 32 namespace soplex 33 { 34 template <class R> 35 void SPxBasisBase<R>::reDim() 36 { 37 38 assert(theLP != 0); 39 40 SPxOut::debug(this, "DCHBAS01 SPxBasisBase<R>::reDim(): matrixIsSetup={}, factorized={}\n", 41 matrixIsSetup, factorized); 42 43 thedesc.reSize(theLP->nRows(), theLP->nCols()); 44 45 if(theLP->dim() != matrix.size()) 46 { 47 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 48 "ICHBAS02 basis redimensioning invalidates factorization" 49 << std::endl;) 50 51 matrix.reSize(theLP->dim()); 52 theBaseId.reSize(theLP->dim()); 53 matrixIsSetup = false; 54 factorized = false; 55 } 56 57 SPxOut::debug(this, "DCHBAS03 SPxBasisBase<R>::reDim(): --> matrixIsSetup={}, factorized={}\n", 58 matrixIsSetup, factorized); 59 60 assert(matrix.size() >= theLP->dim()); 61 assert(theBaseId.size() >= theLP->dim()); 62 } 63 64 /* adapt basis and basis descriptor to added rows */ 65 template <class R> 66 void SPxBasisBase<R>::addedRows(int n) 67 { 68 assert(theLP != 0); 69 70 if(n > 0) 71 { 72 reDim(); 73 74 if(theLP->rep() == SPxSolverBase<R>::COLUMN) 75 { 76 /* after adding rows in column representation, reDim() should set these bools to false. */ 77 assert(!matrixIsSetup && !factorized); 78 79 for(int i = theLP->nRows() - n; i < theLP->nRows(); ++i) 80 { 81 thedesc.rowStatus(i) = dualRowStatus(i); 82 baseId(i) = theLP->SPxLPBase<R>::rId(i); 83 } 84 } 85 else 86 { 87 assert(theLP->rep() == SPxSolverBase<R>::ROW); 88 89 for(int i = theLP->nRows() - n; i < theLP->nRows(); ++i) 90 thedesc.rowStatus(i) = dualRowStatus(i); 91 } 92 93 /* If matrix was set up, load new basis vectors to the matrix. 94 * In the row case, the basis is not effected by adding rows. However, 95 * since @c matrix stores references to the rows in the LP (SPxLPBase<R>), a realloc 96 * in SPxLPBase<R> (e.g. due to space requirements) might invalidate these references. 97 * We therefore have to "reload" the matrix if it is set up. Note that reDim() 98 * leaves @c matrixIsSetup untouched if only row have been added, since the basis 99 * matrix already has the correct size. */ 100 if(status() > NO_PROBLEM && matrixIsSetup) 101 loadMatrixVecs(); 102 103 /* update basis status */ 104 switch(status()) 105 { 106 case PRIMAL: 107 case UNBOUNDED: 108 setStatus(REGULAR); 109 break; 110 111 case OPTIMAL: 112 case INFEASIBLE: 113 setStatus(DUAL); 114 break; 115 116 case NO_PROBLEM: 117 case SINGULAR: 118 case REGULAR: 119 case DUAL: 120 break; 121 122 default: 123 SPX_MSG_ERROR(std::cerr << "ECHBAS04 Unknown basis status!" << std::endl;) 124 throw SPxInternalCodeException("XCHBAS01 This should never happen."); 125 } 126 } 127 } 128 129 template <class R> 130 void SPxBasisBase<R>::removedRow(int i) 131 { 132 133 assert(status() > NO_PROBLEM); 134 assert(theLP != 0); 135 136 if(theLP->rep() == SPxSolverBase<R>::ROW) 137 { 138 if(theLP->isBasic(thedesc.rowStatus(i))) 139 { 140 setStatus(NO_PROBLEM); 141 factorized = false; 142 143 SPxOut::debug(this, "DCHBAS05 Warning: deleting basic row!\n"); 144 } 145 } 146 else 147 { 148 assert(theLP->rep() == SPxSolverBase<R>::COLUMN); 149 factorized = false; 150 151 if(!theLP->isBasic(thedesc.rowStatus(i))) 152 { 153 setStatus(NO_PROBLEM); 154 SPxOut::debug(this, "DCHBAS06 Warning: deleting nonbasic row!\n"); 155 } 156 else if(status() > NO_PROBLEM && matrixIsSetup) 157 { 158 for(int j = theLP->dim(); j >= 0; --j) 159 { 160 SPxId id = baseId(j); 161 162 if(id.isSPxRowId() && !theLP->has(SPxRowId(id))) 163 { 164 baseId(j) = baseId(theLP->dim()); 165 166 if(j < theLP->dim()) 167 matrix[j] = &theLP->vector(baseId(j)); 168 169 break; 170 } 171 } 172 } 173 } 174 175 thedesc.rowStatus(i) = thedesc.rowStatus(theLP->nRows()); 176 reDim(); 177 } 178 179 template <class R> 180 void SPxBasisBase<R>::removedRows(const int perm[]) 181 { 182 assert(status() > NO_PROBLEM); 183 assert(theLP != 0); 184 185 int i; 186 int n = thedesc.nRows(); 187 188 if(theLP->rep() == SPxSolverBase<R>::ROW) 189 { 190 for(i = 0; i < n; ++i) 191 { 192 if(perm[i] != i) 193 { 194 if(perm[i] < 0) // row got removed 195 { 196 if(theLP->isBasic(thedesc.rowStatus(i))) 197 { 198 setStatus(NO_PROBLEM); 199 factorized = matrixIsSetup = false; 200 SPxOut::debug(this, "DCHBAS07 Warning: deleting basic row!\n"); 201 } 202 } 203 else // row was moved 204 thedesc.rowStatus(perm[i]) = thedesc.rowStatus(i); 205 } 206 } 207 } 208 else 209 { 210 assert(theLP->rep() == SPxSolverBase<R>::COLUMN); 211 212 factorized = false; 213 matrixIsSetup = false; 214 215 for(i = 0; i < n; ++i) 216 { 217 if(perm[i] != i) 218 { 219 if(perm[i] < 0) // row got removed 220 { 221 if(!theLP->isBasic(thedesc.rowStatus(i))) 222 setStatus(NO_PROBLEM); 223 } 224 else // row was moved 225 thedesc.rowStatus(perm[i]) = thedesc.rowStatus(i); 226 } 227 } 228 } 229 230 reDim(); 231 } 232 233 template <class R> 234 static typename SPxBasisBase<R>::Desc::Status 235 primalColStatus(int i, const SPxLPBase<R>* theLP) 236 { 237 assert(theLP != 0); 238 239 if(theLP->upper(i) < R(infinity)) 240 { 241 if(theLP->lower(i) > R(-infinity)) 242 { 243 if(theLP->lower(i) == theLP->SPxLPBase<R>::upper(i)) 244 return SPxBasisBase<R>::Desc::P_FIXED; 245 /* 246 else 247 return (-theLP->lower(i) < theLP->upper(i)) 248 ? SPxBasisBase<R>::Desc::P_ON_LOWER 249 : SPxBasisBase<R>::Desc::P_ON_UPPER; 250 */ 251 else if(theLP->maxObj(i) == 0) 252 return (-theLP->lower(i) < theLP->upper(i)) 253 ? SPxBasisBase<R>::Desc::P_ON_LOWER 254 : SPxBasisBase<R>::Desc::P_ON_UPPER; 255 else 256 return (theLP->maxObj(i) < 0) 257 ? SPxBasisBase<R>::Desc::P_ON_LOWER 258 : SPxBasisBase<R>::Desc::P_ON_UPPER; 259 } 260 else 261 return SPxBasisBase<R>::Desc::P_ON_UPPER; 262 } 263 else if(theLP->lower(i) > R(-infinity)) 264 return SPxBasisBase<R>::Desc::P_ON_LOWER; 265 else 266 return SPxBasisBase<R>::Desc::P_FREE; 267 } 268 269 270 /* adapt basis and basis descriptor to added columns */ 271 template <class R> 272 void SPxBasisBase<R>::addedCols(int n) 273 { 274 assert(theLP != 0); 275 276 if(n > 0) 277 { 278 reDim(); 279 280 if(theLP->rep() == SPxSolverBase<R>::ROW) 281 { 282 /* after adding columns in row representation, reDim() should set these bools to false. */ 283 assert(!matrixIsSetup && !factorized); 284 285 for(int i = theLP->nCols() - n; i < theLP->nCols(); ++i) 286 { 287 thedesc.colStatus(i) = primalColStatus(i, theLP); 288 baseId(i) = theLP->SPxLPBase<R>::cId(i); 289 } 290 } 291 else 292 { 293 assert(theLP->rep() == SPxSolverBase<R>::COLUMN); 294 295 for(int i = theLP->nCols() - n; i < theLP->nCols(); ++i) 296 thedesc.colStatus(i) = primalColStatus(i, theLP); 297 } 298 299 /* If matrix was set up, load new basis vectors to the matrix 300 * In the column case, the basis is not effected by adding columns. However, 301 * since @c matrix stores references to the columns in the LP (SPxLPBase<R>), a realloc 302 * in SPxLPBase<R> (e.g. due to space requirements) might invalidate these references. 303 * We therefore have to "reload" the matrix if it is set up. Note that reDim() 304 * leaves @c matrixIsSetup untouched if only columns have been added, since the 305 * basis matrix already has the correct size. */ 306 if(status() > NO_PROBLEM && matrixIsSetup) 307 loadMatrixVecs(); 308 309 switch(status()) 310 { 311 case DUAL: 312 case INFEASIBLE: 313 setStatus(REGULAR); 314 break; 315 316 case OPTIMAL: 317 case UNBOUNDED: 318 setStatus(PRIMAL); 319 break; 320 321 case NO_PROBLEM: 322 case SINGULAR: 323 case REGULAR: 324 case PRIMAL: 325 break; 326 327 default: 328 SPX_MSG_ERROR(std::cerr << "ECHBAS08 Unknown basis status!" << std::endl;) 329 throw SPxInternalCodeException("XCHBAS02 This should never happen."); 330 } 331 } 332 } 333 334 template <class R> 335 void SPxBasisBase<R>::removedCol(int i) 336 { 337 assert(status() > NO_PROBLEM); 338 assert(theLP != 0); 339 340 if(theLP->rep() == SPxSolverBase<R>::COLUMN) 341 { 342 if(theLP->isBasic(thedesc.colStatus(i))) 343 setStatus(NO_PROBLEM); 344 } 345 else 346 { 347 assert(theLP->rep() == SPxSolverBase<R>::ROW); 348 factorized = false; 349 350 if(!theLP->isBasic(thedesc.colStatus(i))) 351 setStatus(NO_PROBLEM); 352 else if(status() > NO_PROBLEM) 353 { 354 for(int j = theLP->dim(); j >= 0; --j) 355 { 356 SPxId id = baseId(j); 357 358 if(id.isSPxColId() && !theLP->has(SPxColId(id))) 359 { 360 baseId(j) = baseId(theLP->dim()); 361 362 if(matrixIsSetup && 363 j < theLP->dim()) 364 matrix[j] = &theLP->vector(baseId(j)); 365 366 break; 367 } 368 } 369 } 370 } 371 372 thedesc.colStatus(i) = thedesc.colStatus(theLP->nCols()); 373 reDim(); 374 } 375 376 template <class R> 377 void SPxBasisBase<R>::removedCols(const int perm[]) 378 { 379 assert(status() > NO_PROBLEM); 380 assert(theLP != 0); 381 382 int i; 383 int n = thedesc.nCols(); 384 385 if(theLP->rep() == SPxSolverBase<R>::COLUMN) 386 { 387 for(i = 0; i < n; ++i) 388 { 389 if(perm[i] < 0) // column got removed 390 { 391 if(theLP->isBasic(thedesc.colStatus(i))) 392 setStatus(NO_PROBLEM); 393 } 394 else // column was potentially moved 395 thedesc.colStatus(perm[i]) = thedesc.colStatus(i); 396 } 397 } 398 else 399 { 400 assert(theLP->rep() == SPxSolverBase<R>::ROW); 401 factorized = matrixIsSetup = false; 402 403 for(i = 0; i < n; ++i) 404 { 405 if(perm[i] != i) 406 { 407 if(perm[i] < 0) // column got removed 408 { 409 if(!theLP->isBasic(thedesc.colStatus(i))) 410 setStatus(NO_PROBLEM); 411 } 412 else // column was moved 413 thedesc.colStatus(perm[i]) = thedesc.colStatus(i); 414 } 415 } 416 } 417 418 reDim(); 419 } 420 421 422 /** 423 * mark the basis as not factorized 424 */ 425 template <class R> 426 void SPxBasisBase<R>::invalidate() 427 { 428 if(factorized || matrixIsSetup) 429 { 430 SPX_MSG_INFO3((*this->spxout), 431 (*this->spxout) << "ICHBAS09 explicit invalidation of factorization" << 432 std::endl;) 433 } 434 435 factorized = false; 436 matrixIsSetup = false; 437 } 438 439 /** 440 * Create the initial slack basis descriptor and set up the basis matrix accordingly. 441 * This code has been adapted from SPxBasisBase<R>::addedRows() and SPxBasisBase<R>::addedCols(). 442 */ 443 template <class R> 444 void SPxBasisBase<R>::restoreInitialBasis() 445 { 446 assert(!factorized); 447 448 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "ICHBAS10 setup slack basis" << std::endl;) 449 450 if(theLP->rep() == SPxSolverBase<R>::COLUMN) 451 { 452 for(int i = 0; i < theLP->nRows(); ++i) 453 { 454 thedesc.rowStatus(i) = dualRowStatus(i); 455 baseId(i) = theLP->SPxLPBase<R>::rId(i); 456 } 457 458 for(int i = 0; i < theLP->nCols(); ++i) 459 thedesc.colStatus(i) = primalColStatus(i, theLP); 460 } 461 else 462 { 463 assert(theLP->rep() == SPxSolverBase<R>::ROW); 464 465 for(int i = 0; i < theLP->nRows(); ++i) 466 thedesc.rowStatus(i) = dualRowStatus(i); 467 468 for(int i = 0; i < theLP->nCols(); ++i) 469 { 470 thedesc.colStatus(i) = primalColStatus(i, theLP); 471 baseId(i) = theLP->SPxLPBase<R>::cId(i); 472 } 473 } 474 475 /* if matrix was set up, load new basis vectors to the matrix */ 476 if(status() > NO_PROBLEM && matrixIsSetup) 477 loadMatrixVecs(); 478 479 /* update basis status */ 480 setStatus(REGULAR); 481 } 482 483 /** 484 * @todo Implement changedRow(), changedCol(), changedElement() in a more clever 485 * way. For instance, the basis won't be singular (but maybe infeasible) if the 486 * change doesn't affect the basis rows/columns. 487 * 488 * The following methods (changedRow(), changedCol(), changedElement()) radically 489 * change the current basis to the original (slack) basis also present after 490 * loading the LP. The reason is that through the changes, the current basis may 491 * become singular. Going back to the initial basis is quite inefficient, but 492 * correct. 493 */ 494 495 /**@todo is this correctly implemented? 496 */ 497 template <class R> 498 void SPxBasisBase<R>::changedRow(int /*row*/) 499 { 500 invalidate(); 501 restoreInitialBasis(); 502 } 503 504 /**@todo is this correctly implemented? 505 */ 506 template <class R> 507 void SPxBasisBase<R>::changedCol(int /*col*/) 508 { 509 invalidate(); 510 restoreInitialBasis(); 511 } 512 513 /**@todo is this correctly implemented? 514 */ 515 template <class R> 516 void SPxBasisBase<R>::changedElement(int /*row*/, int /*col*/) 517 { 518 invalidate(); 519 restoreInitialBasis(); 520 } 521 } // namespace soplex 522