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_rational.hpp 26 * @todo SLUFactorRational seems to be partly an wrapper for CLUFactorRational (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 #include "soplex/cring.h" 35 #include "soplex/spxdefines.h" 36 37 #ifdef SOPLEX_DEBUG 38 #include <stdio.h> 39 #endif 40 41 namespace soplex 42 { 43 #define SOPLEX_MINSTABILITYRAT SOPLEX_REAL(4e-2) 44 45 inline void SLUFactorRational::solveRight(VectorRational& x, const VectorRational& b) //const 46 { 47 48 solveTime->start(); 49 50 vec = b; 51 CLUFactorRational::solveRight(x.get_ptr(), vec.get_ptr()); 52 53 solveCount++; 54 solveTime->stop(); 55 } 56 57 inline void SLUFactorRational::solveRight(SSVectorRational& x, const SVectorRational& b) //const 58 { 59 60 solveTime->start(); 61 62 vec.assign(b); 63 x.clear(); 64 CLUFactorRational::solveRight(x.altValues(), vec.get_ptr()); 65 66 solveCount++; 67 solveTime->stop(); 68 } 69 70 inline void SLUFactorRational::solveRight4update(SSVectorRational& x, const SVectorRational& b) 71 { 72 73 solveTime->start(); 74 75 int m; 76 int n; 77 int f = 0; 78 79 x.clear(); 80 ssvec = b; 81 n = ssvec.size(); 82 83 if(l.updateType == ETA) 84 { 85 m = vSolveRight4update(x.altValues(), x.altIndexMem(), 86 ssvec.altValues(), ssvec.altIndexMem(), n, 0, 0, 0); 87 x.setSize(m); 88 //x.forceSetup(); 89 x.unSetup(); 90 eta.setup_and_assign(x); 91 } 92 else 93 { 94 forest.clear(); 95 m = vSolveRight4update(x.altValues(), x.altIndexMem(), 96 ssvec.altValues(), ssvec.altIndexMem(), n, 97 forest.altValues(), &f, forest.altIndexMem()); 98 forest.setSize(f); 99 forest.forceSetup(); 100 x.setSize(m); 101 x.forceSetup(); 102 } 103 104 usetup = true; 105 106 solveCount++; 107 solveTime->stop(); 108 } 109 110 inline void SLUFactorRational::solve2right4update( 111 SSVectorRational& x, 112 VectorRational& y, 113 const SVectorRational& b, 114 SSVectorRational& rhs) 115 { 116 117 solveTime->start(); 118 119 int m; 120 int n; 121 int f = 0; 122 int* sidx = ssvec.altIndexMem(); 123 int rsize = rhs.size(); 124 int* ridx = rhs.altIndexMem(); 125 126 x.clear(); 127 y.clear(); 128 usetup = true; 129 ssvec = b; 130 131 if(l.updateType == ETA) 132 { 133 n = ssvec.size(); 134 m = vSolveRight4update2(x.altValues(), x.altIndexMem(), 135 ssvec.get_ptr(), sidx, n, y.get_ptr(), 136 rhs.altValues(), ridx, rsize, 0, 0, 0); 137 x.setSize(m); 138 // x.forceSetup(); 139 x.unSetup(); 140 eta.setup_and_assign(x); 141 } 142 else 143 { 144 forest.clear(); 145 n = ssvec.size(); 146 m = vSolveRight4update2(x.altValues(), x.altIndexMem(), 147 ssvec.get_ptr(), sidx, n, y.get_ptr(), 148 rhs.altValues(), ridx, rsize, 149 forest.altValues(), &f, forest.altIndexMem()); 150 x.setSize(m); 151 x.forceSetup(); 152 forest.setSize(f); 153 forest.forceSetup(); 154 } 155 156 solveCount++; 157 solveTime->stop(); 158 } 159 160 inline void SLUFactorRational::solve3right4update( 161 SSVectorRational& x, 162 VectorRational& y, 163 VectorRational& y2, 164 const SVectorRational& b, 165 SSVectorRational& rhs, 166 SSVectorRational& rhs2) 167 { 168 169 solveTime->start(); 170 171 int m; 172 int n; 173 int f; 174 int* sidx = ssvec.altIndexMem(); 175 int rsize = rhs.size(); 176 int* ridx = rhs.altIndexMem(); 177 int rsize2 = rhs2.size(); 178 int* ridx2 = rhs2.altIndexMem(); 179 180 x.clear(); 181 y.clear(); 182 y2.clear(); 183 usetup = true; 184 ssvec = b; 185 186 if(l.updateType == ETA) 187 { 188 n = ssvec.size(); 189 m = vSolveRight4update3(x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n, 190 y.get_ptr(), rhs.altValues(), ridx, rsize, 191 y2.get_ptr(), rhs2.altValues(), ridx2, rsize2, 192 0, 0, 0); 193 x.setSize(m); 194 // x.forceSetup(); 195 x.unSetup(); 196 eta.setup_and_assign(x); 197 } 198 else 199 { 200 forest.clear(); 201 n = ssvec.size(); 202 m = vSolveRight4update3(x.altValues(), x.altIndexMem(), ssvec.get_ptr(), sidx, n, 203 y.get_ptr(), rhs.altValues(), ridx, rsize, 204 y2.get_ptr(), rhs2.altValues(), ridx2, rsize2, 205 forest.altValues(), &f, forest.altIndexMem()); 206 x.setSize(m); 207 x.forceSetup(); 208 forest.setSize(f); 209 forest.forceSetup(); 210 } 211 212 solveCount++; 213 solveTime->stop(); 214 } 215 216 inline void SLUFactorRational::solveLeft(VectorRational& x, const VectorRational& b) //const 217 { 218 219 solveTime->start(); 220 221 vec = b; 222 ///@todo Why is x.clear() here used and not with solveRight() ? 223 x.clear(); 224 CLUFactorRational::solveLeft(x.get_ptr(), vec.get_ptr()); 225 226 solveCount++; 227 solveTime->stop(); 228 } 229 230 inline void SLUFactorRational::solveLeft(SSVectorRational& x, const SVectorRational& b) //const 231 { 232 233 solveTime->start(); 234 235 ssvec.assign(b); 236 237 x.clear(); 238 int sz = ssvec.size(); // see .altValues() 239 int n = vSolveLeft(x.altValues(), x.altIndexMem(), 240 ssvec.altValues(), ssvec.altIndexMem(), sz); 241 242 if(n > 0) 243 { 244 x.setSize(n); 245 x.forceSetup(); 246 } 247 else 248 x.unSetup(); 249 250 ssvec.setSize(0); 251 ssvec.forceSetup(); 252 253 solveCount++; 254 solveTime->stop(); 255 } 256 257 inline void SLUFactorRational::solveLeft( 258 SSVectorRational& x, 259 VectorRational& y, 260 const SVectorRational& rhs1, 261 SSVectorRational& rhs2) //const 262 { 263 264 solveTime->start(); 265 266 int n; 267 Rational* svec = ssvec.altValues(); 268 int* sidx = ssvec.altIndexMem(); 269 int rn = rhs2.size(); 270 int* ridx = rhs2.altIndexMem(); 271 272 x.clear(); 273 y.clear(); 274 ssvec.assign(rhs1); 275 n = ssvec.size(); // see altValues(); 276 n = vSolveLeft2(x.altValues(), x.altIndexMem(), svec, sidx, n, 277 y.get_ptr(), rhs2.altValues(), ridx, rn); 278 279 x.setSize(n); 280 281 if(n > 0) 282 x.forceSetup(); 283 else 284 x.unSetup(); 285 286 rhs2.setSize(0); 287 rhs2.forceSetup(); 288 ssvec.setSize(0); 289 ssvec.forceSetup(); 290 291 solveCount++; 292 solveTime->stop(); 293 } 294 295 inline void SLUFactorRational::solveLeft( 296 SSVectorRational& x, 297 VectorRational& y, 298 VectorRational& z, 299 const SVectorRational& rhs1, 300 SSVectorRational& rhs2, 301 SSVectorRational& rhs3) 302 { 303 304 solveTime->start(); 305 306 int n; 307 Rational* svec = ssvec.altValues(); 308 int* sidx = ssvec.altIndexMem(); 309 310 x.clear(); 311 y.clear(); 312 z.clear(); 313 ssvec.assign(rhs1); 314 n = ssvec.size(); // see altValues(); 315 n = vSolveLeft3(x.altValues(), x.altIndexMem(), svec, sidx, n, 316 y.get_ptr(), rhs2.altValues(), rhs2.altIndexMem(), rhs2.size(), 317 z.get_ptr(), rhs3.altValues(), rhs3.altIndexMem(), rhs3.size()); 318 319 x.setSize(n); 320 321 if(n > 0) 322 x.forceSetup(); 323 else 324 x.unSetup(); 325 326 ssvec.setSize(0); 327 ssvec.forceSetup(); 328 329 solveCount++; 330 solveTime->stop(); 331 } 332 333 inline Rational SLUFactorRational::stability() const 334 { 335 336 if(status() != OK) 337 return 0; 338 339 if(maxabs < initMaxabs) 340 return 1; 341 342 return initMaxabs / maxabs; 343 } 344 345 inline std::string SLUFactorRational::statistics() const 346 { 347 std::stringstream s; 348 s << "Factorizations : " << std::setw(10) << getFactorCount() << std::endl 349 << " Time spent : " << std::setw(10) << std::fixed << std::setprecision( 350 2) << getFactorTime() << std::endl 351 << "Solves : " << std::setw(10) << getSolveCount() << std::endl 352 << " Time spent : " << std::setw(10) << getSolveTime() << std::endl; 353 354 return s.str(); 355 } 356 357 inline void SLUFactorRational::changeEta(int idx, SSVectorRational& et) 358 { 359 360 int es = et.size(); // see altValues() 361 update(idx, et.altValues(), et.altIndexMem(), es); 362 et.setSize(0); 363 et.forceSetup(); 364 } 365 366 inline SLUFactorRational::Status SLUFactorRational::change( 367 int idx, 368 const SVectorRational& subst, 369 const SSVectorRational* e) 370 { 371 372 // BH 2005-08-23: The boolean usetup indicates that an "update vector" 373 // has been set up. I suppose that SSVectorRational forest is this 374 // update vector, which is set up by solveRight4update() and 375 // solve2right4update() in order to optimize the basis update. 376 377 if(usetup) 378 { 379 if(l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates 380 { 381 // BH 2005-08-19: The size of a SSVectorRational is the size of the 382 // index set, i.e. the number of nonzeros which is only 383 // defined if the SSVectorRational is set up. Since 384 // SSVectorRational::altValues() calls unSetup() the size needs to be 385 // stored before the following call. 386 int fsize = forest.size(); // see altValues() 387 forestUpdate(idx, forest.altValues(), fsize, forest.altIndexMem()); 388 forest.setSize(0); 389 forest.forceSetup(); 390 } 391 else 392 { 393 assert(l.updateType == ETA); 394 changeEta(idx, eta); 395 } 396 } 397 else if(e != 0) // ETA updates 398 { 399 l.updateType = ETA; 400 updateNoClear(idx, e->values(), e->indexMem(), e->size()); 401 l.updateType = uptype; 402 } 403 else if(l.updateType == FOREST_TOMLIN) // FOREST_TOMLIN updates 404 { 405 assert(0); // probably this part is never called. 406 // forestUpdate() with the last parameter set to NULL should fail. 407 forest = subst; 408 CLUFactorRational::solveLright(forest.altValues()); 409 forestUpdate(idx, forest.altValues(), 0, nullptr); 410 forest.setSize(0); 411 forest.forceSetup(); 412 } 413 else // ETA updates 414 { 415 assert(l.updateType == ETA); 416 vec = subst; 417 eta.clear(); 418 CLUFactorRational::solveRight(eta.altValues(), vec.get_ptr()); 419 changeEta(idx, eta); 420 } 421 422 usetup = false; 423 424 SPxOut::debug(this, "DSLUFA01\tupdated\t\tstability = {}\n", stability()); 425 426 return status(); 427 } 428 429 inline void SLUFactorRational::init() 430 { 431 rowMemMult = 5; /* factor of minimum Memory * #of nonzeros */ 432 colMemMult = 5; /* factor of minimum Memory * #of nonzeros */ 433 lMemMult = 1; /* factor of minimum Memory * #of nonzeros */ 434 435 l.firstUpdate = 0; 436 l.firstUnused = 0; 437 thedim = 0; 438 439 usetup = false; 440 maxabs = 1; 441 initMaxabs = 1; 442 lastThreshold = minThreshold; 443 minStability = SOPLEX_MINSTABILITYRAT; 444 stat = UNLOADED; 445 446 vec.clear(); 447 eta.clear(); 448 ssvec.clear(); 449 forest.clear(); 450 451 u.col.size = 100; 452 l.startSize = 100; 453 454 l.rval.reDim(0); 455 456 if(l.ridx) 457 spx_free(l.ridx); 458 459 if(l.rbeg) 460 spx_free(l.rbeg); 461 462 if(l.rorig) 463 spx_free(l.rorig); 464 465 if(l.rperm) 466 spx_free(l.rperm); 467 468 if(u.row.idx) 469 spx_free(u.row.idx); 470 471 if(u.col.idx) 472 spx_free(u.col.idx); 473 474 if(l.idx) 475 spx_free(l.idx); 476 477 if(l.start) 478 spx_free(l.start); 479 480 if(l.row) 481 spx_free(l.row); 482 483 // init() is used in constructor of SLUFactorRational so we have to 484 // clean up if anything goes wrong here 485 try 486 { 487 u.row.val.reDim(100); 488 spx_alloc(u.row.idx, u.row.val.dim()); 489 spx_alloc(u.col.idx, u.col.size); 490 491 l.val.reDim(100); 492 spx_alloc(l.idx, l.val.dim()); 493 spx_alloc(l.start, l.startSize); 494 spx_alloc(l.row, l.startSize); 495 } 496 catch(const SPxMemoryException& x) 497 { 498 freeAll(); 499 throw x; 500 } 501 } 502 503 inline void SLUFactorRational::clear() 504 { 505 if(stat != UNLOADED) 506 init(); 507 } 508 509 /** assignment used to implement operator=() and copy constructor. 510 * If this is initialised, freeAll() has to be called before. 511 * Class objects from SLUFactorRational are not copied here. 512 */ 513 inline void SLUFactorRational::assign(const SLUFactorRational& old) 514 { 515 unsigned int thediminc; 516 517 // slufactor_rational 518 uptype = old.uptype; 519 minThreshold = old.minThreshold; 520 minStability = old.minStability; 521 lastThreshold = old.lastThreshold; 522 523 // clufactor_rational 524 stat = old.stat; 525 thedim = old.thedim; 526 nzCnt = old.nzCnt; 527 initMaxabs = old.initMaxabs; 528 maxabs = old.maxabs; 529 rowMemMult = old.rowMemMult; 530 colMemMult = old.colMemMult; 531 lMemMult = old.lMemMult; 532 factorCount = old.factorCount; 533 factorTime = TimerFactory::createTimer(old.factorTime->type()); 534 solveTime = TimerFactory::createTimer(old.solveTime->type()); 535 timeLimit = old.timeLimit; 536 537 spx_alloc(row.perm, thedim); 538 spx_alloc(row.orig, thedim); 539 spx_alloc(col.perm, thedim); 540 spx_alloc(col.orig, thedim); 541 542 memcpy(row.perm, old.row.perm, (unsigned int)thedim * sizeof(*row.perm)); 543 memcpy(row.orig, old.row.orig, (unsigned int)thedim * sizeof(*row.orig)); 544 memcpy(col.perm, old.col.perm, (unsigned int)thedim * sizeof(*col.perm)); 545 memcpy(col.orig, old.col.orig, (unsigned int)thedim * sizeof(*col.orig)); 546 diag = old.diag; 547 548 work = vec.get_ptr(); 549 550 /* setup U 551 */ 552 thediminc = (unsigned int)(thedim + 1); 553 u.row.used = old.u.row.used; 554 u.row.val = old.u.row.val; 555 556 spx_alloc(u.row.elem, thedim); 557 spx_alloc(u.row.idx, u.row.val.dim()); 558 spx_alloc(u.row.start, thediminc); 559 spx_alloc(u.row.len, thediminc); 560 spx_alloc(u.row.max, thediminc); 561 562 memcpy(u.row.elem, old.u.row.elem, (unsigned int)thedim * sizeof(*u.row.elem)); 563 memcpy(u.row.idx, old.u.row.idx, (unsigned int)u.row.val.dim() * sizeof(*u.row.idx)); 564 memcpy(u.row.start, old.u.row.start, thediminc * sizeof(*u.row.start)); 565 memcpy(u.row.len, old.u.row.len, thediminc * sizeof(*u.row.len)); 566 memcpy(u.row.max, old.u.row.max, thediminc * sizeof(*u.row.max)); 567 568 // need to make row list ok ? 569 if(thedim > 0 && stat == OK) 570 { 571 u.row.list.idx = old.u.row.list.idx; // .idx neu 572 573 const Dring* oring = &old.u.row.list; 574 Dring* ring = &u.row.list; 575 576 while(oring->next != &old.u.row.list) 577 { 578 ring->next = &u.row.elem[oring->next->idx]; 579 ring->next->prev = ring; 580 oring = oring->next; 581 ring = ring->next; 582 } 583 584 ring->next = &u.row.list; 585 ring->next->prev = ring; 586 } 587 588 u.col.size = old.u.col.size; 589 u.col.used = old.u.col.used; 590 591 spx_alloc(u.col.elem, thedim); 592 spx_alloc(u.col.idx, u.col.size); 593 spx_alloc(u.col.start, thediminc); 594 spx_alloc(u.col.len, thediminc); 595 spx_alloc(u.col.max, thediminc); 596 u.col.val = old.u.col.val; 597 598 memcpy(u.col.elem, old.u.col.elem, (unsigned int)thedim * sizeof(*u.col.elem)); 599 memcpy(u.col.idx, old.u.col.idx, (unsigned int)u.col.size * sizeof(*u.col.idx)); 600 memcpy(u.col.start, old.u.col.start, thediminc * sizeof(*u.col.start)); 601 memcpy(u.col.len, old.u.col.len, thediminc * sizeof(*u.col.len)); 602 memcpy(u.col.max, old.u.col.max, thediminc * sizeof(*u.col.max)); 603 604 // need to make col list ok 605 if(thedim > 0 && stat == OK) 606 { 607 u.col.list.idx = old.u.col.list.idx; // .idx neu 608 609 const Dring* oring = &old.u.col.list; 610 Dring* ring = &u.col.list; 611 612 while(oring->next != &old.u.col.list) 613 { 614 ring->next = &u.col.elem[oring->next->idx]; 615 ring->next->prev = ring; 616 oring = oring->next; 617 ring = ring->next; 618 } 619 620 ring->next = &u.col.list; 621 ring->next->prev = ring; 622 } 623 624 /* Setup L 625 */ 626 l.startSize = old.l.startSize; 627 l.firstUpdate = old.l.firstUpdate; 628 l.firstUnused = old.l.firstUnused; 629 l.updateType = old.l.updateType; 630 631 l.val = old.l.val; 632 spx_alloc(l.idx, l.val.dim()); 633 spx_alloc(l.start, l.startSize); 634 spx_alloc(l.row, l.startSize); 635 636 memcpy(l.idx, old.l.idx, (unsigned int)l.val.dim() * sizeof(*l.idx)); 637 memcpy(l.start, old.l.start, (unsigned int)l.startSize * sizeof(*l.start)); 638 memcpy(l.row, old.l.row, (unsigned int)l.startSize * sizeof(*l.row)); 639 640 if(old.l.rval.dim() != 0) 641 { 642 assert(old.l.ridx != 0); 643 assert(old.l.rbeg != 0); 644 assert(old.l.rorig != 0); 645 assert(old.l.rperm != 0); 646 647 int memsize = l.start[l.firstUpdate]; 648 649 l.rval = old.l.rval; 650 spx_alloc(l.ridx, memsize); 651 spx_alloc(l.rbeg, thediminc); 652 spx_alloc(l.rorig, thedim); 653 spx_alloc(l.rperm, thedim); 654 655 memcpy(l.ridx, old.l.ridx, (unsigned int)memsize * sizeof(*l.ridx)); 656 memcpy(l.rbeg, old.l.rbeg, thediminc * sizeof(*l.rbeg)); 657 memcpy(l.rorig, old.l.rorig, (unsigned int)thedim * sizeof(*l.rorig)); 658 memcpy(l.rperm, old.l.rperm, (unsigned int)thedim * sizeof(*l.rperm)); 659 } 660 else 661 { 662 assert(old.l.ridx == 0); 663 assert(old.l.rbeg == 0); 664 assert(old.l.rorig == 0); 665 assert(old.l.rperm == 0); 666 667 l.rval.reDim(0); 668 l.ridx = 0; 669 l.rbeg = 0; 670 l.rorig = 0; 671 l.rperm = 0; 672 } 673 674 assert(row.perm != 0); 675 assert(row.orig != 0); 676 assert(col.perm != 0); 677 assert(col.orig != 0); 678 679 assert(u.row.elem != 0); 680 assert(u.row.idx != 0); 681 assert(u.row.start != 0); 682 assert(u.row.len != 0); 683 assert(u.row.max != 0); 684 685 assert(u.col.elem != 0); 686 assert(u.col.idx != 0); 687 assert(u.col.start != 0); 688 assert(u.col.len != 0); 689 assert(u.col.max != 0); 690 691 assert(l.idx != 0); 692 assert(l.start != 0); 693 assert(l.row != 0); 694 695 } 696 697 inline void SLUFactorRational::freeAll() 698 { 699 700 if(row.perm) 701 spx_free(row.perm); 702 703 if(row.orig) 704 spx_free(row.orig); 705 706 if(col.perm) 707 spx_free(col.perm); 708 709 if(col.orig) 710 spx_free(col.orig); 711 712 if(u.row.elem) 713 spx_free(u.row.elem); 714 715 if(u.row.idx) 716 spx_free(u.row.idx); 717 718 if(u.row.start) 719 spx_free(u.row.start); 720 721 if(u.row.len) 722 spx_free(u.row.len); 723 724 if(u.row.max) 725 spx_free(u.row.max); 726 727 if(u.col.elem) 728 spx_free(u.col.elem); 729 730 if(u.col.idx) 731 spx_free(u.col.idx); 732 733 if(u.col.start) 734 spx_free(u.col.start); 735 736 if(u.col.len) 737 spx_free(u.col.len); 738 739 if(u.col.max) 740 spx_free(u.col.max); 741 742 if(l.idx) 743 spx_free(l.idx); 744 745 if(l.start) 746 spx_free(l.start); 747 748 if(l.row) 749 spx_free(l.row); 750 751 if(l.ridx) 752 spx_free(l.ridx); 753 754 if(l.rbeg) 755 spx_free(l.rbeg); 756 757 if(l.rorig) 758 spx_free(l.rorig); 759 760 if(l.rperm) 761 spx_free(l.rperm); 762 763 spx_free(solveTime); 764 spx_free(factorTime); 765 } 766 767 inline SLUFactorRational::~SLUFactorRational() 768 { 769 freeAll(); 770 } 771 772 static Rational betterThreshold(Rational th) 773 { 774 assert(th < 1); 775 776 if(10 * th < 1) 777 th *= 10; 778 else if(10 * th < 8) 779 th = (th + 1) / 2; 780 else if(th < 0.999) 781 th = 0.99999; 782 783 assert(th < 1); 784 785 return th; 786 } 787 788 inline SLUFactorRational::Status SLUFactorRational::load(const SVectorRational* matrix[], int dm) 789 { 790 assert(dm >= 0); 791 assert(matrix != 0); 792 793 Rational lastStability = stability(); 794 795 initDR(u.row.list); 796 initDR(u.col.list); 797 798 usetup = false; 799 l.updateType = uptype; 800 l.firstUpdate = 0; 801 l.firstUnused = 0; 802 803 if(dm != thedim) 804 { 805 clear(); 806 807 thedim = dm; 808 vec.reDim(thedim); 809 ssvec.reDim(thedim); 810 eta.reDim(thedim); 811 forest.reDim(thedim); 812 work = vec.get_ptr(); 813 814 spx_realloc(row.perm, thedim); 815 spx_realloc(row.orig, thedim); 816 spx_realloc(col.perm, thedim); 817 spx_realloc(col.orig, thedim); 818 diag.reDim(thedim); 819 820 spx_realloc(u.row.elem, thedim); 821 spx_realloc(u.row.len, thedim + 1); 822 spx_realloc(u.row.max, thedim + 1); 823 spx_realloc(u.row.start, thedim + 1); 824 825 spx_realloc(u.col.elem, thedim); 826 spx_realloc(u.col.len, thedim + 1); 827 spx_realloc(u.col.max, thedim + 1); 828 spx_realloc(u.col.start, thedim + 1); 829 830 l.startSize = thedim + SOPLEX_MAXUPDATES; 831 832 spx_realloc(l.row, l.startSize); 833 spx_realloc(l.start, l.startSize); 834 } 835 // the last factorization was reasonably stable, so we decrease the Markowitz threshold (stored in lastThreshold) in 836 // order favour sparsity 837 else if(lastStability > 2.0 * SOPLEX_MINSTABILITYRAT) 838 { 839 // we reset lastThreshold to its previous value in the sequence minThreshold, betterThreshold(minThreshold), 840 // betterThreshold(betterThreshold(minThreshold)), ... 841 Rational last = minThreshold; 842 Rational better = betterThreshold(last); 843 844 while(better < lastThreshold) 845 { 846 last = better; 847 better = betterThreshold(last); 848 } 849 850 lastThreshold = last; 851 852 // we reset the minimum stability (which might have been decreased below) to ensure that the increased sparsity 853 // does not hurt the stability 854 minStability = 2 * SOPLEX_MINSTABILITYRAT; 855 } 856 857 u.row.list.idx = thedim; 858 u.row.start[thedim] = 0; 859 u.row.max[thedim] = 0; 860 u.row.len[thedim] = 0; 861 862 u.col.list.idx = thedim; 863 u.col.start[thedim] = 0; 864 u.col.max[thedim] = 0; 865 u.col.len[thedim] = 0; 866 867 stat = OK; 868 factor(matrix, lastThreshold); 869 870 871 SPxOut::debug(this, "DSLUFA02 threshold = {}\tstability = {}\tSOPLEX_MINSTABILITYRAT = {}\n", 872 lastThreshold, stability(), SOPLEX_MINSTABILITYRAT); 873 SPX_DEBUG( 874 int i; 875 FILE* fl = fopen("dump.lp", "w"); 876 std::cout << "DSLUFA03 Basis:\n"; 877 int j = 0; 878 879 for(i = 0; i < dim(); ++i) 880 j += matrix[i]->size(); 881 for(i = 0; i < dim(); ++i) 882 { 883 for(j = 0; j < matrix[i]->size(); ++j) 884 fprintf(fl, "%8d %8d %s\n", 885 i + 1, matrix[i]->index(j) + 1, matrix[i]->value(j).str()); 886 } 887 fclose(fl); 888 std::cout << "DSLUFA04 LU-Factors:" << std::endl; 889 dump(); 890 891 std::cout << "DSLUFA05 threshold = " << lastThreshold 892 << "\tstability = " << stability() << std::endl; 893 ) 894 895 assert(isConsistent()); 896 return Status(stat); 897 } 898 899 900 inline bool SLUFactorRational::isConsistent() const 901 { 902 #ifdef ENABLE_CONSISTENCY_CHECKS 903 return CLUFactorRational::isConsistent(); 904 #else 905 return true; 906 #endif 907 } 908 909 inline void SLUFactorRational::dump() const 910 { 911 CLUFactorRational::dump(); 912 } 913 } // namespace soplex 914