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 <assert.h> 26 #include <cstdio> 27 #include <iostream> 28 #include <iomanip> 29 #include <sstream> 30 31 #include "soplex/spxdefines.h" 32 #include "soplex/didxset.h" 33 #include "soplex/mpsinput.h" 34 #include "soplex/spxout.h" 35 #include "soplex/exceptions.h" 36 37 namespace soplex 38 { 39 40 template <class R> typename SPxBasisBase<R>::Desc::Status 41 SPxBasisBase<R>::dualStatus(const SPxColId& id) const 42 { 43 return dualColStatus(static_cast<SPxLPBase<R>*>(theLP)->number(id)); 44 } 45 46 template <class R> 47 typename SPxBasisBase<R>::Desc::Status 48 SPxBasisBase<R>::dualStatus(const SPxRowId& id) const 49 { 50 return dualRowStatus((static_cast<SPxLPBase<R>*>(theLP))->number(id)); 51 } 52 53 template <class R> 54 typename SPxBasisBase<R>::Desc::Status 55 SPxBasisBase<R>::dualRowStatus(int i) const 56 { 57 assert(theLP != 0); 58 59 if(theLP->rhs(i) < R(infinity)) 60 { 61 if(theLP->lhs(i) > R(-infinity)) 62 { 63 if(theLP->lhs(i) == theLP->rhs(i)) 64 return Desc::D_FREE; 65 else 66 return Desc::D_ON_BOTH; 67 } 68 else 69 return Desc::D_ON_LOWER; 70 } 71 else if(theLP->lhs(i) > R(-infinity)) 72 return Desc::D_ON_UPPER; 73 else 74 return Desc::D_UNDEFINED; 75 } 76 77 template <class R> 78 typename SPxBasisBase<R>::Desc::Status 79 SPxBasisBase<R>::dualColStatus(int i) const 80 { 81 assert(theLP != 0); 82 83 if(theLP->SPxLPBase<R>::upper(i) < R(infinity)) 84 { 85 if(theLP->SPxLPBase<R>::lower(i) > R(-infinity)) 86 { 87 if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i)) 88 return Desc::D_FREE; 89 else 90 return Desc::D_ON_BOTH; 91 } 92 else 93 return Desc::D_ON_LOWER; 94 } 95 else if(theLP->SPxLPBase<R>::lower(i) > R(-infinity)) 96 return Desc::D_ON_UPPER; 97 else 98 return Desc::D_UNDEFINED; 99 } 100 101 template <class R> 102 void SPxBasisBase<R>::loadMatrixVecs() 103 { 104 assert(theLP != 0); 105 assert(theLP->dim() == matrix.size()); 106 107 SPX_MSG_INFO3((*this->spxout), 108 (*this->spxout) << "IBASIS01 loadMatrixVecs() invalidates factorization" 109 << std::endl;) 110 111 int i; 112 nzCount = 0; 113 114 for(i = theLP->dim() - 1; i >= 0; --i) 115 { 116 matrix[i] = &theLP->vector(baseId(i)); 117 nzCount += matrix[i]->size(); 118 } 119 120 matrixIsSetup = true; 121 factorized = false; 122 123 if(factor != 0) 124 factor->clear(); 125 } 126 127 template <class R> 128 bool SPxBasisBase<R>::isDescValid(const Desc& ds) 129 { 130 131 assert(status() > NO_PROBLEM); 132 assert(theLP != 0); 133 134 int basisdim; 135 136 if(ds.nRows() != theLP->nRows() || ds.nCols() != theLP->nCols()) 137 { 138 SPxOut::debug(this, "IBASIS20 Dimension mismatch\n"); 139 return false; 140 } 141 142 basisdim = 0; 143 144 for(int row = ds.nRows() - 1; row >= 0; --row) 145 { 146 if(ds.rowstat[row] >= 0) 147 { 148 if(ds.rowstat[row] != dualRowStatus(row)) 149 { 150 SPxOut::debug(this, "IBASIS21 Basic row {} with incorrect dual status {}\n", row, 151 ds.rowstat[row]); 152 return false; 153 } 154 } 155 else 156 { 157 basisdim++; 158 159 if((ds.rowstat[row] == Desc::P_FIXED 160 && theLP->SPxLPBase<R>::lhs(row) != theLP->SPxLPBase<R>::rhs(row)) 161 || (ds.rowstat[row] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::rhs(row) >= R(infinity)) 162 || (ds.rowstat[row] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lhs(row) <= R(-infinity))) 163 { 164 SPxOut::debug(this, "IBASIS22 Nonbasic row with incorrect status: lhs={}, rhs={}, stat={}\n", 165 theLP->SPxLPBase<R>::lhs(row), theLP->SPxLPBase<R>::rhs(row), 166 ds.rowstat[row]); 167 return false; 168 } 169 } 170 } 171 172 for(int col = ds.nCols() - 1; col >= 0; --col) 173 { 174 if(ds.colstat[col] >= 0) 175 { 176 if(ds.colstat[col] != dualColStatus(col)) 177 { 178 SPxOut::debug(this, "IBASIS23 Basic column {} with incorrect dual status {}\n", col, 179 ds.colstat[col]); 180 return false; 181 } 182 } 183 else 184 { 185 basisdim++; 186 187 if((ds.colstat[col] == Desc::P_FIXED 188 && theLP->SPxLPBase<R>::lower(col) != theLP->SPxLPBase<R>::upper(col)) 189 || (ds.colstat[col] == Desc::P_ON_UPPER && theLP->SPxLPBase<R>::upper(col) >= R(infinity)) 190 || (ds.colstat[col] == Desc::P_ON_LOWER && theLP->SPxLPBase<R>::lower(col) <= R(-infinity))) 191 { 192 SPxOut::debug(this, 193 "IBASIS24 Nonbasic column {} with incorrect status: lower={}, upper={}, stat={}\n", 194 col, theLP->SPxLPBase<R>::lower(col), theLP->SPxLPBase<R>::upper(col), 195 ds.colstat[col]); 196 return false; 197 } 198 } 199 } 200 201 if(basisdim != theLP->nCols()) 202 { 203 SPxOut::debug(this, "IBASIS25 Incorrect basis dimension {} != {}\n", basisdim, theLP->nCols()); 204 return false; 205 } 206 207 // basis descriptor valid 208 return true; 209 } 210 211 212 /* 213 Loading a #Desc# into the basis can be done more efficiently, by 214 explicitely programming both cases, for the rowwise and for the columnwise 215 representation. This implementation hides this distinction in the use of 216 methods #isBasic()# and #vector()#. 217 */ 218 template <class R> 219 void SPxBasisBase<R>::loadDesc(const Desc& ds) 220 { 221 assert(status() > NO_PROBLEM); 222 assert(theLP != 0); 223 assert(ds.nRows() == theLP->nRows()); 224 assert(ds.nCols() == theLP->nCols()); 225 226 SPxId none; 227 int i; 228 int j; 229 bool consistent = true; 230 231 SPX_MSG_INFO3((*this->spxout), 232 (*this->spxout) << "IBASIS02 loading of Basis invalidates factorization" 233 << std::endl;) 234 235 lastin = none; 236 lastout = none; 237 lastidx = -1; 238 iterCount = 0; 239 updateCount = 0; 240 241 if(&ds != &thedesc) 242 { 243 thedesc = ds; 244 setRep(); 245 } 246 247 assert(theLP->dim() == matrix.size()); 248 249 nzCount = 0; 250 251 for(j = i = 0; i < theLP->nRows(); ++i) 252 { 253 /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis 254 * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them 255 */ 256 if(thedesc.rowStatus(i) >= 0) 257 thedesc.rowStatus(i) = dualRowStatus(i); 258 else if(thedesc.rowStatus(i) == SPxBasisBase<R>::Desc::P_FIXED 259 && theLP->SPxLPBase<R>::lhs(i) != theLP->SPxLPBase<R>::rhs(i)) 260 { 261 if(theLP->SPxLPBase<R>::lhs(i) > R(-infinity) && theLP->SPxLPBase<R>::maxRowObj(i) < 0.0) 262 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER; 263 else if(theLP->SPxLPBase<R>::rhs(i) < R(infinity)) 264 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER; 265 else 266 thedesc.rowStatus(i) = SPxBasisBase<R>::Desc::P_FREE; 267 } 268 269 if(theLP->isBasic(thedesc.rowStatus(i))) 270 { 271 assert(theLP->dim() == matrix.size()); 272 assert(j <= matrix.size()); 273 274 if(j == matrix.size()) 275 { 276 // too many basic variables 277 consistent = false; 278 break; 279 } 280 281 SPxRowId id = theLP->SPxLPBase<R>::rId(i); 282 theBaseId[j] = id; 283 matrix[j] = &theLP->vector(id); 284 nzCount += matrix[j++]->size(); 285 } 286 } 287 288 for(i = 0; i < theLP->nCols(); ++i) 289 { 290 /* for columns and rows with D_... status, the correct D_... status depends on bounds and sides; if a basis 291 * descriptor is loaded after changing bounds or sides, e.g. in the refine() method, we have to correct them 292 */ 293 if(thedesc.colStatus(i) >= 0) 294 thedesc.colStatus(i) = dualColStatus(i); 295 else if(thedesc.colStatus(i) == SPxBasisBase<R>::Desc::P_FIXED 296 && theLP->SPxLPBase<R>::lower(i) != theLP->SPxLPBase<R>::upper(i)) 297 { 298 if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity) && theLP->SPxLPBase<R>::upper(i) >= R(infinity)) 299 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_FREE; 300 else if(theLP->SPxLPBase<R>::upper(i) >= R(infinity) 301 || (theLP->SPxLPBase<R>::lower(i) > R(-infinity) && theLP->SPxLPBase<R>::maxObj(i) < 0.0)) 302 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER; 303 else 304 thedesc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER; 305 } 306 307 if(theLP->isBasic(thedesc.colStatus(i))) 308 { 309 assert(theLP->dim() == matrix.size()); 310 assert(j <= matrix.size()); 311 312 if(j == matrix.size()) 313 { 314 // too many basic variables 315 consistent = false; 316 break; 317 } 318 319 SPxColId id = theLP->SPxLPBase<R>::cId(i); 320 theBaseId[j] = id; 321 matrix[j] = &theLP->vector(id); 322 nzCount += matrix[j++]->size(); 323 } 324 } 325 326 if(j < matrix.size()) 327 consistent = false; // not enough basic variables 328 329 /* if dimensions are inconsistent, restore slack basis 330 * if dimensions are consistent, then we have setup a correct matrix 331 */ 332 if(!consistent) 333 restoreInitialBasis(); 334 else 335 matrixIsSetup = true; 336 337 assert(isDescValid(thedesc)); 338 339 factorized = false; 340 341 if(factor != 0) 342 factor->clear(); 343 } 344 345 template <class R> 346 void SPxBasisBase<R>::setRep() 347 { 348 assert(theLP != 0); 349 350 reDim(); 351 minStab = 0.0; 352 353 if(theLP->rep() == SPxSolverBase<R>::ROW) 354 { 355 thedesc.stat = &thedesc.rowstat; 356 thedesc.costat = &thedesc.colstat; 357 } 358 else 359 { 360 thedesc.stat = &thedesc.colstat; 361 thedesc.costat = &thedesc.rowstat; 362 } 363 } 364 365 template <class R> 366 void SPxBasisBase<R>::load(SPxSolverBase<R>* lp, bool initSlackBasis) 367 { 368 assert(lp != 0); 369 theLP = lp; 370 371 setOutstream(*theLP->spxout); 372 373 setRep(); 374 375 if(initSlackBasis) 376 { 377 restoreInitialBasis(); 378 loadDesc(thedesc); 379 } 380 } 381 382 template <class R> 383 void SPxBasisBase<R>::loadBasisSolver(SLinSolver<R>* p_solver, const bool destroy) 384 { 385 assert(!freeSlinSolver || factor != 0); 386 387 setOutstream(*p_solver->spxout); 388 389 SPX_MSG_INFO3((*this->spxout), 390 (*this->spxout) << "IBASIS03 loading of Solver invalidates factorization" 391 << std::endl;) 392 393 if(freeSlinSolver) 394 { 395 delete factor; 396 factor = 0; 397 } 398 399 factor = p_solver; 400 factorized = false; 401 factor->clear(); 402 freeSlinSolver = destroy; 403 } 404 405 /** 406 * The specification is taken from the 407 * 408 * ILOG CPLEX 7.0 Reference Manual, Appendix E, Page 543. 409 * 410 * This routine should read valid BAS format files. 411 * 412 * @return true if the file was read correctly. 413 * 414 * Here is a very brief outline of the format: 415 * 416 * The format is in a form similar to an MPS file. The basic assumption is that all (column) 417 * variables are nonbasic at their lower bound and all row (variables) are basic; only the 418 * differences to this rule are given. Each data line contains an indicator, a variable name and 419 * possibly a row/constraint name. The following meaning applies with respect to the indicators: 420 * 421 * - XU: the variable is basic, the row is nonbasic at its upper bound 422 * - XL: the variable is basic, the row is nonbasic at its lower bound 423 * - UL: the variable is nonbasic and at its upper bound 424 * - LL: the variable is nonbasic and at its lower bound 425 * 426 * The CPLEX format contains an additional indicator 'BS', but this is unsupported here. 427 * 428 * Nonbasic variables without lower bound have the following default status for SoPlex: 429 * - at their upper bound if finite, 430 * - at zero if free. 431 */ 432 template <class R> 433 bool SPxBasisBase<R>::readBasis( 434 std::istream& is, 435 const NameSet* rowNames, 436 const NameSet* colNames) 437 { 438 assert(theLP != 0); 439 440 /* prepare names */ 441 const NameSet* rNames = rowNames; 442 const NameSet* cNames = colNames; 443 444 NameSet* p_colNames = 0; 445 NameSet* p_rowNames = 0; 446 447 if(colNames == 0) 448 { 449 int nCols = theLP->nCols(); 450 std::stringstream name; 451 452 spx_alloc(p_colNames); 453 p_colNames = new(p_colNames) NameSet(); 454 p_colNames->reMax(nCols); 455 456 for(int j = 0; j < nCols; ++j) 457 { 458 name << "x" << j; 459 DataKey key = theLP->colId(j); 460 p_colNames->add(key, name.str().c_str()); 461 } 462 463 cNames = p_colNames; 464 } 465 466 if(rNames == 0) 467 { 468 int nRows = theLP->nRows(); 469 std::stringstream name; 470 471 spx_alloc(p_rowNames); 472 p_rowNames = new(p_rowNames) NameSet(); 473 p_rowNames->reMax(nRows); 474 475 for(int i = 0; i < nRows; ++i) 476 { 477 name << "C" << i; 478 DataKey key = theLP->rowId(i); 479 p_rowNames->add(key, name.str().c_str()); 480 } 481 482 rNames = p_rowNames; 483 } 484 485 /* load default basis if necessary */ 486 if(status() == NO_PROBLEM) 487 load(theLP, false); 488 489 /* initialize with standard settings */ 490 Desc l_desc(thedesc); 491 492 for(int i = 0; i < theLP->nRows(); i++) 493 l_desc.rowstat[i] = dualRowStatus(i); 494 495 for(int i = 0; i < theLP->nCols(); i++) 496 { 497 if(theLP->SPxLPBase<R>::lower(i) == theLP->SPxLPBase<R>::upper(i)) 498 l_desc.colstat[i] = Desc::P_FIXED; 499 else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity) 500 && theLP->SPxLPBase<R>::upper(i) >= R(infinity)) 501 l_desc.colstat[i] = Desc::P_FREE; 502 else if(theLP->SPxLPBase<R>::lower(i) <= R(-infinity)) 503 l_desc.colstat[i] = Desc::P_ON_UPPER; 504 else 505 l_desc.colstat[i] = Desc::P_ON_LOWER; 506 } 507 508 MPSInput mps(is); 509 510 if(mps.readLine() && (mps.field0() != 0) && !strcmp(mps.field0(), "NAME")) 511 { 512 while(mps.readLine()) 513 { 514 int c = -1; 515 int r = -1; 516 517 if((mps.field0() != 0) && !strcmp(mps.field0(), "ENDATA")) 518 { 519 mps.setSection(MPSInput::ENDATA); 520 break; 521 } 522 523 if((mps.field1() == 0) || (mps.field2() == 0)) 524 break; 525 526 if((c = cNames->number(mps.field2())) < 0) 527 break; 528 529 if(*mps.field1() == 'X') 530 if(mps.field3() == 0 || (r = rNames->number(mps.field3())) < 0) 531 break; 532 533 if(!strcmp(mps.field1(), "XU")) 534 { 535 l_desc.colstat[c] = dualColStatus(c); 536 537 if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::GREATER_EQUAL) 538 l_desc.rowstat[r] = Desc::P_ON_LOWER; 539 else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL) 540 l_desc.rowstat[r] = Desc::P_FIXED; 541 else 542 l_desc.rowstat[r] = Desc::P_ON_UPPER; 543 } 544 else if(!strcmp(mps.field1(), "XL")) 545 { 546 l_desc.colstat[c] = dualColStatus(c); 547 548 if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::LESS_EQUAL) 549 l_desc.rowstat[r] = Desc::P_ON_UPPER; 550 else if(theLP->LPRowSetBase<R>::type(r) == LPRowBase<R>::EQUAL) 551 l_desc.rowstat[r] = Desc::P_FIXED; 552 else 553 l_desc.rowstat[r] = Desc::P_ON_LOWER; 554 } 555 else if(!strcmp(mps.field1(), "UL")) 556 { 557 l_desc.colstat[c] = Desc::P_ON_UPPER; 558 } 559 else if(!strcmp(mps.field1(), "LL")) 560 { 561 l_desc.colstat[c] = Desc::P_ON_LOWER; 562 } 563 else 564 { 565 mps.syntaxError(); 566 break; 567 } 568 } 569 } 570 571 if(!mps.hasError()) 572 { 573 if(mps.section() == MPSInput::ENDATA) 574 { 575 // force basis to be different from NO_PROBLEM 576 // otherwise the basis will be overwritten at later stages. 577 setStatus(REGULAR); 578 loadDesc(l_desc); 579 } 580 else 581 mps.syntaxError(); 582 } 583 584 if(rowNames == 0) 585 { 586 p_rowNames->~NameSet(); 587 spx_free(p_rowNames); 588 } 589 590 if(colNames == 0) 591 { 592 p_colNames->~NameSet(); 593 spx_free(p_colNames); 594 } 595 596 #ifndef NDEBUG 597 SPX_DEBUG(thedesc.dump()); 598 #endif 599 600 return !mps.hasError(); 601 } 602 603 604 /* Get row name - copied from spxmpswrite.cpp 605 * 606 * @todo put this in a common file and unify accross different formats (mps, lp, basis). 607 */ 608 template <class R> 609 static const char* getRowName( 610 const SPxLPBase<R>* lp, 611 int idx, 612 const NameSet* rnames, 613 char* buf) 614 { 615 assert(buf != 0); 616 assert(idx >= 0); 617 assert(idx < lp->nRows()); 618 619 if(rnames != 0) 620 { 621 DataKey key = lp->rId(idx); 622 623 if(rnames->has(key)) 624 return (*rnames)[key]; 625 } 626 627 spxSnprintf(buf, 16, "C%d", idx); 628 629 return buf; 630 } 631 632 /* Get column name - copied from spxmpswrite.cpp 633 * 634 * @todo put this in a common file and unify accross different formats (mps, lp, basis). 635 */ 636 template <class R> 637 static const char* getColName( 638 const SPxLPBase<R>* lp, 639 int idx, 640 const NameSet* cnames, 641 char* buf) 642 { 643 assert(buf != 0); 644 assert(idx >= 0); 645 assert(idx < lp->nCols()); 646 647 if(cnames != 0) 648 { 649 DataKey key = lp->cId(idx); 650 651 if(cnames->has(key)) 652 return (*cnames)[key]; 653 } 654 655 spxSnprintf(buf, 16, "x%d", idx); 656 657 return buf; 658 } 659 660 /* writes a file in MPS basis format to \p os. 661 * 662 * See SPxBasisBase<R>::readBasis() for a short description of the format. 663 */ 664 template <class R> 665 void SPxBasisBase<R>::writeBasis( 666 std::ostream& os, 667 const NameSet* rowNames, 668 const NameSet* colNames, 669 const bool cpxFormat 670 ) const 671 { 672 assert(theLP != 0); 673 674 os.setf(std::ios::left); 675 os << "NAME soplex.bas\n"; 676 677 /* do not write basis if there is none */ 678 if(status() == NO_PROBLEM) 679 { 680 os << "ENDATA" << std::endl; 681 return; 682 } 683 684 /* start writing */ 685 char buf[255]; 686 int row = 0; 687 688 for(int col = 0; col < theLP->nCols(); col++) 689 { 690 if(thedesc.colStatus(col) > 0) 691 { 692 /* Find non basic row */ 693 for(; row < theLP->nRows(); row++) 694 { 695 if(thedesc.rowStatus(row) < 0) 696 break; 697 } 698 699 assert(row != theLP->nRows()); 700 701 if(thedesc.rowStatus(row) == Desc::P_ON_UPPER && (!cpxFormat 702 || theLP->LPRowSetBase<R>::type(row) == LPRowBase<R>::RANGE)) 703 os << " XU "; 704 else 705 os << " XL "; 706 707 os << std::setw(8) << getColName(theLP, col, colNames, buf); 708 709 /* break in two parts since buf is reused */ 710 os << " " 711 << getRowName(theLP, row, rowNames, buf) 712 << std::endl; 713 714 row++; 715 } 716 else 717 { 718 if(thedesc.colStatus(col) == Desc::P_ON_UPPER) 719 { 720 os << " UL " 721 << getColName(theLP, col, colNames, buf) 722 << std::endl; 723 } 724 else 725 { 726 /* Default is all non-basic variables on lower bound (if finite) or at zero (if free). 727 * nothing to do in this case. 728 */ 729 assert(theLP->lower(col) <= R(-infinity) || thedesc.colStatus(col) == Desc::P_ON_LOWER 730 || thedesc.colStatus(col) == Desc::P_FIXED); 731 assert(theLP->lower(col) > R(-infinity) || theLP->upper(col) < R(infinity) 732 || thedesc.colStatus(col) == Desc::P_FREE); 733 } 734 } 735 } 736 737 #ifndef NDEBUG 738 SPX_DEBUG(thedesc.dump()); 739 740 // Check that we covered all nonbasic rows - the remaining should be basic. 741 for(; row < theLP->nRows(); row++) 742 { 743 if(thedesc.rowStatus(row) < 0) 744 break; 745 } 746 747 assert(row == theLP->nRows()); 748 749 #endif // NDEBUG 750 751 os << "ENDATA" << std::endl; 752 } 753 754 template <class R> 755 void SPxBasisBase<R>::printMatrix() const 756 { 757 758 assert(matrixIsSetup); 759 760 for(int i = 0; i < matrix.size(); i++) 761 { 762 std::cout << "C" << i << "=" << *matrix[i] << std::endl; 763 } 764 } 765 766 template <class R> 767 void SPxBasisBase<R>::printMatrixMTX(int number) 768 { 769 int dim; 770 int nnz; 771 char filename[SPX_MAXSTRLEN]; 772 773 dim = matrix.size(); 774 nnz = nzCount; 775 spxSnprintf(filename, SPX_MAXSTRLEN, "basis/basis%d.mtx", number); 776 std::cout << "printing basis matrix to file " << filename << "\n"; 777 FILE* basisfile; 778 basisfile = fopen(filename, "w"); 779 // print marker necessary for reading the file in Matlab 780 fprintf(basisfile, "%%%%MatrixMarket matrix coordinate Real general\n"); 781 // print matrix information 782 fprintf(basisfile, "%d %d %d\n", dim, dim, nnz); 783 784 // print matrix data 785 for(int i = 0; i < matrix.size(); ++i) 786 { 787 for(int j = 0; j < baseVec(i).size(); ++j) 788 { 789 int idx = baseVec(i).index(j); 790 R val = baseVec(i).value(j); 791 fprintf(basisfile, "%d %d %.13" SOPLEX_REAL_FORMAT "\n", i + 1, idx + 1, val); 792 } 793 } 794 795 fclose(basisfile); 796 797 return; 798 } 799 800 template <class R> 801 void SPxBasisBase<R>::change( 802 int i, 803 SPxId& id, 804 const SVectorBase<R>* enterVec, 805 const SSVectorBase<R>* eta) 806 { 807 808 assert(matrixIsSetup); 809 assert(!id.isValid() || (enterVec != 0)); 810 assert(factor != 0); 811 812 lastidx = i; 813 lastin = id; 814 815 if(id.isValid() && i >= 0) 816 { 817 assert(enterVec != 0); 818 819 // update the counter for nonzeros in the basis matrix 820 nzCount = nzCount - matrix[i]->size() + enterVec->size(); 821 // let the new id enter the basis 822 matrix[i] = enterVec; 823 lastout = theBaseId[i]; 824 theBaseId[i] = id; 825 826 ++iterCount; 827 ++updateCount; 828 829 SPxOut::debug(this, 830 "factor_stats: iteration= {}, update= {}, total_update= {}, nonzero_B= {}, nonzero_LU= {}, factor_fill= {}, time= {}\n", 831 this->iteration(), updateCount, totalUpdateCount, nzCount, factor->memory(), lastFill, 832 theLP->time()); 833 834 // never factorize? Just do it ! 835 if(!factorized) 836 factorize(); 837 838 // too much memory growth ? 839 else if(R(factor->memory()) > 1000 + factor->dim() + lastMem * memFactor) 840 { 841 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 842 "IBASIS04 memory growth factor triggers refactorization" 843 << " memory= " << factor->memory() 844 << " lastMem= " << lastMem 845 << " memFactor= " << memFactor 846 << std::endl;) 847 factorize(); 848 } 849 850 // relative fill too high ? 851 else if(R(factor->memory()) > lastFill * R(nzCount)) 852 { 853 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS04 fill factor triggers refactorization" 854 << " memory= " << factor->memory() 855 << " nzCount= " << nzCount 856 << " lastFill= " << lastFill 857 << std::endl;) 858 859 factorize(); 860 } 861 // absolute fill in basis matrix too high ? 862 else if(nzCount > lastNzCount) 863 { 864 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS05 nonzero factor triggers refactorization" 865 << " nzCount= " << nzCount 866 << " lastNzCount= " << lastNzCount 867 << " nonzeroFactor= " << nonzeroFactor 868 << std::endl;) 869 factorize(); 870 } 871 // too many updates ? 872 else if(updateCount >= maxUpdates) 873 { 874 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS06 update count triggers refactorization" 875 << " updateCount= " << updateCount 876 << " maxUpdates= " << maxUpdates 877 << std::endl;) 878 factorize(); 879 } 880 else 881 { 882 try 883 { 884 #ifdef MEASUREUPDATETIME 885 theTime.start(); 886 #endif 887 factor->change(i, *enterVec, eta); 888 totalUpdateCount++; 889 #ifdef MEASUREUPDATETIME 890 theTime.stop(); 891 #endif 892 } 893 catch(...) 894 { 895 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 896 "IBASIS13 problems updating factorization; refactorizing basis" 897 << std::endl;) 898 899 #ifdef MEASUREUPDATETIME 900 theTime.stop(); 901 #endif 902 903 // singularity was detected in update; we refactorize 904 factorize(); 905 906 // if factorize() detects singularity, an exception is thrown, hence at this point we have a regular basis 907 // and can try the update again 908 assert(status() >= SPxBasisBase<R>::REGULAR); 909 910 try 911 { 912 #ifdef MEASUREUPDATETIME 913 theTime.start(); 914 #endif 915 factor->change(i, *enterVec, eta); 916 totalUpdateCount++; 917 #ifdef MEASUREUPDATETIME 918 theTime.stop(); 919 #endif 920 } 921 // with a freshly factorized, regular basis, the update is unlikely to fail; if this happens nevertheless, 922 // we have to invalidate the basis to have the statuses correct 923 catch(const SPxException& F) 924 { 925 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 926 "IBASIS14 problems updating factorization; invalidating factorization" 927 << std::endl;) 928 929 #ifdef MEASUREUPDATETIME 930 theTime.stop(); 931 #endif 932 933 factorized = false; 934 throw F; 935 } 936 } 937 938 assert(minStab > 0.0); 939 940 if(factor->status() != SLinSolver<R>::OK || factor->stability() < minStab) 941 { 942 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << "IBASIS07 stability triggers refactorization" 943 << " stability= " << factor->stability() 944 << " minStab= " << minStab 945 << std::endl;) 946 factorize(); 947 } 948 } 949 } 950 else 951 lastout = id; 952 } 953 954 template <class R> 955 void SPxBasisBase<R>::factorize() 956 { 957 958 assert(factor != 0); 959 960 if(!matrixIsSetup) 961 loadDesc(thedesc); 962 963 assert(matrixIsSetup); 964 965 updateCount = 0; 966 967 switch(factor->load(matrix.get_ptr(), matrix.size())) 968 { 969 case SLinSolver<R>::OK : 970 if(status() == SINGULAR) 971 setStatus(REGULAR); 972 973 factorized = true; 974 minStab = factor->stability(); 975 976 // This seems always to be about 1e-7 977 if(minStab > 1e-4) 978 minStab *= 0.001; 979 980 if(minStab > 1e-5) 981 minStab *= 0.01; 982 983 if(minStab > 1e-6) 984 minStab *= 0.1; 985 986 break; 987 988 case SLinSolver<R>::SINGULAR : 989 setStatus(SINGULAR); 990 factorized = false; 991 break; 992 993 default : 994 SPX_MSG_ERROR(std::cerr << "EBASIS08 error: unknown status of factorization.\n";) 995 factorized = false; 996 throw SPxInternalCodeException("XBASIS01 This should never happen."); 997 } 998 999 // get nonzero count of factorization 1000 lastMem = factor->memory(); 1001 // compute fill ratio between factorization and basis matrix (multiplied with tolerance) 1002 lastFill = fillFactor * R(lastMem) / R(nzCount > 0 ? nzCount : 1); 1003 lastNzCount = int(nonzeroFactor * R(nzCount > 0 ? nzCount : 1)); 1004 1005 if(status() == SINGULAR) 1006 { 1007 throw SPxStatusException("Cannot factorize singular matrix"); 1008 } 1009 } 1010 1011 template <class R> 1012 VectorBase<R>& SPxBasisBase<R>::multWithBase(VectorBase<R>& x) const 1013 { 1014 assert(status() > SINGULAR); 1015 assert(theLP->dim() == x.dim()); 1016 1017 int i; 1018 VectorBase<R> tmp(x); 1019 1020 if(!matrixIsSetup) 1021 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc); 1022 1023 assert(matrixIsSetup); 1024 1025 for(i = x.dim() - 1; i >= 0; --i) 1026 x[i] = *(matrix[i]) * tmp; 1027 1028 return x; 1029 } 1030 1031 template <class R> 1032 void SPxBasisBase<R>::multWithBase(SSVectorBase<R>& x, SSVectorBase<R>& result) const 1033 { 1034 assert(status() > SINGULAR); 1035 assert(theLP->dim() == x.dim()); 1036 assert(x.dim() == result.dim()); 1037 1038 if(!matrixIsSetup) 1039 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc); 1040 1041 result.clear(); 1042 1043 assert(matrixIsSetup); 1044 1045 for(int i = 0; i < x.dim(); ++i) 1046 result.add(i, (*matrix[i]) * x); 1047 1048 return; 1049 } 1050 1051 template <class R> 1052 VectorBase<R>& SPxBasisBase<R>::multBaseWith(VectorBase<R>& x) const 1053 { 1054 assert(status() > SINGULAR); 1055 assert(theLP->dim() == x.dim()); 1056 1057 int i; 1058 VectorBase<R> tmp(x); 1059 1060 if(!matrixIsSetup) 1061 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc); 1062 1063 assert(matrixIsSetup); 1064 1065 x.clear(); 1066 1067 for(i = x.dim() - 1; i >= 0; --i) 1068 { 1069 if(tmp[i] != 0.0) 1070 x.multAdd(tmp[i], *(matrix[i])); 1071 } 1072 1073 return x; 1074 } 1075 1076 template <class R> 1077 void SPxBasisBase<R>::multBaseWith(SSVectorBase<R>& x, SSVectorBase<R>& result) const 1078 { 1079 assert(status() > SINGULAR); 1080 assert(theLP->dim() == x.dim()); 1081 assert(x.dim() == result.dim()); 1082 1083 if(!matrixIsSetup) 1084 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc); 1085 1086 assert(matrixIsSetup); 1087 1088 result.clear(); 1089 1090 if(x.isSetup()) 1091 { 1092 for(int i = 0; i < x.size(); ++i) 1093 { 1094 int idx = x.index(i); 1095 result.multAdd(x[idx], (*matrix[idx])); 1096 } 1097 } 1098 else 1099 { 1100 for(int i = 0; i < x.dim(); ++i) 1101 result.multAdd(x[i], (*matrix[i])); 1102 } 1103 1104 return; 1105 } 1106 1107 template <class R> 1108 /* compute an estimated condition number for the current basis matrix 1109 * by computing estimates of the norms of B and B^-1 using the power method. 1110 * maxiters and tolerance control the accuracy of the estimate. 1111 */ 1112 R SPxBasisBase<R>::condition(int maxiters, R tolerance) 1113 { 1114 int dimension = matrix.size(); 1115 int miniters = 3; // minimum number of power method iterations 1116 int i; 1117 int c; 1118 R norm; 1119 R norminv; 1120 R norm1; 1121 R norm2; 1122 1123 // catch corner case of empty matrix 1124 if(dimension <= 0) 1125 return 1.0; 1126 1127 SSVectorBase<R> x(dimension, theLP->tolerances()); 1128 SSVectorBase<R> y(dimension, theLP->tolerances()); 1129 1130 // check whether a regular basis matrix is available 1131 if(status() < REGULAR) 1132 return 0; 1133 1134 if(!matrixIsSetup) 1135 (const_cast<SPxBasisBase<R>*>(this))->loadDesc(thedesc); 1136 1137 if(!factorized) 1138 factorize(); 1139 1140 // initialize vectors 1141 norm1 = 1.0 / (R) dimension; 1142 1143 for(i = 0; i < dimension; i++) 1144 x.add(i, norm1); 1145 1146 y = x; 1147 1148 // compute norm of B 1149 for(c = 0; c < maxiters; ++c) 1150 { 1151 norm2 = norm1; 1152 1153 // y = B*x 1154 multBaseWith(x, y); 1155 norm1 = y.length(); 1156 1157 // stop if converged 1158 if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1) 1159 break; 1160 1161 // x = B^T*y and normalize 1162 multWithBase(y, x); 1163 norm2 = 1.0 / x.length(); 1164 x *= norm2; 1165 } 1166 1167 norm = norm1; 1168 1169 // reinitialize vectors 1170 x.clear(); 1171 y.clear(); 1172 norm1 = 1.0 / (R) dimension; 1173 1174 for(i = 0; i < dimension; i++) 1175 x.add(i, norm1); 1176 1177 y = x; 1178 1179 // compute norm of B^-1 1180 for(c = 0; c < maxiters; ++c) 1181 { 1182 norm2 = norm1; 1183 1184 // x = B^-1*y 1185 factor->solveRight(x, y); 1186 x.setup(); 1187 norm1 = x.length(); 1188 1189 // stop if converged 1190 if(c >= miniters && spxAbs(norm1 - norm2) < tolerance * norm1) 1191 break; 1192 1193 // y = B^-T*x and normalize 1194 factor->solveLeft(y, x); 1195 y.setup(); 1196 norm2 = 1.0 / y.length(); 1197 y *= norm2; 1198 } 1199 1200 norminv = norm1; 1201 1202 return norm * norminv; 1203 } 1204 1205 /* compute one of several matrix metrics based on the diagonal of the LU factorization */ 1206 template <class R> 1207 R SPxBasisBase<R>::getMatrixMetric(int type) 1208 { 1209 R metric = R(infinity); 1210 1211 if(factorized) 1212 metric = factor->matrixMetric(type); 1213 1214 return metric; 1215 } 1216 1217 template <class R> 1218 void SPxBasisBase<R>::dump() 1219 { 1220 assert(status() > NO_PROBLEM); 1221 assert(theLP != 0); 1222 assert(thedesc.nRows() == theLP->nRows()); 1223 assert(thedesc.nCols() == theLP->nCols()); 1224 assert(theLP->dim() == matrix.size()); 1225 1226 int i, basesize; 1227 1228 // Dump regardless of the verbosity level if this method is called. 1229 1230 std::cout << "DBASIS09 Basis entries:"; 1231 basesize = 0; 1232 1233 for(i = 0; i < theLP->nRows(); ++i) 1234 { 1235 if(theLP->isBasic(thedesc.rowStatus(i))) 1236 { 1237 if(basesize % 10 == 0) 1238 std::cout << std::endl << "DBASIS10 "; 1239 1240 SPxRowId id = theLP->SPxLPBase<R>::rId(i); 1241 std::cout << "\tR" << theLP->number(id); 1242 basesize++; 1243 } 1244 } 1245 1246 for(i = 0; i < theLP->nCols(); ++i) 1247 { 1248 if(theLP->isBasic(thedesc.colStatus(i))) 1249 { 1250 if(basesize % 10 == 0) 1251 std::cout << std::endl << "DBASIS11 "; 1252 1253 SPxColId id = theLP->SPxLPBase<R>::cId(i); 1254 std::cout << "\tC" << theLP->number(id); 1255 basesize++; 1256 } 1257 } 1258 1259 std::cout << std::endl; 1260 1261 assert(basesize == matrix.size()); 1262 } 1263 1264 template <class R> 1265 1266 bool SPxBasisBase<R>::isConsistent() const 1267 { 1268 #ifdef ENABLE_CONSISTENCY_CHECKS 1269 int primals = 0; 1270 int i; 1271 1272 if(status() > NO_PROBLEM) 1273 { 1274 if(theLP == 0) 1275 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1276 1277 if(theBaseId.size() != theLP->dim() || matrix.size() != theLP->dim()) 1278 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1279 1280 if(thedesc.nCols() != theLP->nCols() || thedesc.nRows() != theLP->nRows()) 1281 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1282 1283 for(i = 0; i < thedesc.nRows(); ++i) 1284 { 1285 if(thedesc.rowStatus(i) >= 0) 1286 { 1287 if(thedesc.rowStatus(i) != dualRowStatus(i)) 1288 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1289 } 1290 else 1291 ++primals; 1292 } 1293 1294 for(i = 0; i < thedesc.nCols(); ++i) 1295 { 1296 if(thedesc.colStatus(i) >= 0) 1297 { 1298 if(thedesc.colStatus(i) != dualColStatus(i)) 1299 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1300 } 1301 else 1302 ++primals; 1303 } 1304 1305 if(primals != thedesc.nCols()) 1306 return SPX_MSG_INCONSISTENT("SPxBasisBase<R>"); 1307 } 1308 1309 return thedesc.isConsistent() && theBaseId.isConsistent() 1310 && matrix.isConsistent() && factor->isConsistent(); 1311 #else 1312 return true; 1313 #endif // CONSISTENCY_CHECKS 1314 } 1315 1316 template <class R> 1317 SPxBasisBase<R>::SPxBasisBase(Timer::TYPE ttype) 1318 : theLP(0) 1319 , matrixIsSetup(false) 1320 , factor(0) 1321 , factorized(false) 1322 , maxUpdates(200) 1323 , nonzeroFactor(10.0) 1324 , fillFactor(5.0) 1325 , memFactor(1.5) 1326 , iterCount(0) 1327 , lastIterCount(0) 1328 , iterDegenCheck(0) 1329 , updateCount(0) 1330 , totalUpdateCount(0) 1331 , nzCount(1) 1332 , lastMem(0) 1333 , lastFill(0) 1334 , lastNzCount(0) 1335 , theTime(0) 1336 , timerType(ttype) 1337 , lastidx(0) 1338 , minStab(0.0) 1339 , thestatus(NO_PROBLEM) 1340 , freeSlinSolver(false) 1341 , spxout(0) 1342 { 1343 // info: is not consistent at this moment, e.g. because theLP == 0 1344 1345 theTime = TimerFactory::createTimer(timerType); 1346 } 1347 1348 1349 /**@warning Do not change the LP object. 1350 * Only pointer to that object is copied. 1351 * Hint: no problem, we use this function for copy 1352 * constructor of SPxSolverBase<R> 1353 */ 1354 template <class R> 1355 SPxBasisBase<R>::SPxBasisBase(const SPxBasisBase<R>& old) 1356 : theLP(old.theLP) 1357 , theBaseId(old.theBaseId) 1358 , matrix(old.matrix) 1359 , matrixIsSetup(old.matrixIsSetup) 1360 , factor(old.factor) 1361 , factorized(old.factorized) 1362 , maxUpdates(old.maxUpdates) 1363 , nonzeroFactor(old.nonzeroFactor) 1364 , fillFactor(old.fillFactor) 1365 , memFactor(old.memFactor) 1366 , iterCount(old.iterCount) 1367 , lastIterCount(old.lastIterCount) 1368 , iterDegenCheck(old.iterDegenCheck) 1369 , updateCount(old.updateCount) 1370 , totalUpdateCount(old.totalUpdateCount) 1371 , nzCount(old.nzCount) 1372 , lastMem(old.lastMem) 1373 , lastFill(old.lastFill) 1374 , lastNzCount(old.lastNzCount) 1375 , theTime(old.theTime) 1376 , lastin(old.lastin) 1377 , lastout(old.lastout) 1378 , lastidx(old.lastidx) 1379 , minStab(old.minStab) 1380 , thestatus(old.thestatus) 1381 , thedesc(old.thedesc) 1382 , spxout(old.spxout) 1383 { 1384 theTime = TimerFactory::createTimer(old.theTime->type()); 1385 1386 this->factor = old.factor->clone(); 1387 freeSlinSolver = true; 1388 1389 assert(SPxBasisBase<R>::isConsistent()); 1390 } 1391 1392 template <class R> 1393 SPxBasisBase<R>::~SPxBasisBase() 1394 { 1395 1396 assert(!freeSlinSolver || factor != 0); 1397 1398 if(freeSlinSolver) 1399 { 1400 delete factor; 1401 factor = 0; 1402 } 1403 1404 theTime->~Timer(); 1405 spx_free(theTime); 1406 } 1407 1408 template <class R> 1409 1410 /**@warning Note that we do not create a deep copy of the corresponding SPxSolverBase<R> object. 1411 * Only the reference to that object is copied. 1412 */ 1413 SPxBasisBase<R>& SPxBasisBase<R>::operator=(const SPxBasisBase<R>& rhs) 1414 { 1415 1416 assert(!freeSlinSolver || factor != 0); 1417 1418 if(this != &rhs) 1419 { 1420 theLP = rhs.theLP; 1421 theBaseId = rhs.theBaseId; 1422 matrix = rhs.matrix; 1423 matrixIsSetup = rhs.matrixIsSetup; 1424 1425 if(freeSlinSolver) 1426 { 1427 delete factor; 1428 factor = 0; 1429 } 1430 1431 factor = rhs.factor->clone(); 1432 freeSlinSolver = true; 1433 1434 factorized = rhs.factorized; 1435 maxUpdates = rhs.maxUpdates; 1436 nonzeroFactor = rhs.nonzeroFactor; 1437 fillFactor = rhs.fillFactor; 1438 memFactor = rhs.memFactor; 1439 iterCount = rhs.iterCount; 1440 nzCount = rhs.nzCount; 1441 lastFill = rhs.lastFill; 1442 lastNzCount = rhs.lastNzCount; 1443 lastin = rhs.lastin; 1444 lastout = rhs.lastout; 1445 lastidx = rhs.lastidx; 1446 minStab = rhs.minStab; 1447 thestatus = rhs.thestatus; 1448 thedesc = rhs.thedesc; 1449 1450 assert(SPxBasisBase<R>::isConsistent()); 1451 } 1452 1453 return *this; 1454 } 1455 1456 1457 1458 // 1459 // Auxiliary functions. 1460 // 1461 1462 // Pretty-printing of basis status. 1463 template <class R> // Why can't we remove the class R and make it empty? 1464 std::ostream& operator<<(std::ostream& os, 1465 const typename SPxBasisBase<R>::SPxStatus& status) 1466 { 1467 switch(status) 1468 { 1469 case SPxBasisBase<R>::NO_PROBLEM: 1470 os << "NO_PROBLEM"; 1471 break; 1472 1473 case SPxBasisBase<R>::SINGULAR: 1474 os << "SINGULAR"; 1475 break; 1476 1477 case SPxBasisBase<R>::REGULAR: 1478 os << "REGULAR"; 1479 break; 1480 1481 case SPxBasisBase<R>::DUAL: 1482 os << "DUAL"; 1483 break; 1484 1485 case SPxBasisBase<R>::PRIMAL: 1486 os << "PRIMAL"; 1487 break; 1488 1489 case SPxBasisBase<R>::OPTIMAL: 1490 os << "OPTIMAL"; 1491 break; 1492 1493 case SPxBasisBase<R>::UNBOUNDED: 1494 os << "UNBOUNDED"; 1495 break; 1496 1497 case SPxBasisBase<R>::INFEASIBLE: 1498 os << "INFEASIBLE"; 1499 break; 1500 1501 default: 1502 os << "?other?"; 1503 break; 1504 } 1505 1506 return os; 1507 } 1508 1509 1510 } // namespace soplex 1511