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 26 #include <assert.h> 27 #include <iostream> 28 29 #include "soplex/spxdefines.h" 30 #include "soplex/svset.h" 31 #include "soplex/sorter.h" 32 33 namespace soplex 34 { 35 #define SOPLEX_STABLE 1e-3 // the sparsest row/column may only have a pivot of size STABLE*maxEntry 36 37 38 template <class R> 39 bool SPxWeightST<R>::isConsistent() const 40 { 41 #ifdef ENABLE_CONSISTENCY_CHECKS 42 return rowWeight.isConsistent() 43 && colWeight.isConsistent() 44 && rowRight.isConsistent() 45 && colUp.isConsistent(); 46 // && SPxStarter::isConsistent(); // not yet implemented 47 #else 48 return true; 49 #endif 50 } 51 52 /* Generating the Starting Basis 53 The generation of a starting basis works as follows: After setting up the 54 preference arrays #weight# and #coWeight#, #Id#s are selected to be dual in 55 a greedy manner. Initially, the first #Id# is taken. Then the next #Id# is 56 checked wheter its vector is linearly dependend of the vectors of the #Id#s 57 allready selected. If not, it is added to them. This is iterated until a 58 full matrix has been constructed. 59 60 Testing for linear independence is done very much along the lines of LU 61 factorization. A vector is taken, and updated with all previous L-vectors. 62 Then the maximal absolut element is selected as pivot element for computing 63 the L-vector. This is stored for further processing. 64 */ 65 66 /* 67 The following two functions set the status of |id| to primal or dual, 68 respectively. 69 */ 70 template <class R> 71 void SPxWeightST<R>::setPrimalStatus( 72 typename SPxBasisBase<R>::Desc& desc, 73 const SPxSolverBase<R>& base, 74 const SPxId& id) 75 { 76 if(id.isSPxRowId()) 77 { 78 int n = base.number(SPxRowId(id)); 79 80 if(base.rhs(n) >= R(infinity)) 81 { 82 if(base.lhs(n) <= R(-infinity)) 83 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_FREE; 84 else 85 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER; 86 } 87 else 88 { 89 if(base.lhs(n) <= R(-infinity)) 90 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER; 91 else if(base.lhs(n) >= base.rhs(n) - base.epsilon()) 92 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_FIXED; 93 else if(rowRight[n]) 94 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER; 95 else 96 desc.rowStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER; 97 } 98 } 99 else 100 { 101 int n = base.number(SPxColId(id)); 102 103 if(base.SPxLPBase<R>::upper(n) >= R(infinity)) 104 { 105 if(base.SPxLPBase<R>::lower(n) <= R(-infinity)) 106 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_FREE; 107 else 108 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER; 109 } 110 else 111 { 112 if(base.SPxLPBase<R>::lower(n) <= R(-infinity)) 113 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER; 114 else if(base.SPxLPBase<R>::lower(n) >= base.SPxLPBase<R>::upper(n) - base.epsilon()) 115 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_FIXED; 116 else if(colUp[n]) 117 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_UPPER; 118 else 119 desc.colStatus(n) = SPxBasisBase<R>::Desc::P_ON_LOWER; 120 } 121 } 122 } 123 124 // ---------------------------------------------------------------- 125 template <class R> 126 static void setDualStatus( 127 typename SPxBasisBase<R>::Desc& desc, 128 const SPxSolverBase<R>& base, 129 const SPxId& id) 130 { 131 if(id.isSPxRowId()) 132 { 133 int n = base.number(SPxRowId(id)); 134 desc.rowStatus(n) = base.basis().dualRowStatus(n); 135 } 136 else 137 { 138 int n = base.number(SPxColId(id)); 139 desc.colStatus(n) = base.basis().dualColStatus(n); 140 } 141 } 142 // ---------------------------------------------------------------- 143 144 /// Compare class for row weights, used for sorting. 145 template <typename T> 146 struct Compare 147 { 148 public: 149 /// constructor 150 Compare() : weight(0) {} 151 // const SPxSolverBase* base; ///< the solver 152 const T* weight; ///< the weights to compare 153 154 /// compares the weights 155 T operator()(int i1, int i2) const 156 { 157 return weight[i1] - weight[i2]; 158 } 159 }; 160 161 // ---------------------------------------------------------------- 162 /** 163 The following method initializes \p pref such that it contains the 164 set of \ref soplex::SPxId "SPxIds" ordered following \p rowWeight and 165 \p colWeight. For the sorting we take the following approach: first 166 we sort the rows, then the columns. Finally we perform a mergesort 167 of both. 168 */ 169 template <class R> 170 static void initPrefs( 171 DataArray < SPxId >& pref, 172 const SPxSolverBase<R>& base, 173 const Array<R>& rowWeight, 174 const Array<R>& colWeight) 175 { 176 DataArray<int> row(base.nRows()); 177 DataArray<int> col(base.nCols()); 178 int i; 179 int j; 180 int k; 181 182 Compare<R> compare; 183 // compare.base = &base; 184 185 for(i = 0; i < base.nRows(); ++i) 186 row[i] = i; 187 188 compare.weight = rowWeight.get_const_ptr(); 189 190 SPxQuicksort(row.get_ptr(), row.size(), compare); // Sort rows 191 192 for(i = 0; i < base.nCols(); ++i) 193 col[i] = i; 194 195 compare.weight = colWeight.get_const_ptr(); 196 197 SPxQuicksort(col.get_ptr(), col.size(), compare); // Sort column 198 199 i = 0; 200 j = 0; 201 k = 0; 202 203 while(k < pref.size()) // merge sort 204 { 205 if(rowWeight[row[i]] < colWeight[col[j]]) 206 { 207 pref[k++] = base.rId(row[i++]); 208 209 if(i >= base.nRows()) 210 while(k < pref.size()) 211 pref[k++] = base.cId(col[j++]); 212 } 213 else 214 { 215 pref[k++] = base.cId(col[j++]); 216 217 if(j >= base.nCols()) 218 while(k < pref.size()) 219 pref[k++] = base.rId(row[i++]); 220 } 221 } 222 223 assert(i == base.nRows()); 224 assert(j == base.nCols()); 225 } 226 227 // ---------------------------------------------------------------- 228 template <class R> 229 void SPxWeightST<R>::generate(SPxSolverBase<R>& base) 230 { 231 SPxId tmpId; 232 233 forbidden.reSize(base.dim()); 234 rowWeight.reSize(base.nRows()); 235 colWeight.reSize(base.nCols()); 236 rowRight.reSize(base.nRows()); 237 colUp.reSize(base.nCols()); 238 239 if(base.rep() == SPxSolverBase<R>::COLUMN) 240 { 241 weight = &colWeight; 242 coWeight = &rowWeight; 243 } 244 else 245 { 246 weight = &rowWeight; 247 coWeight = &colWeight; 248 } 249 250 assert(weight->size() == base.coDim()); 251 assert(coWeight->size() == base.dim()); 252 253 setupWeights(base); 254 255 typename SPxBasisBase<R>::Desc desc(base); 256 // desc.reSize(base.nRows(), base.nCols()); 257 258 DataArray < SPxId > pref(base.nRows() + base.nCols()); 259 initPrefs(pref, base, rowWeight, colWeight); 260 261 int i; 262 int stepi; 263 int j; 264 int sel; 265 266 for(i = 0; i < base.dim(); ++i) 267 forbidden[i] = 0; 268 269 if(base.rep() == SPxSolverBase<R>::COLUMN) 270 { 271 // in COLUMN rep we scan from beginning to end 272 i = 0; 273 stepi = 1; 274 } 275 else 276 { 277 // in ROW rep we scan from end to beginning 278 i = pref.size() - 1; 279 stepi = -1; 280 } 281 282 int dim = base.dim(); 283 R maxEntry = 0; 284 285 for(; i >= 0 && i < pref.size(); i += stepi) 286 { 287 tmpId = pref[i]; 288 const SVectorBase<R>& vec = base.vector(tmpId); 289 sel = -1; 290 291 // column or row singleton ? 292 if(vec.size() == 1) 293 { 294 int idx = vec.index(0); 295 296 if(forbidden[idx] < 2) 297 { 298 sel = idx; 299 dim += (forbidden[idx] > 0) ? 1 : 0; 300 } 301 } 302 else 303 { 304 maxEntry = vec.maxAbs(); 305 306 // initialize the nonzero counter 307 int minRowEntries = base.nRows(); 308 309 // find a stable index with a sparse row/column 310 for(j = vec.size(); --j >= 0;) 311 { 312 R x = vec.value(j); 313 int k = vec.index(j); 314 int nRowEntries = base.coVector(k).size(); 315 316 if(!forbidden[k] 317 && (spxAbs(x) > this->tolerances()->scaleAccordingToEpsilon(SOPLEX_STABLE) * maxEntry) 318 && (nRowEntries < minRowEntries)) 319 { 320 minRowEntries = nRowEntries; 321 sel = k; 322 } 323 } 324 } 325 326 // we found a valid index 327 if(sel >= 0) 328 { 329 SPX_DEBUG( 330 331 if(pref[i].type() == SPxId::ROW_ID) 332 std::cout << "DWEIST01 r" << base.number(pref[i]); 333 else 334 std::cout << "DWEIST02 c" << base.number(pref[i]); 335 ) 336 337 forbidden[sel] = 2; 338 339 // put current column/row into basis 340 if(base.rep() == SPxSolverBase<R>::COLUMN) 341 setDualStatus(desc, base, pref[i]); 342 else 343 setPrimalStatus(desc, base, pref[i]); 344 345 for(j = vec.size(); --j >= 0;) 346 { 347 R x = vec.value(j); 348 int k = vec.index(j); 349 R feastol = this->tolerances()->floatingPointFeastol(); 350 351 if(!forbidden[k] && (x > feastol * maxEntry || -x > feastol * maxEntry)) 352 { 353 forbidden[k] = 1; 354 --dim; 355 } 356 } 357 358 if(--dim == 0) 359 { 360 //@ for(++i; i < pref.size(); ++i) 361 if(base.rep() == SPxSolverBase<R>::COLUMN) 362 { 363 // set all remaining indeces to nonbasic status 364 for(i += stepi; i >= 0 && i < pref.size(); i += stepi) 365 setPrimalStatus(desc, base, pref[i]); 366 367 // fill up the basis wherever linear independence is assured 368 for(i = forbidden.size(); --i >= 0;) 369 { 370 if(forbidden[i] < 2) 371 setDualStatus(desc, base, base.coId(i)); 372 } 373 } 374 else 375 { 376 for(i += stepi; i >= 0 && i < pref.size(); i += stepi) 377 setDualStatus(desc, base, pref[i]); 378 379 for(i = forbidden.size(); --i >= 0;) 380 { 381 if(forbidden[i] < 2) 382 setPrimalStatus(desc, base, base.coId(i)); 383 } 384 } 385 386 break; 387 } 388 } 389 // sel == -1 390 else if(base.rep() == SPxSolverBase<R>::COLUMN) 391 setPrimalStatus(desc, base, pref[i]); 392 else 393 setDualStatus(desc, base, pref[i]); 394 395 #ifndef NDEBUG 396 { 397 int n, m; 398 399 for(n = 0, m = forbidden.size(); n < forbidden.size(); ++n) 400 m -= (forbidden[n] != 0) ? 1 : 0; 401 402 assert(m == dim); 403 } 404 #endif // NDEBUG 405 } 406 407 assert(dim == 0); 408 409 base.loadBasis(desc); 410 #ifdef TEST 411 base.init(); 412 413 int changed = 0; 414 const VectorBase<R>& pvec = base.pVec(); 415 416 for(i = pvec.dim() - 1; i >= 0; --i) 417 { 418 if(desc.colStatus(i) == SPxBasisBase<R>::Desc::P_ON_UPPER 419 && base.lower(i) > R(-infinity) && pvec[i] > base.maxObj(i)) 420 { 421 changed = 1; 422 desc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_LOWER; 423 } 424 else if(desc.colStatus(i) == SPxBasisBase<R>::Desc::P_ON_LOWER 425 && base.upper(i) < R(infinity) && pvec[i] < base.maxObj(i)) 426 { 427 changed = 1; 428 desc.colStatus(i) = SPxBasisBase<R>::Desc::P_ON_UPPER; 429 } 430 } 431 432 if(changed) 433 { 434 std::cout << "changed basis\n"; 435 base.loadBasis(desc); 436 } 437 else 438 std::cout << "nothing changed\n"; 439 440 #endif // TEST 441 } 442 443 // ---------------------------------------------------------------- 444 445 /* Computation of Weights 446 */ 447 template <class R> 448 void SPxWeightST<R>::setupWeights(SPxSolverBase<R>& base) 449 { 450 const VectorBase<R>& obj = base.maxObj(); 451 const VectorBase<R>& low = base.lower(); 452 const VectorBase<R>& up = base.upper(); 453 const VectorBase<R>& rhs = base.rhs(); 454 const VectorBase<R>& lhs = base.lhs(); 455 int i; 456 457 R eps = base.epsilon(); 458 R maxabs = 1.0; 459 460 // find absolut biggest entry in bounds and left-/right hand side 461 for(i = 0; i < base.nCols(); i++) 462 { 463 if((up[i] < R(infinity)) && (spxAbs(up[i]) > maxabs)) 464 maxabs = spxAbs(up[i]); 465 466 if((low[i] > R(-infinity)) && (spxAbs(low[i]) > maxabs)) 467 maxabs = spxAbs(low[i]); 468 } 469 470 for(i = 0; i < base.nRows(); i++) 471 { 472 if((rhs[i] < R(infinity)) && (spxAbs(rhs[i]) > maxabs)) 473 maxabs = spxAbs(rhs[i]); 474 475 if((lhs[i] > R(-infinity)) && (spxAbs(lhs[i]) > maxabs)) 476 maxabs = spxAbs(lhs[i]); 477 } 478 479 /**@todo The comments are wrong. The first is for dual simplex and 480 * the secound for primal one. Is anything else wrong? 481 * Also the values are nearly the same for both cases. 482 * Should this be ? Changed the values for 483 * r_fixed to 0 because of maros-r7. It is not clear why 484 * this makes a difference because all constraints in that 485 * instance are of equality type. 486 * Why is rowRight sometimes not set? 487 */ 488 if(base.rep() * base.type() > 0) // primal simplex 489 { 490 const R ax = 1e-3 / obj.maxAbs(); 491 const R bx = 1.0 / maxabs; 492 const R nne = ax / base.nRows(); // 1e-4 * ax; 493 const R c_fixed = 1e+5; 494 const R r_fixed = 0; // TK20010103: was 1e+4 (maros-r7) 495 const R c_dbl_bounded = 1e+1; 496 const R r_dbl_bounded = 0; 497 const R c_bounded = 1e+1; 498 const R r_bounded = 0; 499 const R c_free = -1e+4; 500 const R r_free = -1e+5; 501 502 for(i = base.nCols() - 1; i >= 0; i--) 503 { 504 R n = nne * (base.colVector(i).size() - 1); // very small value that is zero for col singletons 505 R x = ax * obj[i]; // this is at most 1e-3, probably a lot smaller 506 R u = bx * up [i]; // this is at most 1, probably a lot smaller 507 R l = bx * low[i]; // this is at most 1, probably a lot smaller 508 509 if(up[i] < R(infinity)) 510 { 511 if(spxAbs(low[i] - up[i]) < eps) 512 colWeight[i] = c_fixed + n + spxAbs(x); 513 else if(low[i] > R(-infinity)) 514 { 515 colWeight[i] = c_dbl_bounded + l - u + n; 516 517 l = spxAbs(l); 518 u = spxAbs(u); 519 520 if(u < l) 521 { 522 colUp[i] = true; 523 colWeight[i] += x; 524 } 525 else 526 { 527 colUp[i] = false; 528 colWeight[i] -= x; 529 } 530 } 531 else 532 { 533 colWeight[i] = c_bounded - u + x + n; 534 colUp[i] = true; 535 } 536 } 537 else 538 { 539 if(low[i] > R(-infinity)) 540 { 541 colWeight[i] = c_bounded + l + n - x; 542 colUp[i] = false; 543 } 544 else 545 { 546 colWeight[i] = c_free + n - spxAbs(x); 547 } 548 } 549 } 550 551 for(i = base.nRows() - 1; i >= 0; i--) 552 { 553 if(rhs[i] < R(infinity)) 554 { 555 if(spxAbs(lhs[i] - rhs[i]) < eps) 556 { 557 rowWeight[i] = r_fixed; 558 } 559 else if(lhs[i] > R(-infinity)) 560 { 561 R u = bx * rhs[i]; 562 R l = bx * lhs[i]; 563 564 rowWeight[i] = r_dbl_bounded + l - u; 565 rowRight[i] = spxAbs(u) < spxAbs(l); 566 } 567 else 568 { 569 rowWeight[i] = r_bounded - bx * rhs[i]; 570 rowRight[i] = true; 571 } 572 } 573 else 574 { 575 if(lhs[i] > R(-infinity)) 576 { 577 rowWeight[i] = r_bounded + bx * lhs[i]; 578 rowRight[i] = false; 579 } 580 else 581 { 582 rowWeight[i] = r_free; 583 } 584 } 585 } 586 } 587 else 588 { 589 assert(base.rep() * base.type() < 0); // dual simplex 590 591 const R ax = 1.0 / obj.maxAbs(); 592 const R bx = 1e-2 / maxabs; 593 const R nne = 1e-4 * bx; 594 const R c_fixed = 1e+5; 595 const R r_fixed = 1e+4; 596 const R c_dbl_bounded = 1; 597 const R r_dbl_bounded = 0; 598 const R c_bounded = 0; 599 const R r_bounded = 0; 600 const R c_free = -1e+4; 601 const R r_free = -1e+5; 602 603 for(i = base.nCols() - 1; i >= 0; i--) 604 { 605 R n = nne * (base.colVector(i).size() - 1); 606 R x = ax * obj[i]; 607 R u = bx * up [i]; 608 R l = bx * low[i]; 609 610 if(up[i] < R(infinity)) 611 { 612 if(spxAbs(low[i] - up[i]) < eps) 613 colWeight[i] = c_fixed + n + spxAbs(x); 614 else if(low[i] > R(-infinity)) 615 { 616 if(x > 0) 617 { 618 colWeight[i] = c_dbl_bounded + x - u + n; 619 colUp[i] = true; 620 } 621 else 622 { 623 colWeight[i] = c_dbl_bounded - x + l + n; 624 colUp[i] = false; 625 } 626 } 627 else 628 { 629 colWeight[i] = c_bounded - u + x + n; 630 colUp[i] = true; 631 } 632 } 633 else 634 { 635 if(low[i] > R(-infinity)) 636 { 637 colWeight[i] = c_bounded - x + l + n; 638 colUp[i] = false; 639 } 640 else 641 colWeight[i] = c_free + n - spxAbs(x); 642 } 643 } 644 645 for(i = base.nRows() - 1; i >= 0; i--) 646 { 647 const R len1 = 1; // (base.rowVector(i).length() + base.epsilon()); 648 R n = 0; // nne * (base.rowVector(i).size() - 1); 649 R u = bx * len1 * rhs[i]; 650 R l = bx * len1 * lhs[i]; 651 R x = ax * len1 * (obj * base.rowVector(i)); 652 653 if(rhs[i] < R(infinity)) 654 { 655 if(spxAbs(lhs[i] - rhs[i]) < eps) 656 rowWeight[i] = r_fixed + n + spxAbs(x); 657 else if(lhs[i] > R(-infinity)) 658 { 659 if(x > 0) 660 { 661 rowWeight[i] = r_dbl_bounded + x - u + n; 662 rowRight[i] = true; 663 } 664 else 665 { 666 rowWeight[i] = r_dbl_bounded - x + l + n; 667 rowRight[i] = false; 668 } 669 } 670 else 671 { 672 rowWeight[i] = r_bounded - u + n + x; 673 rowRight[i] = true; 674 } 675 } 676 else 677 { 678 if(lhs[i] > R(-infinity)) 679 { 680 rowWeight[i] = r_bounded + l + n - x; 681 rowRight[i] = false; 682 } 683 else 684 { 685 rowWeight[i] = r_free + n - spxAbs(x); 686 } 687 } 688 } 689 } 690 691 SPX_DEBUG( 692 { 693 for(i = 0; i < base.nCols(); i++) 694 std::cout << "C i= " << i 695 << " up= " << colUp[i] 696 << " w= " << colWeight[i] 697 << std::endl; 698 for(i = 0; i < base.nRows(); i++) 699 std::cout << "R i= " << i 700 << " rr= " << rowRight[i] 701 << " w= " << rowWeight[i] 702 << std::endl; 703 }) 704 } 705 } // namespace soplex 706