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 /**@file slufactor.hpp 26 * @todo SLUfactor seems to be partly an wrapper for CLUFactor (was C). 27 * This should be properly integrated and demangled. 28 * @todo Does is make sense, to call x.clear() when next x.altValues() 29 * is called. 30 */ 31 32 #include <assert.h> 33 #include <sstream> 34 35 #ifdef SOPLEX_DEBUG 36 #include <stdio.h> 37 #endif 38 39 namespace soplex 40 { 41 /// note: we keep this constant since it is just a tradeoff between sparsity and stability and does 42 /// not need to be changed when precisions are decreased 43 #define SOPLEX_MINSTABILITY R(4e-2) 44 45 template <class R> 46 void SLUFactor<R>::solveRight(VectorBase<R>& x, const VectorBase<R>& b) //const 47 { 48 49 this->solveTime->start(); 50 51 this->vec = b; 52 x.clear(); 53 CLUFactor<R>::solveRight(x.get_ptr(), vec.get_ptr()); 54 55 solveCount++; 56 solveTime->stop(); 57 } 58 59 template <class R> 60 void SLUFactor<R>::solveRight(SSVectorBase<R>& x, const SVectorBase<R>& b) //const 61 { 62 63 solveTime->start(); 64 65 vec.assign(b); 66 x.clear(); 67 CLUFactor<R>::solveRight(x.altValues(), vec.get_ptr()); 68 69 solveCount++; 70 solveTime->stop(); 71 } 72 73 template <class R> 74 void SLUFactor<R>::solveRight4update(SSVectorBase<R>& x, const SVectorBase<R>& b) 75 { 76 77 solveTime->start(); 78 79 int m; 80 int n; 81 int f; 82 83 x.clear(); 84 ssvec = b; 85 n = ssvec.size(); 86 R epsilon = this->tolerances()->epsilon(); 87 88 if(this->l.updateType == ETA) 89 { 90 m = this->vSolveRight4update(epsilon, x.altValues(), x.altIndexMem(), 91 ssvec.altValues(), ssvec.altIndexMem(), n, 0, 0, 0); 92 x.setSize(m); 93 //x.forceSetup(); 94 x.unSetup(); 95 eta.setup_and_assign(x); 96 } 97 else 98 { 99 forest.clear(); 100 m = this->vSolveRight4update(epsilon, x.altValues(), x.altIndexMem(), 101 ssvec.altValues(), ssvec.altIndexMem(), n, 102 forest.altValues(), &f, forest.altIndexMem()); 103 forest.setSize(f); 104 forest.forceSetup(); 105 x.setSize(m); 106 x.forceSetup(); 107 } 108 109 usetup = true; 110 ssvec.setSize(0); 111 ssvec.forceSetup(); 112 113 solveCount++; 114 solveTime->stop(); 115 } 116 117 template <class R> 118 void SLUFactor<R>::solve2right4update( 119 SSVectorBase<R>& x, 120 VectorBase<R>& y, 121 const SVectorBase<R>& b, 122 SSVectorBase<R>& rhs) 123 { 124 125 solveTime->start(); 126 127 int m; 128 int n; 129 int f; 130 int* sidx = ssvec.altIndexMem(); 131 ssvec.setSize(0); 132 ssvec.forceSetup(); 133 int rsize = rhs.size(); 134 int* ridx = rhs.altIndexMem(); 135 R epsilon = this->tolerances()->epsilon(); 136 137 x.clear(); 138 y.clear(); 139 usetup = true; 140 ssvec = b; 141 142 if(this->l.updateType == ETA) 143 { 144 n = ssvec.size(); 145 m = this->vSolveRight4update2(epsilon, x.altValues(), x.altIndexMem(), 146 ssvec.get_ptr(), sidx, n, y.get_ptr(), 147 epsilon, rhs.altValues(), ridx, rsize, 0, 0, 0); 148 x.setSize(m); 149 // x.forceSetup(); 150 x.unSetup(); 151 eta.setup_and_assign(x); 152 } 153 else 154 { 155 forest.clear(); 156 n = ssvec.size(); 157 m = this->vSolveRight4update2(epsilon, x.altValues(), x.altIndexMem(), 158 ssvec.get_ptr(), sidx, n, y.get_ptr(), 159 epsilon, rhs.altValues(), ridx, rsize, 160 forest.altValues(), &f, forest.altIndexMem()); 161 x.setSize(m); 162 x.forceSetup(); 163 forest.setSize(f); 164 forest.forceSetup(); 165 } 166 167 rhs.forceSetup(); 168 ssvec.setSize(0); 169 ssvec.forceSetup(); 170 171 solveCount += 2; 172 solveTime->stop(); 173 } 174 175 template <class R> 176 void SLUFactor<R>::solve2right4update( 177 SSVectorBase<R>& x, 178 SSVectorBase<R>& y, 179 const SVectorBase<R>& b, 180 SSVectorBase<R>& rhs) 181 { 182 183 solveTime->start(); 184 185 int n; 186 int f; 187 int* sidx = ssvec.altIndexMem(); 188 ssvec.setSize(0); 189 ssvec.forceSetup(); 190 int rsize = rhs.size(); 191 int* ridx = rhs.altIndexMem(); 192 R epsilon = this->tolerances()->epsilon(); 193 194 x.clear(); 195 y.clear(); 196 usetup = true; 197 ssvec = b; 198 199 if(this->l.updateType == ETA) 200 { 201 n = ssvec.size(); 202 this->vSolveRight4update2sparse(epsilon, x.altValues(), x.altIndexMem(), 203 ssvec.get_ptr(), sidx, n, 204 epsilon, y.altValues(), y.altIndexMem(), 205 rhs.altValues(), ridx, rsize, 206 0, 0, 0); 207 x.setSize(n); 208 // x.forceSetup(); 209 x.unSetup(); 210 y.setSize(rsize); 211 y.unSetup(); 212 eta.setup_and_assign(x); 213 } 214 else 215 { 216 forest.clear(); 217 n = ssvec.size(); 218 this->vSolveRight4update2sparse(epsilon, x.altValues(), x.altIndexMem(), 219 ssvec.get_ptr(), sidx, n, 220 epsilon, y.altValues(), y.altIndexMem(), 221 rhs.altValues(), ridx, rsize, 222 forest.altValues(), &f, forest.altIndexMem()); 223 x.setSize(n); 224 x.forceSetup(); 225 y.setSize(rsize); 226 y.forceSetup(); 227 forest.setSize(f); 228 forest.forceSetup(); 229 } 230 231 rhs.forceSetup(); 232 ssvec.setSize(0); 233 ssvec.forceSetup(); 234 235 solveCount += 2; 236 solveTime->stop(); 237 } 238 239 240 template <class R> 241 void SLUFactor<R>::solve3right4update( 242 SSVectorBase<R>& x, 243 VectorBase<R>& y, 244 VectorBase<R>& y2, 245 const SVectorBase<R>& b, 246 SSVectorBase<R>& rhs, 247 SSVectorBase<R>& rhs2) 248 { 249 250 solveTime->start(); 251 252 int m; 253 int n; 254 int f; 255 int* sidx = ssvec.altIndexMem(); 256 ssvec.setSize(0); 257 ssvec.forceSetup(); 258 int rsize = rhs.size(); 259 int* ridx = rhs.altIndexMem(); 260 int rsize2 = rhs2.size(); 261 int* ridx2 = rhs2.altIndexMem(); 262 R epsilon = this->tolerances()->epsilon(); 263 264 x.clear(); 265 y.clear(); 266 y2.clear(); 267 usetup = true; 268 ssvec = b; 269 270 if(this->l.updateType == ETA) 271 { 272 n = ssvec.size(); 273 m = this->vSolveRight4update3(epsilon, 274 x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n, 275 y.get_ptr(), epsilon, rhs.altValues(), ridx, rsize, 276 y2.get_ptr(), epsilon, rhs2.altValues(), ridx2, rsize2, 277 0, 0, 0); 278 x.setSize(m); 279 // x.forceSetup(); 280 x.unSetup(); 281 eta.setup_and_assign(x); 282 } 283 else 284 { 285 forest.clear(); 286 n = ssvec.size(); 287 m = this->vSolveRight4update3(epsilon, 288 x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n, 289 y.get_ptr(), epsilon, rhs.altValues(), ridx, rsize, 290 y2.get_ptr(), epsilon, rhs2.altValues(), ridx2, rsize2, 291 forest.altValues(), &f, forest.altIndexMem()); 292 x.setSize(m); 293 x.forceSetup(); 294 forest.setSize(f); 295 forest.forceSetup(); 296 } 297 298 rhs.forceSetup(); 299 rhs2.forceSetup(); 300 ssvec.setSize(0); 301 ssvec.forceSetup(); 302 303 solveCount += 3; 304 solveTime->stop(); 305 } 306 307 template <class R> 308 void SLUFactor<R>::solve3right4update( 309 SSVectorBase<R>& x, 310 SSVectorBase<R>& y, 311 SSVectorBase<R>& y2, 312 const SVectorBase<R>& b, 313 SSVectorBase<R>& rhs, 314 SSVectorBase<R>& rhs2) 315 { 316 solveTime->start(); 317 318 int n; 319 int f; 320 int* sidx = ssvec.altIndexMem(); 321 ssvec.setSize(0); 322 ssvec.forceSetup(); 323 int rsize = rhs.size(); 324 int* ridx = rhs.altIndexMem(); 325 int rsize2 = rhs2.size(); 326 int* ridx2 = rhs2.altIndexMem(); 327 R epsilon = this->tolerances()->epsilon(); 328 329 x.clear(); 330 y.clear(); 331 y2.clear(); 332 usetup = true; 333 ssvec = b; 334 335 if(this->l.updateType == ETA) 336 { 337 n = ssvec.size(); 338 this->vSolveRight4update3sparse(epsilon, x.altValues(), x.altIndexMem(), 339 ssvec.get_ptr(), sidx, n, 340 epsilon, y.altValues(), y.altIndexMem(), 341 rhs.altValues(), ridx, rsize, 342 epsilon, y2.altValues(), y2.altIndexMem(), 343 rhs2.altValues(), ridx2, rsize2, 344 0, 0, 0); 345 x.setSize(n); 346 // x.forceSetup(); 347 x.unSetup(); 348 y.setSize(rsize); 349 y.unSetup(); 350 y2.setSize(rsize2); 351 y2.unSetup(); 352 eta.setup_and_assign(x); 353 } 354 else 355 { 356 forest.clear(); 357 n = ssvec.size(); 358 this->vSolveRight4update3sparse(epsilon, x.altValues(), x.altIndexMem(), 359 ssvec.get_ptr(), sidx, n, 360 epsilon, y.altValues(), y.altIndexMem(), 361 rhs.altValues(), ridx, rsize, 362 epsilon, y2.altValues(), y2.altIndexMem(), 363 rhs2.altValues(), ridx2, rsize2, 364 forest.altValues(), &f, forest.altIndexMem()); 365 x.setSize(n); 366 x.forceSetup(); 367 y.setSize(rsize); 368 y.forceSetup(); 369 y2.setSize(rsize2); 370 y2.forceSetup(); 371 372 forest.setSize(f); 373 forest.forceSetup(); 374 } 375 376 rhs.forceSetup(); 377 rhs2.forceSetup(); 378 ssvec.setSize(0); 379 ssvec.forceSetup(); 380 381 solveCount += 3; 382 solveTime->stop(); 383 } 384 385 386 template <class R> 387 void SLUFactor<R>::solveLeft(VectorBase<R>& x, const VectorBase<R>& b) //const 388 { 389 390 solveTime->start(); 391 392 vec = b; 393 x.clear(); 394 CLUFactor<R>::solveLeft(x.get_ptr(), vec.get_ptr()); 395 396 solveCount++; 397 solveTime->stop(); 398 } 399 400 template <class R> 401 void SLUFactor<R>::solveLeft(SSVectorBase<R>& x, const SVectorBase<R>& b) //const 402 { 403 R epsilon = this->tolerances()->epsilon(); 404 405 solveTime->start(); 406 407 // copy to SSVec is done to avoid having to deal with the Nonzero datatype 408 // TODO change SVec to standard sparse format 409 ssvec.assign(b); 410 411 x.clear(); 412 int sz = ssvec.size(); // see .altValues() 413 int n = this->vSolveLeft(epsilon, x.altValues(), x.altIndexMem(), 414 ssvec.altValues(), ssvec.altIndexMem(), sz); 415 416 if(n > 0) 417 { 418 x.setSize(n); 419 x.forceSetup(); 420 } 421 else 422 x.unSetup(); 423 424 ssvec.setSize(0); 425 ssvec.forceSetup(); 426 427 solveCount++; 428 solveTime->stop(); 429 } 430 431 template <class R> 432 void SLUFactor<R>::solveLeft( 433 SSVectorBase<R>& x, 434 VectorBase<R>& y, 435 const SVectorBase<R>& rhs1, 436 SSVectorBase<R>& rhs2) //const 437 { 438 439 solveTime->start(); 440 441 int n; 442 R* svec = ssvec.altValues(); 443 int* sidx = ssvec.altIndexMem(); 444 int rn = rhs2.size(); 445 int* ridx = rhs2.altIndexMem(); 446 R epsilon = this->tolerances()->epsilon(); 447 448 x.clear(); 449 y.clear(); 450 ssvec.assign(rhs1); 451 n = ssvec.size(); // see altValues(); 452 n = this->vSolveLeft2(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n, 453 y.get_ptr(), rhs2.altValues(), ridx, rn); 454 455 // this will unsetup x 456 x.setSize(n); 457 458 if(n > 0) 459 x.forceSetup(); 460 461 ssvec.setSize(0); 462 ssvec.forceSetup(); 463 464 solveCount += 2; 465 solveTime->stop(); 466 } 467 468 template <class R> 469 void SLUFactor<R>::solveLeft( 470 SSVectorBase<R>& x, 471 SSVectorBase<R>& y, 472 const SVectorBase<R>& rhs1, 473 SSVectorBase<R>& rhs2) //const 474 { 475 476 solveTime->start(); 477 478 int n1, n2; 479 R* svec = ssvec.altValues(); 480 int* sidx = ssvec.altIndexMem(); 481 R epsilon = this->tolerances()->epsilon(); 482 483 x.clear(); 484 y.clear(); 485 ssvec.assign(rhs1); 486 n1 = ssvec.size(); // see altValues(); 487 n2 = rhs2.size(); 488 489 if(n2 < 10) 490 { 491 this->vSolveLeft2sparse(epsilon, 492 x.altValues(), x.altIndexMem(), 493 svec, sidx, n1, 494 y.altValues(), y.altIndexMem(), 495 rhs2.altValues(), rhs2.altIndexMem(), n2); 496 y.setSize(n2); 497 498 if(n2 > 0) 499 y.forceSetup(); 500 } 501 else 502 { 503 n1 = this->vSolveLeft2(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n1, 504 y.altValues(), rhs2.altValues(), rhs2.altIndexMem(), n2); 505 // y.setup(); 506 } 507 508 x.setSize(n1); 509 510 if(n1 > 0) 511 x.forceSetup(); 512 513 // if (n2 > 0) 514 // y.forceSetup(); 515 516 ssvec.setSize(0); 517 ssvec.forceSetup(); 518 519 solveCount += 2; 520 solveTime->stop(); 521 } 522 523 524 template <class R> 525 void SLUFactor<R>::solveLeft( 526 SSVectorBase<R>& x, 527 VectorBase<R>& y, 528 VectorBase<R>& z, 529 const SVectorBase<R>& rhs1, 530 SSVectorBase<R>& rhs2, 531 SSVectorBase<R>& rhs3) 532 { 533 534 solveTime->start(); 535 536 int n, n2, n3; 537 R* svec = ssvec.altValues(); 538 int* sidx = ssvec.altIndexMem(); 539 R epsilon = this->tolerances()->epsilon(); 540 541 x.clear(); 542 y.clear(); 543 z.clear(); 544 ssvec.assign(rhs1); 545 n = ssvec.size(); // see altValues(); 546 n2 = rhs2.size(); 547 n3 = rhs3.size(); 548 549 n = this->vSolveLeft3(epsilon, x.altValues(), x.altIndexMem(), svec, sidx, n, 550 y.get_ptr(), rhs2.altValues(), rhs2.altIndexMem(), n2, 551 z.get_ptr(), rhs3.altValues(), rhs3.altIndexMem(), n3); 552 553 x.setSize(n); 554 555 if(n > 0) 556 x.forceSetup(); 557 558 ssvec.setSize(0); 559 ssvec.forceSetup(); 560 561 solveCount += 3; 562 solveTime->stop(); 563 } 564 565 template <class R> 566 void SLUFactor<R>::solveLeft( 567 SSVectorBase<R>& x, 568 SSVectorBase<R>& y, 569 SSVectorBase<R>& z, 570 const SVectorBase<R>& rhs1, 571 SSVectorBase<R>& rhs2, 572 SSVectorBase<R>& rhs3) 573 { 574 575 solveTime->start(); 576 577 int n1, n2, n3; 578 R* svec = ssvec.altValues(); 579 int* sidx = ssvec.altIndexMem(); 580 R epsilon = this->tolerances()->epsilon(); 581 582 x.clear(); 583 y.clear(); 584 z.clear(); 585 ssvec.assign(rhs1); 586 n1 = ssvec.size(); // see altValues(); 587 n2 = rhs2.size(); 588 n3 = rhs3.size(); 589 this->vSolveLeft3sparse(epsilon, 590 x.altValues(), x.altIndexMem(), 591 svec, sidx, n1, 592 y.altValues(), y.altIndexMem(), 593 rhs2.altValues(), rhs2.altIndexMem(), n2, 594 z.altValues(), z.altIndexMem(), 595 rhs3.altValues(), rhs3.altIndexMem(), n3); 596 x.setSize(n1); 597 y.setSize(n2); 598 z.setSize(n3); 599 600 if(n1 > 0) 601 x.forceSetup(); 602 603 if(n2 > 0) 604 y.forceSetup(); 605 606 if(n3 > 0) 607 z.forceSetup(); 608 609 ssvec.setSize(0); 610 ssvec.forceSetup(); 611 612 solveCount += 3; 613 solveTime->stop(); 614 } 615 616 617 template <class R> 618 R SLUFactor<R>::stability() const 619 { 620 if(status() != this->OK) 621 return 0; 622 623 if(this->maxabs < this->initMaxabs) 624 return 1; 625 626 assert(this->maxabs != 0.0); 627 return this->initMaxabs / this->maxabs; 628 } 629 630 template <class R> 631 R SLUFactor<R>::matrixMetric(int type) const 632 { 633 R result = 0.0; 634 635 // catch corner case of empty matrix 636 if(dim() == 0) 637 return 1.0; 638 639 switch(type) 640 { 641 // compute condition estimate by ratio of max/min of elements on the diagonal 642 case 0: 643 { 644 R mindiag = spxAbs(this->diag[0]); 645 R maxdiag = spxAbs(this->diag[0]); 646 647 for(int i = 1; i < dim(); ++i) 648 { 649 R absdiag = spxAbs(this->diag[i]); 650 651 if(absdiag < mindiag) 652 mindiag = absdiag; 653 else if(absdiag > maxdiag) 654 maxdiag = absdiag; 655 } 656 657 result = maxdiag / mindiag; 658 break; 659 } 660 661 // compute sum of inverses of all elements on the diagonal 662 case 1: 663 result = 0.0; 664 665 for(int i = 0; i < dim(); ++i) 666 result += 1.0 / this->diag[i]; 667 668 break; 669 670 // compute determinant (product of all diagonal elements of U) 671 case 2: 672 result = 1.0; 673 674 for(int i = 0; i < dim(); ++i) 675 result *= this->diag[i]; 676 677 result = 1.0 / result; 678 break; 679 } 680 681 return result; 682 } 683 684 template <class R> 685 std::string SLUFactor<R>::statistics() const 686 { 687 std::stringstream s; 688 s << "Factorizations : " << std::setw(10) << getFactorCount() << std::endl 689 << " Time spent : " << std::setw(10) << std::fixed << std::setprecision( 690 2) << getFactorTime() << std::endl 691 << "Solves : " << std::setw(10) << getSolveCount() << std::endl 692 << " Time spent : " << std::setw(10) << getSolveTime() << std::endl; 693 694 return s.str(); 695 } 696 697 template <class R> 698 void SLUFactor<R>::changeEta(int idx, SSVectorBase<R>& et) 699 { 700 701 int es = et.size(); // see altValues() 702 this->update(idx, et.altValues(), et.altIndexMem(), es); 703 et.setSize(0); 704 et.forceSetup(); 705 } 706 707 template <class R> 708 typename SLUFactor<R>::Status SLUFactor<R>::change( 709 int idx, 710 const SVectorBase<R>& subst, 711 const SSVectorBase<R>* e) 712 { 713 714 // BH 2005-08-23: The boolean usetup indicates that an "update VectorBase<R>" 715 // has been set up. I suppose that SSVectorBase<R> forest is this 716 // update VectorBase<R>, which is set up by solveRight4update() and 717 // solve2right4update() in order to optimize the basis update. 718 719 if(usetup) 720 { 721 if(this->l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates 722 { 723 // BH 2005-08-19: The size of a SSVectorBase<R> is the size of the 724 // index set, i.e. the number of nonzeros which is only 725 // defined if the SSVectorBase<R> is set up. Since 726 // SSVectorBase<R> ::altValues() calls unSetup() the size needs to be 727 // stored before the following call. 728 int fsize = forest.size(); // see altValues() 729 this->forestUpdate(idx, forest.altValues(), fsize, forest.altIndexMem()); 730 forest.setSize(0); 731 forest.forceSetup(); 732 } 733 else 734 { 735 assert(this->l.updateType == ETA); 736 changeEta(idx, eta); 737 } 738 } 739 else if(e != 0) // ETA updates 740 { 741 this->l.updateType = ETA; 742 this->updateNoClear(idx, e->values(), e->indexMem(), e->size()); 743 this->l.updateType = uptype; 744 } 745 else if(this->l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates 746 { 747 assert(0); // probably this part is never called. 748 // forestUpdate() with the last parameter set to NULL should fail. 749 forest = subst; 750 CLUFactor<R>::solveLright(forest.altValues()); 751 this->forestUpdate(idx, forest.altValues(), 0, nullptr); 752 forest.setSize(0); 753 forest.forceSetup(); 754 } 755 else // ETA updates 756 { 757 assert(this->l.updateType == ETA); 758 vec = subst; 759 eta.clear(); 760 CLUFactor<R>::solveRight(eta.altValues(), vec.get_ptr()); 761 changeEta(idx, eta); 762 } 763 764 usetup = false; 765 766 SPxOut::debug(this, "DSLUFA01\tupdated\t\tstability = {}\n", stability()); 767 768 return status(); 769 } 770 771 template <class R> 772 void SLUFactor<R>::clear() 773 { 774 775 this->rowMemMult = 5; /* factor of minimum Memory * #of nonzeros */ 776 this->colMemMult = 5; /* factor of minimum Memory * #of nonzeros */ 777 this->lMemMult = 1; /* factor of minimum Memory * #of nonzeros */ 778 779 this->l.firstUpdate = 0; 780 this->l.firstUnused = 0; 781 this->thedim = 0; 782 783 usetup = false; 784 this->maxabs = 1; 785 this->initMaxabs = 1; 786 lastThreshold = minThreshold; 787 minStability = SOPLEX_MINSTABILITY; 788 this->stat = this->UNLOADED; 789 790 vec.clear(); 791 eta.clear(); 792 ssvec.clear(); 793 forest.clear(); 794 795 this->u.row.size = 100; 796 this->u.col.size = 100; 797 this->l.size = 100; 798 this->l.startSize = 100; 799 800 if(this->l.ridx) 801 spx_free(this->l.ridx); 802 803 if(this->l.rbeg) 804 spx_free(this->l.rbeg); 805 806 if(this->l.rorig) 807 spx_free(this->l.rorig); 808 809 if(this->l.rperm) 810 spx_free(this->l.rperm); 811 812 if(!this->u.row.val.empty()) 813 this->u.row.val.clear(); 814 815 if(this->u.row.idx) 816 spx_free(this->u.row.idx); 817 818 if(this->u.col.idx) 819 spx_free(this->u.col.idx); 820 821 if(this->l.val.empty()) 822 this->l.val.clear(); 823 824 if(this->l.idx) 825 spx_free(this->l.idx); 826 827 if(this->l.start) 828 spx_free(this->l.start); 829 830 if(this->l.row) 831 spx_free(this->l.row); 832 833 // G clear() is used in constructor of SLUFactor<R> so we have to 834 // G clean up if anything goes wrong here 835 try 836 { 837 this->u.row.val.resize(this->u.row.size); 838 spx_alloc(this->u.row.idx, this->u.row.size); 839 spx_alloc(this->u.col.idx, this->u.col.size); 840 841 this->l.val.resize(this->l.size); 842 spx_alloc(this->l.idx, this->l.size); 843 spx_alloc(this->l.start, this->l.startSize); 844 spx_alloc(this->l.row, this->l.startSize); 845 } 846 catch(const SPxMemoryException& x) 847 { 848 freeAll(); 849 throw x; 850 } 851 } 852 853 /** assignment used to implement operator=() and copy constructor. 854 * If this is initialised, freeAll() has to be called before. 855 * Class objects from SLUFactor<R> are not copied here. 856 */ 857 template <class R> 858 void SLUFactor<R>::assign(const SLUFactor<R>& old) 859 { 860 this->spxout = old.spxout; 861 862 solveTime = TimerFactory::createTimer(old.solveTime->type()); 863 this->factorTime = TimerFactory::createTimer(old.factorTime->type()); 864 865 // slufactor 866 uptype = old.uptype; 867 minThreshold = old.minThreshold; 868 minStability = old.minStability; 869 lastThreshold = old.lastThreshold; 870 871 // clufactor 872 this->stat = old.stat; 873 this->thedim = old.thedim; 874 this->nzCnt = old.nzCnt; 875 this->initMaxabs = old.initMaxabs; 876 this->maxabs = old.maxabs; 877 this->rowMemMult = old.rowMemMult; 878 this->colMemMult = old.colMemMult; 879 this->lMemMult = old.lMemMult; 880 881 spx_alloc(this->row.perm, this->thedim); 882 spx_alloc(this->row.orig, this->thedim); 883 spx_alloc(this->col.perm, this->thedim); 884 spx_alloc(this->col.orig, this->thedim); 885 // spx_alloc(this->diag, this->thedim); 886 this->diag.reserve(this->thedim); // small performance improvement before copying 887 888 memcpy(this->row.perm, old.row.perm, (unsigned int)this->thedim * sizeof(*this->row.perm)); 889 memcpy(this->row.orig, old.row.orig, (unsigned int)this->thedim * sizeof(*this->row.orig)); 890 memcpy(this->col.perm, old.col.perm, (unsigned int)this->thedim * sizeof(*this->col.perm)); 891 memcpy(this->col.orig, old.col.orig, (unsigned int)this->thedim * sizeof(*this->col.orig)); 892 893 this->diag = old.diag; 894 895 this->work = vec.get_ptr(); 896 897 /* setup U 898 */ 899 this->u.row.size = old.u.row.size; 900 this->u.row.used = old.u.row.used; 901 902 spx_alloc(this->u.row.elem, this->thedim); 903 this->u.row.val.reserve(this->u.row.size); // small performance improvement 904 spx_alloc(this->u.row.idx, this->u.row.size); 905 spx_alloc(this->u.row.start, this->thedim + 1); 906 spx_alloc(this->u.row.len, this->thedim + 1); 907 spx_alloc(this->u.row.max, this->thedim + 1); 908 909 memcpy(this->u.row.elem, old.u.row.elem, 910 (unsigned int)this->thedim * sizeof(*this->u.row.elem)); 911 this->u.row.val = old.u.row.val; 912 memcpy(this->u.row.idx, old.u.row.idx, 913 (unsigned int)this->u.row.size * sizeof(*this->u.row.idx)); 914 memcpy(this->u.row.start, old.u.row.start, 915 (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.start)); 916 memcpy(this->u.row.len, old.u.row.len, 917 (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.len)); 918 memcpy(this->u.row.max, old.u.row.max, 919 (unsigned int)(this->thedim + 1) * sizeof(*this->u.row.max)); 920 921 // need to make row list ok ? 922 if(this->thedim > 0 && this->stat == this->OK) 923 { 924 this->u.row.list.idx = old.u.row.list.idx; // .idx neu 925 926 const typename CLUFactor<R>::Dring* oring = &old.u.row.list; 927 typename CLUFactor<R>::Dring* ring = &this->u.row.list; 928 929 while(oring->next != &old.u.row.list) 930 { 931 ring->next = &this->u.row.elem[oring->next->idx]; 932 ring->next->prev = ring; 933 oring = oring->next; 934 ring = ring->next; 935 } 936 937 ring->next = &this->u.row.list; 938 ring->next->prev = ring; 939 } 940 941 this->u.col.size = old.u.col.size; 942 this->u.col.used = old.u.col.used; 943 944 spx_alloc(this->u.col.elem, this->thedim); 945 spx_alloc(this->u.col.idx, this->u.col.size); 946 spx_alloc(this->u.col.start, this->thedim + 1); 947 spx_alloc(this->u.col.len, this->thedim + 1); 948 spx_alloc(this->u.col.max, this->thedim + 1); 949 950 if(!old.u.col.val.empty()) 951 { 952 this->u.col.val.reserve(this->u.col.size); // small performance improvement before copying 953 this->u.col.val = old.u.col.val; 954 } 955 else 956 { 957 this->u.col.val.clear(); 958 } 959 960 memcpy(this->u.col.elem, old.u.col.elem, 961 (unsigned int)this->thedim * sizeof(*this->u.col.elem)); 962 memcpy(this->u.col.idx, old.u.col.idx, 963 (unsigned int)this->u.col.size * sizeof(*this->u.col.idx)); 964 memcpy(this->u.col.start, old.u.col.start, 965 (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.start)); 966 memcpy(this->u.col.len, old.u.col.len, 967 (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.len)); 968 memcpy(this->u.col.max, old.u.col.max, 969 (unsigned int)(this->thedim + 1) * sizeof(*this->u.col.max)); 970 971 // need to make col list ok 972 if(this->thedim > 0 && this->stat == this->OK) 973 { 974 this->u.col.list.idx = old.u.col.list.idx; // .idx neu 975 976 const typename CLUFactor<R>::Dring* oring = &old.u.col.list; 977 typename CLUFactor<R>::Dring* ring = &this->u.col.list; 978 979 while(oring->next != &old.u.col.list) 980 { 981 ring->next = &this->u.col.elem[oring->next->idx]; 982 ring->next->prev = ring; 983 oring = oring->next; 984 ring = ring->next; 985 } 986 987 ring->next = &this->u.col.list; 988 ring->next->prev = ring; 989 } 990 991 /* Setup L 992 */ 993 this->l.size = old.l.size; 994 this->l.startSize = old.l.startSize; 995 this->l.firstUpdate = old.l.firstUpdate; 996 this->l.firstUnused = old.l.firstUnused; 997 this->l.updateType = old.l.updateType; 998 999 this->l.val.reserve(this->l.size); // small performance improvement for copying 1000 spx_alloc(this->l.idx, this->l.size); 1001 spx_alloc(this->l.start, this->l.startSize); 1002 spx_alloc(this->l.row, this->l.startSize); 1003 1004 this->l.val = old.l.val; 1005 memcpy(this->l.idx, old.l.idx, (unsigned int)this->l.size * sizeof(*this->l.idx)); 1006 memcpy(this->l.start, old.l.start, (unsigned int)this->l.startSize * sizeof(*this->l.start)); 1007 memcpy(this->l.row, old.l.row, (unsigned int)this->l.startSize * sizeof(*this->l.row)); 1008 1009 if(!this->l.rval.empty()) 1010 { 1011 assert(old.l.ridx != 0); 1012 assert(old.l.rbeg != 0); 1013 assert(old.l.rorig != 0); 1014 assert(old.l.rperm != 0); 1015 1016 int memsize = this->l.start[this->l.firstUpdate]; 1017 1018 this->l.rval.reserve(memsize); // small performance improvement for copying 1019 spx_alloc(this->l.ridx, memsize); 1020 spx_alloc(this->l.rbeg, this->thedim + 1); 1021 spx_alloc(this->l.rorig, this->thedim); 1022 spx_alloc(this->l.rperm, this->thedim); 1023 1024 this->l.rval = old.l.rval; 1025 memcpy(this->l.ridx, old.l.ridx, (unsigned int)memsize * sizeof(*this->l.ridx)); 1026 memcpy(this->l.rbeg, old.l.rbeg, (unsigned int)(this->thedim + 1) * sizeof(*this->l.rbeg)); 1027 memcpy(this->l.rorig, old.l.rorig, (unsigned int)this->thedim * sizeof(*this->l.rorig)); 1028 memcpy(this->l.rperm, old.l.rperm, (unsigned int)this->thedim * sizeof(*this->l.rperm)); 1029 } 1030 else 1031 { 1032 assert(old.l.ridx == 0); 1033 assert(old.l.rbeg == 0); 1034 assert(old.l.rorig == 0); 1035 assert(old.l.rperm == 0); 1036 1037 this->l.ridx = 0; 1038 this->l.rbeg = 0; 1039 this->l.rorig = 0; 1040 this->l.rperm = 0; 1041 } 1042 1043 assert(this->row.perm != 0); 1044 assert(this->row.orig != 0); 1045 assert(this->col.perm != 0); 1046 assert(this->col.orig != 0); 1047 1048 assert(this->u.row.elem != 0); 1049 assert(this->u.row.idx != 0); 1050 assert(this->u.row.start != 0); 1051 assert(this->u.row.len != 0); 1052 assert(this->u.row.max != 0); 1053 1054 assert(this->u.col.elem != 0); 1055 assert(this->u.col.idx != 0); 1056 assert(this->u.col.start != 0); 1057 assert(this->u.col.len != 0); 1058 assert(this->u.col.max != 0); 1059 1060 assert(this->l.idx != 0); 1061 assert(this->l.start != 0); 1062 assert(this->l.row != 0); 1063 1064 } 1065 1066 template <class R> 1067 SLUFactor<R>& SLUFactor<R>::operator=(const SLUFactor<R>& old) 1068 { 1069 1070 if(this != &old) 1071 { 1072 // we don't need to copy them, because they are temporary vectors 1073 vec.clear(); 1074 ssvec.clear(); 1075 1076 eta = old.eta; 1077 forest = old.forest; 1078 1079 timerType = old.timerType; 1080 1081 freeAll(); 1082 1083 try 1084 { 1085 assign(old); 1086 } 1087 catch(const SPxMemoryException& x) 1088 { 1089 freeAll(); 1090 throw x; 1091 } 1092 1093 assert(isConsistent()); 1094 } 1095 1096 return *this; 1097 } 1098 1099 template <class R> 1100 SLUFactor<R>::SLUFactor() 1101 : vec(1) 1102 , ssvec(1) 1103 , usetup(false) 1104 , uptype(FOREST_TOMLIN) 1105 , eta(1) 1106 , forest(1) 1107 , minThreshold(0.01) 1108 , timerType(Timer::USER_TIME) 1109 { 1110 this->row.perm = 0; 1111 this->row.orig = 0; 1112 this->col.perm = 0; 1113 this->col.orig = 0; 1114 this->u.row.elem = 0; 1115 this->u.row.idx = 0; 1116 this->u.row.start = 0; 1117 this->u.row.len = 0; 1118 this->u.row.max = 0; 1119 this->u.col.elem = 0; 1120 this->u.col.idx = 0; 1121 this->u.col.start = 0; 1122 this->u.col.len = 0; 1123 this->u.col.max = 0; 1124 this->l.idx = 0; 1125 this->l.start = 0; 1126 this->l.row = 0; 1127 this->l.ridx = 0; 1128 this->l.rbeg = 0; 1129 this->l.rorig = 0; 1130 this->l.rperm = 0; 1131 1132 this->nzCnt = 0; 1133 this->thedim = 0; 1134 1135 try 1136 { 1137 solveTime = TimerFactory::createTimer(timerType); 1138 this->factorTime = TimerFactory::createTimer(timerType); 1139 spx_alloc(this->row.perm, this->thedim); 1140 spx_alloc(this->row.orig, this->thedim); 1141 spx_alloc(this->col.perm, this->thedim); 1142 spx_alloc(this->col.orig, this->thedim); 1143 1144 this->diag.resize(this->thedim); 1145 1146 this->work = vec.get_ptr(); 1147 1148 this->u.row.size = 1; 1149 this->u.row.used = 0; 1150 spx_alloc(this->u.row.elem, this->thedim); 1151 this->u.row.val.resize(this->u.row.size); 1152 spx_alloc(this->u.row.idx, this->u.row.size); 1153 spx_alloc(this->u.row.start, this->thedim + 1); 1154 spx_alloc(this->u.row.len, this->thedim + 1); 1155 spx_alloc(this->u.row.max, this->thedim + 1); 1156 1157 this->u.row.list.idx = this->thedim; 1158 this->u.row.start[this->thedim] = 0; 1159 this->u.row.max [this->thedim] = 0; 1160 this->u.row.len [this->thedim] = 0; 1161 1162 this->u.col.size = 1; 1163 this->u.col.used = 0; 1164 spx_alloc(this->u.col.elem, this->thedim); 1165 spx_alloc(this->u.col.idx, this->u.col.size); 1166 spx_alloc(this->u.col.start, this->thedim + 1); 1167 spx_alloc(this->u.col.len, this->thedim + 1); 1168 spx_alloc(this->u.col.max, this->thedim + 1); 1169 1170 this->u.col.list.idx = this->thedim; 1171 this->u.col.start[this->thedim] = 0; 1172 this->u.col.max[this->thedim] = 0; 1173 this->u.col.len[this->thedim] = 0; 1174 1175 this->l.size = 1; 1176 1177 this->l.val.resize(this->l.size); 1178 spx_alloc(this->l.idx, this->l.size); 1179 1180 this->l.startSize = 1; 1181 this->l.firstUpdate = 0; 1182 this->l.firstUnused = 0; 1183 1184 spx_alloc(this->l.start, this->l.startSize); 1185 spx_alloc(this->l.row, this->l.startSize); 1186 } 1187 catch(const SPxMemoryException& x) 1188 { 1189 freeAll(); 1190 throw x; 1191 } 1192 1193 this->l.ridx = 0; 1194 this->l.rbeg = 0; 1195 this->l.rorig = 0; 1196 this->l.rperm = 0; 1197 1198 SLUFactor<R>::clear(); // clear() is virtual 1199 1200 this->factorCount = 0; 1201 this->hugeValues = 0; 1202 solveCount = 0; 1203 assert(this->row.perm != 0); 1204 1205 assert(this->row.orig != 0); 1206 assert(this->col.perm != 0); 1207 assert(this->col.orig != 0); 1208 1209 assert(this->u.row.elem != 0); 1210 assert(this->u.row.idx != 0); 1211 assert(this->u.row.start != 0); 1212 assert(this->u.row.len != 0); 1213 assert(this->u.row.max != 0); 1214 1215 assert(this->u.col.elem != 0); 1216 assert(this->u.col.idx != 0); 1217 assert(this->u.col.start != 0); 1218 assert(this->u.col.len != 0); 1219 assert(this->u.col.max != 0); 1220 1221 assert(this->l.idx != 0); 1222 assert(this->l.start != 0); 1223 assert(this->l.row != 0); 1224 1225 assert(SLUFactor<R>::isConsistent()); 1226 } 1227 1228 template <class R> 1229 SLUFactor<R>::SLUFactor(const SLUFactor<R>& old) 1230 : SLinSolver<R>(old) 1231 , vec(1) // we don't need to copy it, because they are temporary vectors 1232 , ssvec(1) // we don't need to copy it, because they are temporary vectors 1233 , usetup(old.usetup) 1234 , eta(old.eta) 1235 , forest(old.forest) 1236 , timerType(old.timerType) 1237 { 1238 this->row.perm = 0; 1239 this->row.orig = 0; 1240 this->col.perm = 0; 1241 this->col.orig = 0; 1242 this->u.row.elem = 0; 1243 this->u.row.val.clear(); 1244 this->u.row.idx = 0; 1245 this->u.row.start = 0; 1246 this->u.row.len = 0; 1247 this->u.row.max = 0; 1248 this->u.col.elem = 0; 1249 this->u.col.idx = 0; 1250 this->u.col.start = 0; 1251 this->u.col.len = 0; 1252 this->u.col.max = 0; 1253 this->l.idx = 0; 1254 this->l.start = 0; 1255 this->l.row = 0; 1256 this->l.ridx = 0; 1257 this->l.rbeg = 0; 1258 this->l.rorig = 0; 1259 this->l.rperm = 0; 1260 1261 solveCount = 0; 1262 1263 try 1264 { 1265 assign(old); 1266 } 1267 catch(const SPxMemoryException& x) 1268 { 1269 freeAll(); 1270 throw x; 1271 } 1272 1273 assert(SLUFactor<R>::isConsistent()); 1274 } 1275 1276 template <class R> 1277 void SLUFactor<R>::freeAll() 1278 { 1279 1280 if(this->row.perm) 1281 spx_free(this->row.perm); 1282 1283 if(this->row.orig) 1284 spx_free(this->row.orig); 1285 1286 if(this->col.perm) 1287 spx_free(this->col.perm); 1288 1289 if(this->col.orig) 1290 spx_free(this->col.orig); 1291 1292 if(this->u.row.elem) 1293 spx_free(this->u.row.elem); 1294 1295 if(!this->u.row.val.empty()) 1296 this->u.row.val.clear(); 1297 1298 if(this->u.row.idx) 1299 spx_free(this->u.row.idx); 1300 1301 if(this->u.row.start) 1302 spx_free(this->u.row.start); 1303 1304 if(this->u.row.len) 1305 spx_free(this->u.row.len); 1306 1307 if(this->u.row.max) 1308 spx_free(this->u.row.max); 1309 1310 if(this->u.col.elem) 1311 spx_free(this->u.col.elem); 1312 1313 if(this->u.col.idx) 1314 spx_free(this->u.col.idx); 1315 1316 if(this->u.col.start) 1317 spx_free(this->u.col.start); 1318 1319 if(this->u.col.len) 1320 spx_free(this->u.col.len); 1321 1322 if(this->u.col.max) 1323 spx_free(this->u.col.max); 1324 1325 if(!this->l.val.empty()) 1326 this->l.val.clear(); 1327 1328 if(this->l.idx) 1329 spx_free(this->l.idx); 1330 1331 if(this->l.start) 1332 spx_free(this->l.start); 1333 1334 if(this->l.row) 1335 spx_free(this->l.row); 1336 1337 1338 if(!this->u.col.val.empty()) 1339 this->u.col.val.clear(); 1340 1341 if(this->l.ridx) 1342 spx_free(this->l.ridx); 1343 1344 if(this->l.rbeg) 1345 spx_free(this->l.rbeg); 1346 1347 if(this->l.rorig) 1348 spx_free(this->l.rorig); 1349 1350 if(this->l.rperm) 1351 spx_free(this->l.rperm); 1352 1353 if(solveTime) 1354 { 1355 solveTime->~Timer(); 1356 spx_free(solveTime); 1357 } 1358 1359 if(this->factorTime) 1360 { 1361 this->factorTime->~Timer(); 1362 spx_free(this->factorTime); 1363 } 1364 } 1365 1366 template <class R> 1367 SLUFactor<R>::~SLUFactor() 1368 { 1369 freeAll(); 1370 } 1371 1372 template <class R> 1373 static R betterThreshold(R th, Real epsilon) 1374 { 1375 assert(th < R(1.0)); 1376 1377 if(LT(th, R(0.1), epsilon)) 1378 th *= R(10.0); 1379 else if(LT(th, R(0.9), epsilon)) 1380 th = (th + R(1.0)) / R(2.0); 1381 else if(LT(th, R(0.999), epsilon)) 1382 th = R(0.99999); 1383 1384 assert(th < R(1.0)); 1385 1386 return th; 1387 } 1388 1389 template <class R> 1390 typename SLUFactor<R>::Status SLUFactor<R>::load(const SVectorBase<R>* matrix[], int dm) 1391 { 1392 assert(dm >= 0); 1393 assert(matrix != 0); 1394 1395 R lastStability = stability(); 1396 1397 initDR(this->u.row.list); 1398 initDR(this->u.col.list); 1399 1400 usetup = false; 1401 this->l.updateType = uptype; 1402 this->l.firstUpdate = 0; 1403 this->l.firstUnused = 0; 1404 1405 if(dm != this->thedim) 1406 { 1407 clear(); 1408 1409 this->thedim = dm; 1410 vec.reDim(this->thedim); 1411 ssvec.reDim(this->thedim); 1412 eta.reDim(this->thedim); 1413 forest.reDim(this->thedim); 1414 this->work = vec.get_ptr(); 1415 1416 spx_realloc(this->row.perm, this->thedim); 1417 spx_realloc(this->row.orig, this->thedim); 1418 spx_realloc(this->col.perm, this->thedim); 1419 spx_realloc(this->col.orig, this->thedim); 1420 this->diag.resize(this->thedim); 1421 1422 spx_realloc(this->u.row.elem, this->thedim); 1423 spx_realloc(this->u.row.len, this->thedim + 1); 1424 spx_realloc(this->u.row.max, this->thedim + 1); 1425 spx_realloc(this->u.row.start, this->thedim + 1); 1426 1427 spx_realloc(this->u.col.elem, this->thedim); 1428 spx_realloc(this->u.col.len, this->thedim + 1); 1429 spx_realloc(this->u.col.max, this->thedim + 1); 1430 spx_realloc(this->u.col.start, this->thedim + 1); 1431 1432 this->l.startSize = this->thedim + SOPLEX_MAXUPDATES; 1433 1434 spx_realloc(this->l.row, this->l.startSize); 1435 spx_realloc(this->l.start, this->l.startSize); 1436 } 1437 // the last factorization was reasonably stable, so we decrease the Markowitz threshold (stored in lastThreshold) in 1438 // order to favour sparsity 1439 else if(lastStability > 2.0 * minStability) 1440 { 1441 // we reset lastThreshold to its previous value in the sequence minThreshold, betterThreshold(minThreshold), 1442 // betterThreshold(betterThreshold(minThreshold)), ... 1443 R last = minThreshold; 1444 R better = betterThreshold(last, this->tolerances()->epsilon()); 1445 1446 while(better < lastThreshold) 1447 { 1448 last = better; 1449 better = betterThreshold(last, this->tolerances()->epsilon()); 1450 } 1451 1452 lastThreshold = last; 1453 1454 // we reset the minimum stability (which might have been decreased below) to ensure that the increased sparsity 1455 // does not hurt the stability 1456 minStability = 2 * SOPLEX_MINSTABILITY; 1457 } 1458 1459 this->u.row.list.idx = this->thedim; 1460 this->u.row.start[this->thedim] = 0; 1461 this->u.row.max[this->thedim] = 0; 1462 this->u.row.len[this->thedim] = 0; 1463 1464 this->u.col.list.idx = this->thedim; 1465 this->u.col.start[this->thedim] = 0; 1466 this->u.col.max[this->thedim] = 0; 1467 this->u.col.len[this->thedim] = 0; 1468 1469 for(;;) 1470 { 1471 ///@todo if the factorization fails with stat = SINGULAR, distinuish between proven singularity (e.g., because of 1472 ///an empty column) and singularity due to numerics, that could be avoided by changing minStability and 1473 ///lastThreshold; in the first case, we want to abort, otherwise change the numerics 1474 this->stat = this->OK; 1475 this->factor(matrix, lastThreshold, this->tolerances()->epsilonFactorization()); 1476 1477 // finish if the factorization is stable 1478 if(stability() >= minStability) 1479 break; 1480 1481 // otherwise, we increase the Markowitz threshold 1482 R x = lastThreshold; 1483 lastThreshold = betterThreshold(lastThreshold, this->tolerances()->epsilon()); 1484 1485 // until it doesn't change anymore at its maximum value 1486 if(EQ(x, lastThreshold, this->tolerances()->epsilon())) 1487 break; 1488 1489 // we relax the stability requirement 1490 minStability /= 2.0; 1491 1492 SPX_MSG_INFO3((*this->spxout), (*this->spxout) << 1493 "ISLUFA01 refactorizing with increased Markowitz threshold: " 1494 << lastThreshold << std::endl;) 1495 } 1496 1497 SPxOut::debug(this, "DSLUFA02 threshold = {} \tstability = {}\tminStability = {}\n", lastThreshold, 1498 stability(), minStability); 1499 SPX_DEBUG( 1500 { 1501 int i; 1502 FILE* fl = fopen("dump.lp", "w"); 1503 std::cout << "DSLUFA03 Basis:\n"; 1504 int j = 0; 1505 1506 for(i = 0; i < dim(); ++i) 1507 j += matrix[i]->size(); 1508 for(i = 0; i < dim(); ++i) 1509 { 1510 for(j = 0; j < matrix[i]->size(); ++j) 1511 { 1512 fprintf(fl, "%8d %8d ", 1513 i + 1, matrix[i]->index(j) + 1); 1514 std::cout << matrix[i]->value(j) << std::endl; 1515 } 1516 } 1517 fclose(fl); 1518 std::cout << "DSLUFA04 LU-Factors:" << std::endl; 1519 dump(); 1520 1521 std::cout << "DSLUFA05 threshold = " << lastThreshold 1522 << "\tstability = " << stability() << std::endl; 1523 } 1524 ) 1525 1526 assert(isConsistent()); 1527 return Status(this->stat); 1528 } 1529 1530 1531 template <class R> 1532 bool SLUFactor<R>::isConsistent() const 1533 { 1534 #ifdef ENABLE_CONSISTENCY_CHECKS 1535 return CLUFactor<R>::isConsistent(); 1536 #else 1537 return true; 1538 #endif 1539 } 1540 1541 template <class R> 1542 void SLUFactor<R>::dump() const 1543 { 1544 CLUFactor<R>::dump(); 1545 } 1546 } // namespace soplex 1547