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 "soplex/cring.h" 27 28 namespace soplex 29 { 30 /* This number is used to decide wether a value is zero 31 * or was explicitly set to zero. 32 */ 33 /* 34 #define SOPLEX_MARKER 1e-100 35 */ 36 37 static const Real verySparseFactorRat = 0.001; 38 static const Real verySparseFactor4rightRat = 0.2; 39 40 /* generic heap management */ 41 static void enQueueMaxRat(int* heap, int* size, int elem) 42 { 43 int i, j; 44 45 j = (*size)++; 46 47 while(j > 0) 48 { 49 i = (j - 1) / 2; 50 51 if(elem > heap[i]) 52 { 53 heap[j] = heap[i]; 54 j = i; 55 } 56 else 57 break; 58 } 59 60 heap[j] = elem; 61 62 #ifdef SOPLEX_DEBUG 63 64 // no NDEBUG define, since this block is really expensive 65 for(i = 1; i < *size; ++i) 66 for(j = 0; j < i; ++j) 67 assert(heap[i] != heap[j]); 68 69 #endif /* SOPLEX_DEBUG */ 70 } 71 72 static int deQueueMaxRat(int* heap, int* size) 73 { 74 int e, elem; 75 int i, j, s; 76 int e1, e2; 77 78 elem = *heap; 79 e = heap[s = --(*size)]; 80 --s; 81 82 for(j = 0, i = 1; i < s; i = 2 * j + 1) 83 { 84 e1 = heap[i]; 85 e2 = heap[i + 1]; 86 87 if(e1 > e2) 88 { 89 if(e < e1) 90 { 91 heap[j] = e1; 92 j = i; 93 } 94 else 95 { 96 heap[j] = e; 97 return elem; 98 } 99 } 100 else 101 { 102 if(e < e2) 103 { 104 heap[j] = e2; 105 j = i + 1; 106 } 107 else 108 { 109 heap[j] = e; 110 return elem; 111 } 112 } 113 } 114 115 if(i < *size && e < heap[i]) 116 { 117 heap[j] = heap[i]; 118 j = i; 119 } 120 121 heap[j] = e; 122 123 return elem; 124 } 125 126 static void enQueueMinRat(int* heap, int* size, int elem) 127 { 128 int i, j; 129 130 j = (*size)++; 131 132 while(j > 0) 133 { 134 i = (j - 1) / 2; 135 136 if(elem < heap[i]) 137 { 138 heap[j] = heap[i]; 139 j = i; 140 } 141 else 142 break; 143 } 144 145 heap[j] = elem; 146 147 #ifdef SOPLEX_DEBUG 148 149 // no NDEBUG define, since this block is really expensive 150 for(i = 1; i < *size; ++i) 151 for(j = 0; j < i; ++j) 152 assert(heap[i] != heap[j]); 153 154 #endif /* SOPLEX_DEBUG */ 155 } 156 157 static int deQueueMinRat(int* heap, int* size) 158 { 159 int e, elem; 160 int i, j, s; 161 int e1, e2; 162 163 elem = *heap; 164 e = heap[s = --(*size)]; 165 --s; 166 167 for(j = 0, i = 1; i < s; i = 2 * j + 1) 168 { 169 e1 = heap[i]; 170 e2 = heap[i + 1]; 171 172 if(e1 < e2) 173 { 174 if(e > e1) 175 { 176 heap[j] = e1; 177 j = i; 178 } 179 else 180 { 181 heap[j] = e; 182 return elem; 183 } 184 } 185 else 186 { 187 if(e > e2) 188 { 189 heap[j] = e2; 190 j = i + 1; 191 } 192 else 193 { 194 heap[j] = e; 195 return elem; 196 } 197 } 198 } 199 200 if(i < *size && e > heap[i]) 201 { 202 heap[j] = heap[i]; 203 j = i; 204 } 205 206 heap[j] = e; 207 208 return elem; 209 } 210 211 /************************************************************/ 212 inline CLUFactorRational::Temp::Temp() 213 : s_mark(0) 214 , s_cact(0) 215 , stage(0) 216 , pivot_col(0) 217 , pivot_colNZ(0) 218 , pivot_row(0) 219 , pivot_rowNZ(0) 220 {} 221 222 inline void CLUFactorRational::Temp::init(int p_dim) 223 { 224 s_max.reDim(p_dim); 225 spx_realloc(s_cact, p_dim); 226 spx_realloc(s_mark, p_dim); 227 stage = 0; 228 } 229 230 inline void CLUFactorRational::Temp::clear() 231 { 232 if(s_mark != 0) 233 spx_free(s_mark); 234 235 if(s_cact != 0) 236 spx_free(s_cact); 237 238 if(pivot_col != 0) 239 spx_free(pivot_col); 240 241 if(pivot_colNZ != 0) 242 spx_free(pivot_colNZ); 243 244 if(pivot_row != 0) 245 spx_free(pivot_row); 246 247 if(pivot_rowNZ != 0) 248 spx_free(pivot_rowNZ); 249 250 try 251 { 252 s_max.reDim(0); 253 } 254 catch(const SPxMemoryException& x) 255 { 256 throw x; 257 } 258 } 259 260 inline CLUFactorRational::Temp::~Temp() 261 { 262 clear(); 263 } 264 265 /************************************************************/ 266 inline void CLUFactorRational::initPerm() 267 { 268 269 for(int i = 0; i < thedim; ++i) 270 row.orig[i] = row.perm[i] = col.orig[i] = col.perm[i] = -1; 271 } 272 273 /*****************************************************************************/ 274 275 inline void CLUFactorRational::setPivot(const int p_stage, 276 const int p_col, 277 const int p_row, 278 const Rational& val) 279 { 280 assert(row.perm[p_row] < 0); 281 assert(col.perm[p_col] < 0); 282 283 row.orig[p_stage] = p_row; 284 col.orig[p_stage] = p_col; 285 row.perm[p_row] = p_stage; 286 col.perm[p_col] = p_stage; 287 diag[p_row] = 1.0 / val; 288 289 if(spxAbs(diag[p_row]) > maxabs) 290 maxabs = spxAbs(diag[p_row]); 291 } 292 293 /*****************************************************************************/ 294 /* 295 * Perform garbage collection on row file 296 */ 297 inline void CLUFactorRational::packRows() 298 { 299 int n, i, j, l_row; 300 Dring* ring, *list; 301 302 int* l_ridx = u.row.idx; 303 VectorRational& l_rval = u.row.val; 304 int* l_rlen = u.row.len; 305 int* l_rmax = u.row.max; 306 int* l_rbeg = u.row.start; 307 308 n = 0; 309 list = &(u.row.list); 310 311 for(ring = list->next; ring != list; ring = ring->next) 312 { 313 l_row = ring->idx; 314 315 if(l_rbeg[l_row] != n) 316 { 317 do 318 { 319 l_row = ring->idx; 320 i = l_rbeg[l_row]; 321 assert(l_rlen[l_row] <= l_rmax[l_row]); 322 l_rbeg[l_row] = n; 323 l_rmax[l_row] = l_rlen[l_row]; 324 j = i + l_rlen[l_row]; 325 326 for(; i < j; ++i, ++n) 327 { 328 assert(n <= i); 329 l_ridx[n] = l_ridx[i]; 330 l_rval[n] = l_rval[i]; 331 } 332 333 ring = ring->next; 334 } 335 while(ring != list); 336 337 goto terminatePackRows; 338 } 339 340 n += l_rlen[l_row]; 341 342 l_rmax[l_row] = l_rlen[l_row]; 343 } 344 345 terminatePackRows: 346 347 u.row.max[thedim] = 0; 348 u.row.used = n; 349 } 350 351 /*****************************************************************************/ 352 /* 353 * Perform garbage collection on column file 354 */ 355 inline void CLUFactorRational::forestPackColumns() 356 { 357 int n, i, j, colno; 358 Dring* ring, *list; 359 360 VectorRational& cval = u.col.val; 361 int* cidx = u.col.idx; 362 int* clen = u.col.len; 363 int* cmax = u.col.max; 364 int* cbeg = u.col.start; 365 366 n = 0; 367 list = &u.col.list; 368 369 for(ring = list->next; ring != list; ring = ring->next) 370 { 371 colno = ring->idx; 372 373 if(cbeg[colno] != n) 374 { 375 do 376 { 377 colno = ring->idx; 378 i = cbeg[colno]; 379 cbeg[colno] = n; 380 cmax[colno] = clen[colno]; 381 j = i + clen[colno]; 382 383 for(; i < j; ++i) 384 { 385 cval[n] = cval[i]; 386 cidx[n++] = cidx[i]; 387 } 388 389 ring = ring->next; 390 } 391 while(ring != list); 392 393 goto terminatePackColumns; 394 } 395 396 n += clen[colno]; 397 398 cmax[colno] = clen[colno]; 399 } 400 401 terminatePackColumns : 402 403 u.col.used = n; 404 u.col.max[thedim] = 0; 405 } 406 407 /* 408 * Make row of fac large enough to hold len nonzeros. 409 */ 410 inline void CLUFactorRational::remaxRow(int p_row, int len) 411 { 412 assert(u.row.max[p_row] < len); 413 414 if(u.row.elem[p_row].next == &(u.row.list)) /* last in row file */ 415 { 416 int delta = len - u.row.max[p_row]; 417 418 if(delta > u.row.val.dim() - u.row.used) 419 { 420 packRows(); 421 delta = len - u.row.max[p_row]; // packRows() changes u.row.max[] ! 422 423 if(u.row.val.dim() < rowMemMult * u.row.used + len) 424 minRowMem(2 * u.row.used + len); 425 426 /* minRowMem(rowMemMult * u.row.used + len); */ 427 } 428 429 assert(delta <= u.row.val.dim() - u.row.used 430 431 && "ERROR: could not allocate memory for row file"); 432 433 u.row.used += delta; 434 u.row.max[p_row] = len; 435 } 436 else /* row must be moved to end of row file */ 437 { 438 int i, j, k; 439 VectorRational& val = u.row.val; 440 Dring* ring; 441 442 if(len > u.row.val.dim() - u.row.used) 443 { 444 packRows(); 445 446 if(u.row.val.dim() < rowMemMult * u.row.used + len) 447 minRowMem(2 * u.row.used + len); 448 449 /* minRowMem(rowMemMult * u.row.used + len);*/ 450 } 451 452 assert(len <= u.row.val.dim() - u.row.used 453 454 && "ERROR: could not allocate memory for row file"); 455 456 int* idx = u.row.idx; 457 j = u.row.used; 458 i = u.row.start[p_row]; 459 k = u.row.len[p_row] + i; 460 u.row.start[p_row] = j; 461 u.row.used += len; 462 463 u.row.max[u.row.elem[p_row].prev->idx] += u.row.max[p_row]; 464 u.row.max[p_row] = len; 465 removeDR(u.row.elem[p_row]); 466 ring = u.row.list.prev; 467 init2DR(u.row.elem[p_row], *ring); 468 469 for(; i < k; ++i, ++j) 470 { 471 val[j] = val[i]; 472 idx[j] = idx[i]; 473 } 474 } 475 476 assert(u.row.start[u.row.list.prev->idx] + u.row.max[u.row.list.prev->idx] 477 478 == u.row.used); 479 } 480 481 /*************************************************************************/ 482 /* 483 * Perform garbage collection on column file 484 */ 485 inline void CLUFactorRational::packColumns() 486 { 487 int n, i, j, l_col; 488 Dring* ring, *list; 489 490 int* l_cidx = u.col.idx; 491 int* l_clen = u.col.len; 492 int* l_cmax = u.col.max; 493 int* l_cbeg = u.col.start; 494 495 n = 0; 496 list = &(u.col.list); 497 498 for(ring = list->next; ring != list; ring = ring->next) 499 { 500 l_col = ring->idx; 501 502 if(l_cbeg[l_col] != n) 503 { 504 do 505 { 506 l_col = ring->idx; 507 i = l_cbeg[l_col]; 508 l_cbeg[l_col] = n; 509 l_cmax[l_col] = l_clen[l_col]; 510 j = i + l_clen[l_col]; 511 512 for(; i < j; ++i) 513 l_cidx[n++] = l_cidx[i]; 514 515 ring = ring->next; 516 } 517 while(ring != list); 518 519 goto terminatePackColumns; 520 } 521 522 n += l_clen[l_col]; 523 524 l_cmax[l_col] = l_clen[l_col]; 525 } 526 527 terminatePackColumns : 528 529 u.col.used = n; 530 u.col.max[thedim] = 0; 531 } 532 533 /* 534 * Make column col of fac large enough to hold len nonzeros. 535 */ 536 inline void CLUFactorRational::remaxCol(int p_col, int len) 537 { 538 assert(u.col.max[p_col] < len); 539 540 if(u.col.elem[p_col].next == &(u.col.list)) /* last in column file */ 541 { 542 int delta = len - u.col.max[p_col]; 543 544 if(delta > u.col.size - u.col.used) 545 { 546 packColumns(); 547 delta = len - u.col.max[p_col]; 548 549 if(u.col.size < colMemMult * u.col.used + len) 550 minColMem(2 * u.col.used + len); 551 552 /* minColMem(colMemMult * u.col.used + len); */ 553 } 554 555 assert(delta <= u.col.size - u.col.used 556 557 && "ERROR: could not allocate memory for column file"); 558 559 u.col.used += delta; 560 u.col.max[p_col] = len; 561 } 562 563 else /* column must be moved to end of column file */ 564 { 565 int i, j, k; 566 int* idx; 567 Dring* ring; 568 569 if(len > u.col.size - u.col.used) 570 { 571 packColumns(); 572 573 if(u.col.size < colMemMult * u.col.used + len) 574 minColMem(2 * u.col.used + len); 575 576 /* minColMem(colMemMult * u.col.used + len); */ 577 } 578 579 assert(len <= u.col.size - u.col.used 580 581 && "ERROR: could not allocate memory for column file"); 582 583 j = u.col.used; 584 i = u.col.start[p_col]; 585 k = u.col.len[p_col] + i; 586 u.col.start[p_col] = j; 587 u.col.used += len; 588 589 u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col]; 590 u.col.max[p_col] = len; 591 removeDR(u.col.elem[p_col]); 592 ring = u.col.list.prev; 593 init2DR(u.col.elem[p_col], *ring); 594 595 idx = u.col.idx; 596 597 for(; i < k; ++i) 598 idx[j++] = idx[i]; 599 } 600 } 601 602 /* 603 * Make column col of fac large enough to hold len nonzeros. 604 */ 605 inline void CLUFactorRational::forestReMaxCol(int p_col, int len) 606 { 607 assert(u.col.max[p_col] < len); 608 609 if(u.col.elem[p_col].next == &(u.col.list)) /* last in column file */ 610 { 611 int delta = len - u.col.max[p_col]; 612 613 if(delta > u.col.size - u.col.used) 614 { 615 forestPackColumns(); 616 delta = len - u.col.max[p_col]; 617 618 if(u.col.size < colMemMult * u.col.used + len) 619 forestMinColMem(int(colMemMult * u.col.used + len)); 620 } 621 622 assert(delta <= u.col.size - u.col.used 623 624 && "ERROR: could not allocate memory for column file"); 625 626 u.col.used += delta; 627 u.col.max[p_col] = len; 628 } 629 630 else /* column must be moved to end of column file */ 631 { 632 int i, j, k; 633 int* idx = u.col.idx; 634 VectorRational& val = u.col.val; 635 Dring* ring; 636 637 if(len > u.col.size - u.col.used) 638 { 639 forestPackColumns(); 640 641 if(u.col.size < colMemMult * u.col.used + len) 642 forestMinColMem(int(colMemMult * u.col.used + len)); 643 } 644 645 assert(len <= u.col.size - u.col.used 646 647 && "ERROR: could not allocate memory for column file"); 648 649 j = u.col.used; 650 i = u.col.start[p_col]; 651 k = u.col.len[p_col] + i; 652 u.col.start[p_col] = j; 653 u.col.used += len; 654 655 u.col.max[u.col.elem[p_col].prev->idx] += u.col.max[p_col]; 656 u.col.max[p_col] = len; 657 removeDR(u.col.elem[p_col]); 658 ring = u.col.list.prev; 659 init2DR(u.col.elem[p_col], *ring); 660 661 for(; i < k; ++i) 662 { 663 val[j] = val[i]; 664 idx[j++] = idx[i]; 665 } 666 } 667 } 668 669 /*****************************************************************************/ 670 671 /** 672 * \brief Performs the Forrest-Tomlin update of the LU factorization. 673 * 674 * BH: I suppose this is implemented as described in UH Suhl, LM Suhl: A fast LU 675 * update for linear programming, Annals of OR 43, p. 33-47, 1993. 676 * 677 * @param p_col Index of basis column to replace. 678 * @param p_work Dense vector to substitute in the basis. 679 * @param num Number of nonzeros in vector represented by p_work. 680 * @param nonz Indices of nonzero elements in vector p_work. 681 * 682 * The parameters num and nonz are used for the following optimization: If both 683 * are nonzero, indices of the nonzero elements provided in nonz (num giving 684 * their number) allow to access only those nonzero entries. Otherwise we have 685 * to go through the entire dense vector element by element. 686 * 687 * After copying p_work into U, p_work is used to expand the row r, which is 688 * needed to restore the triangular structure of U. 689 * 690 * Also num and nonz are used to maintain a heap if there are only very few 691 * nonzeros to be eliminated. This is plainly wrong if the method is called with 692 * nonz==0, see todo at the corresponding place below. 693 * 694 * @throw SPxStatusException if the loaded matrix is singular 695 * 696 * @todo Use an extra member variable as a buffer for working with the dense 697 * row instead of misusing p_work. I think that should be as efficient and 698 * much cleaner. 699 */ 700 701 inline void CLUFactorRational::forestUpdate(int p_col, Rational* p_work, int num, int* nonz) 702 { 703 int i, j, k, h, m, n; 704 int ll, c, r, rowno; 705 Rational x; 706 707 int* lbeg = l.start; 708 709 VectorRational& cval = u.col.val; 710 int* cidx = u.col.idx; 711 int* cmax = u.col.max; 712 int* clen = u.col.len; 713 int* cbeg = u.col.start; 714 715 VectorRational& rval = u.row.val; 716 int* ridx = u.row.idx; 717 int* rmax = u.row.max; 718 int* rlen = u.row.len; 719 int* rbeg = u.row.start; 720 721 int* rperm = row.perm; 722 int* rorig = row.orig; 723 int* cperm = col.perm; 724 int* corig = col.orig; 725 726 Rational l_maxabs = maxabs; 727 int dim = thedim; 728 729 /* Remove column p_col from U 730 */ 731 j = cbeg[p_col]; 732 i = clen[p_col]; 733 nzCnt -= i; 734 735 for(i += j - 1; i >= j; --i) 736 { 737 m = cidx[i]; // remove column p_col from row m 738 k = rbeg[m]; 739 h = --(rlen[m]) + k; // decrease length of row m 740 741 while(ridx[k] != p_col) 742 ++k; 743 744 assert(k <= h); // k is the position of p_col, h is last position 745 746 ridx[k] = ridx[h]; // store last index at the position of p_col 747 748 rval[k] = rval[h]; 749 } 750 751 /* Insert new vector column p_col thereby determining the highest permuted 752 * row index r. 753 * 754 * Distinguish between optimized call (num > 0, nonz != 0) and 755 * non-optimized one. 756 */ 757 assert(num); // otherwise the assert( nonz != 0 ) below should fail 758 759 if(num) 760 { 761 // Optimized call. 762 assert(nonz != 0); 763 764 clen[p_col] = 0; 765 766 if(num > cmax[p_col]) 767 forestReMaxCol(p_col, num); 768 769 cidx = u.col.idx; 770 771 cval = u.col.val; 772 773 k = cbeg[p_col]; 774 775 r = 0; 776 777 for(j = 0; j < num; ++j) 778 { 779 i = nonz[j]; 780 x = p_work[i]; 781 p_work[i] = 0; 782 783 if(x != 0) 784 { 785 if(spxAbs(x) > l_maxabs) 786 l_maxabs = spxAbs(x); 787 788 /* insert to column file */ 789 assert(k - cbeg[p_col] < cmax[p_col]); 790 791 cval[k] = x; 792 793 cidx[k++] = i; 794 795 /* insert to row file */ 796 if(rmax[i] <= rlen[i]) 797 { 798 remaxRow(i, rlen[i] + 1); 799 rval = u.row.val; 800 ridx = u.row.idx; 801 } 802 803 h = rbeg[i] + (rlen[i])++; 804 805 rval[h] = x; 806 ridx[h] = p_col; 807 808 /* check permuted row index */ 809 810 if(rperm[i] > r) 811 r = rperm[i]; 812 } 813 } 814 815 nzCnt += (clen[p_col] = k - cbeg[p_col]); 816 } 817 else 818 { 819 // Non-optimized call: We have to access all elements of p_work. 820 assert(nonz == 0); 821 822 /* 823 * clen[col] = 0; 824 * reMaxCol(fac, col, dim); 825 */ 826 cidx = u.col.idx; 827 cval = u.col.val; 828 k = cbeg[p_col]; 829 j = k + cmax[p_col]; 830 r = 0; 831 832 for(i = 0; i < dim; ++i) 833 { 834 x = p_work[i]; 835 p_work[i] = 0; 836 837 if(x != 0) 838 { 839 if(spxAbs(x) > l_maxabs) 840 l_maxabs = spxAbs(x); 841 842 /* insert to column file */ 843 if(k >= j) 844 { 845 clen[p_col] = k - cbeg[p_col]; 846 forestReMaxCol(p_col, dim - i); 847 cidx = u.col.idx; 848 cval = u.col.val; 849 k = cbeg[p_col]; 850 j = k + cmax[p_col]; 851 k += clen[p_col]; 852 } 853 854 assert(k - cbeg[p_col] < cmax[p_col]); 855 856 cval[k] = x; 857 cidx[k++] = i; 858 859 /* insert to row file */ 860 861 if(rmax[i] <= rlen[i]) 862 { 863 remaxRow(i, rlen[i] + 1); 864 rval = u.row.val; 865 ridx = u.row.idx; 866 } 867 868 h = rbeg[i] + (rlen[i])++; 869 870 rval[h] = x; 871 ridx[h] = p_col; 872 873 /* check permuted row index */ 874 875 if(rperm[i] > r) 876 r = rperm[i]; 877 } 878 } 879 880 nzCnt += (clen[p_col] = k - cbeg[p_col]); 881 882 if(cbeg[p_col] + cmax[p_col] == u.col.used) 883 { 884 u.col.used -= cmax[p_col]; 885 cmax[p_col] = clen[p_col]; 886 u.col.used += cmax[p_col]; 887 } 888 } 889 890 c = cperm[p_col]; 891 892 if(r > c) /* Forest Tomlin update */ 893 { 894 /* update permutations 895 */ 896 j = rorig[c]; 897 898 // memmove is more efficient than a for loop 899 // for ( i = c; i < r; ++i ) 900 // rorig[i] = rorig[i + 1]; 901 memmove(&rorig[c], &rorig[c + 1], (unsigned int)(r - c) * sizeof(int)); 902 903 rorig[r] = j; 904 905 for(i = c; i <= r; ++i) 906 rperm[rorig[i]] = i; 907 908 j = corig[c]; 909 910 // memmove is more efficient than a for loop 911 // for ( i = c; i < r; ++i ) 912 // corig[i] = corig[i + 1]; 913 memmove(&corig[c], &corig[c + 1], (unsigned int)(r - c) * sizeof(int)); 914 915 corig[r] = j; 916 917 for(i = c; i <= r; ++i) 918 cperm[corig[i]] = i; 919 920 921 rowno = rorig[r]; 922 923 j = rbeg[rowno]; 924 925 i = rlen[rowno]; 926 927 nzCnt -= i; 928 929 if(i < verySparseFactorRat * (dim - c)) // few nonzeros to be eliminated 930 { 931 /** 932 * The following assert is obviously violated if this method is called 933 * with nonzero==0. 934 * 935 * @todo Use an extra member variable as a buffer for the heap instead of 936 * misusing nonz and num. The method enQueueMinRat() seems to 937 * sort the nonzeros or something, for which it only needs 938 * some empty vector of size num. 939 */ 940 assert(nonz != 0); 941 942 /* move row r from U to p_work 943 */ 944 num = 0; 945 946 for(i += j - 1; i >= j; --i) 947 { 948 k = ridx[i]; 949 p_work[k] = rval[i]; 950 enQueueMinRat(nonz, &num, cperm[k]); 951 m = --(clen[k]) + cbeg[k]; 952 953 for(h = m; cidx[h] != rowno; --h) 954 ; 955 956 cidx[h] = cidx[m]; 957 958 cval[h] = cval[m]; 959 } 960 961 962 /* Eliminate row r from U to L file 963 */ 964 ll = makeLvec(r - c, rowno); 965 966 VectorRational& lval = l.val; 967 int* lidx = l.idx; 968 969 assert((num == 0) || (nonz != 0)); 970 971 /* for(i = c; i < r; ++i) */ 972 while(num) 973 { 974 #ifndef NDEBUG 975 // The numbers seem to be often 1e-100, is this ok ? 976 977 for(i = 0; i < num; ++i) 978 assert(p_work[corig[nonz[i]]] != 0); 979 980 #endif // NDEBUG 981 i = deQueueMinRat(nonz, &num); 982 983 if(i == r) 984 break; 985 986 k = corig[i]; 987 988 assert(p_work[k] != 0); 989 990 n = rorig[i]; 991 992 x = p_work[k] * diag[n]; 993 994 lidx[ll] = n; 995 996 lval[ll] = x; 997 998 p_work[k] = 0; 999 1000 ll++; 1001 1002 if(spxAbs(x) > l_maxabs) 1003 l_maxabs = spxAbs(x); 1004 1005 j = rbeg[n]; 1006 1007 m = rlen[n] + j; 1008 1009 for(; j < m; ++j) 1010 { 1011 int jj = ridx[j]; 1012 Rational y = p_work[jj]; 1013 1014 if(y == 0) 1015 enQueueMinRat(nonz, &num, cperm[jj]); 1016 1017 y -= x * rval[j]; 1018 1019 //p_work[jj] = y + (( y == 0 ) ? SOPLEX_MARKER : 0 ); 1020 p_work[jj] = y; 1021 } 1022 } 1023 1024 if(lbeg[l.firstUnused - 1] == ll) 1025 (l.firstUnused)--; 1026 else 1027 lbeg[l.firstUnused] = ll; 1028 1029 1030 /* Set diagonal value 1031 */ 1032 if(i != r) 1033 { 1034 stat = SLinSolverRational::SINGULAR; 1035 throw SPxStatusException("XFORE01 The loaded matrix is singular"); 1036 } 1037 1038 k = corig[r]; 1039 1040 x = p_work[k]; 1041 diag[rowno] = 1 / x; 1042 p_work[k] = 0; 1043 1044 1045 /* make row large enough to fit all nonzeros. 1046 */ 1047 1048 if(rmax[rowno] < num) 1049 { 1050 rlen[rowno] = 0; 1051 remaxRow(rowno, num); 1052 rval = u.row.val; 1053 ridx = u.row.idx; 1054 } 1055 1056 nzCnt += num; 1057 1058 /* Insert work to updated row thereby clearing work; 1059 */ 1060 n = rbeg[rowno]; 1061 1062 for(i = 0; i < num; ++i) 1063 { 1064 j = corig[nonz[i]]; 1065 x = p_work[j]; 1066 1067 // BH 2005-08-24: This if is very important. It may well happen that 1068 // during the elimination of row r a nonzero elements cancels out 1069 // and becomes zero. This would lead to an infinite loop in the 1070 // above elimination code, since the corresponding column index would 1071 // be enqueued for further elimination again and agian. 1072 1073 if(x != 0) 1074 { 1075 if(spxAbs(x) > l_maxabs) 1076 l_maxabs = spxAbs(x); 1077 1078 ridx[n] = j; 1079 1080 rval[n] = x; 1081 1082 p_work[j] = 0; 1083 1084 ++n; 1085 1086 if(clen[j] >= cmax[j]) 1087 { 1088 forestReMaxCol(j, clen[j] + 1); 1089 cidx = u.col.idx; 1090 cval = u.col.val; 1091 } 1092 1093 cval[cbeg[j] + clen[j]] = x; 1094 1095 cidx[cbeg[j] + clen[j]++] = rowno; 1096 } 1097 } 1098 1099 rlen[rowno] = n - rbeg[rowno]; 1100 } 1101 else /* few nonzeros to be eliminated */ 1102 { 1103 /* move row r from U to p_work 1104 */ 1105 for(i += j - 1; i >= j; --i) 1106 { 1107 k = ridx[i]; 1108 p_work[k] = rval[i]; 1109 m = --(clen[k]) + cbeg[k]; 1110 1111 for(h = m; cidx[h] != rowno; --h) 1112 ; 1113 1114 cidx[h] = cidx[m]; 1115 1116 cval[h] = cval[m]; 1117 } 1118 1119 1120 /* Eliminate row r from U to L file 1121 */ 1122 ll = makeLvec(r - c, rowno); 1123 1124 VectorRational& lval = l.val; 1125 int* lidx = l.idx; 1126 1127 for(i = c; i < r; ++i) 1128 { 1129 k = corig[i]; 1130 1131 if(p_work[k] != 0) 1132 { 1133 n = rorig[i]; 1134 x = p_work[k] * diag[n]; 1135 lidx[ll] = n; 1136 lval[ll] = x; 1137 p_work[k] = 0; 1138 ll++; 1139 1140 if(spxAbs(x) > l_maxabs) 1141 l_maxabs = spxAbs(x); 1142 1143 j = rbeg[n]; 1144 1145 m = rlen[n] + j; 1146 1147 for(; j < m; ++j) 1148 p_work[ridx[j]] -= x * rval[j]; 1149 } 1150 } 1151 1152 if(lbeg[l.firstUnused - 1] == ll) 1153 (l.firstUnused)--; 1154 else 1155 lbeg[l.firstUnused] = ll; 1156 1157 1158 /* Set diagonal value 1159 */ 1160 k = corig[r]; 1161 1162 x = p_work[k]; 1163 1164 if(x == 0) 1165 { 1166 stat = SLinSolverRational::SINGULAR; 1167 throw SPxStatusException("XFORE02 The loaded matrix is singular"); 1168 // return; 1169 } 1170 1171 diag[rowno] = 1 / x; 1172 1173 p_work[k] = 0; 1174 1175 1176 /* count remaining nonzeros in work and make row large enough 1177 * to fit them all. 1178 */ 1179 n = 0; 1180 1181 for(i = r + 1; i < dim; ++i) 1182 if(p_work[corig[i]] != 0) 1183 n++; 1184 1185 if(rmax[rowno] < n) 1186 { 1187 rlen[rowno] = 0; 1188 remaxRow(rowno, n); 1189 rval = u.row.val; 1190 ridx = u.row.idx; 1191 } 1192 1193 nzCnt += n; 1194 1195 /* Insert p_work to updated row thereby clearing p_work; 1196 */ 1197 n = rbeg[rowno]; 1198 1199 for(i = r + 1; i < dim; ++i) 1200 { 1201 j = corig[i]; 1202 x = p_work[j]; 1203 1204 if(x != 0) 1205 { 1206 if(spxAbs(x) > l_maxabs) 1207 l_maxabs = spxAbs(x); 1208 1209 ridx[n] = j; 1210 1211 rval[n] = x; 1212 1213 p_work[j] = 0; 1214 1215 ++n; 1216 1217 if(clen[j] >= cmax[j]) 1218 { 1219 forestReMaxCol(j, clen[j] + 1); 1220 cidx = u.col.idx; 1221 cval = u.col.val; 1222 } 1223 1224 cval[cbeg[j] + clen[j]] = x; 1225 1226 cidx[cbeg[j] + clen[j]++] = rowno; 1227 } 1228 } 1229 1230 rlen[rowno] = n - rbeg[rowno]; 1231 } 1232 } 1233 1234 else if(r == c) 1235 { 1236 /* Move diagonal element to diag. Note, that it must be the last 1237 * element, since it has just been inserted above. 1238 */ 1239 rowno = rorig[r]; 1240 i = rbeg[rowno] + --(rlen[rowno]); 1241 diag[rowno] = 1 / rval[i]; 1242 1243 for(j = i = --(clen[p_col]) + cbeg[p_col]; cidx[i] != rowno; --i) 1244 ; 1245 1246 cidx[i] = cidx[j]; 1247 1248 cval[i] = cval[j]; 1249 } 1250 else /* r < c */ 1251 { 1252 stat = SLinSolverRational::SINGULAR; 1253 throw SPxStatusException("XFORE03 The loaded matrix is singular"); 1254 // return; 1255 } 1256 1257 maxabs = l_maxabs; 1258 1259 assert(isConsistent()); 1260 stat = SLinSolverRational::OK; 1261 } 1262 1263 inline void CLUFactorRational::update(int p_col, Rational* p_work, const int* p_idx, int num) 1264 { 1265 int ll, i, j; 1266 Rational x, rezi; 1267 1268 assert(p_work[p_col] != 0); 1269 rezi = 1 / p_work[p_col]; 1270 p_work[p_col] = 0; 1271 1272 ll = makeLvec(num, p_col); 1273 // ll = fac->makeLvec(num, col); 1274 VectorRational& lval = l.val; 1275 int* lidx = l.idx; 1276 1277 for(i = num - 1; (j = p_idx[i]) != p_col; --i) 1278 { 1279 lidx[ll] = j; 1280 lval[ll] = rezi * p_work[j]; 1281 p_work[j] = 0; 1282 ++ll; 1283 } 1284 1285 lidx[ll] = p_col; 1286 1287 lval[ll] = 1 - rezi; 1288 ++ll; 1289 1290 for(--i; i >= 0; --i) 1291 { 1292 j = p_idx[i]; 1293 lidx[ll] = j; 1294 lval[ll] = x = rezi * p_work[j]; 1295 p_work[j] = 0; 1296 ++ll; 1297 1298 if(spxAbs(x) > maxabs) 1299 maxabs = spxAbs(x); 1300 } 1301 1302 stat = SLinSolverRational::OK; 1303 } 1304 1305 inline void CLUFactorRational::updateNoClear( 1306 int p_col, 1307 const Rational* p_work, 1308 const int* p_idx, 1309 int num) 1310 { 1311 int ll, i, j; 1312 Rational x, rezi; 1313 1314 assert(p_work[p_col] != 0); 1315 rezi = 1 / p_work[p_col]; 1316 ll = makeLvec(num, p_col); 1317 //ll = fac->makeLvec(num, col); 1318 VectorRational& lval = l.val; 1319 int* lidx = l.idx; 1320 1321 for(i = num - 1; (j = p_idx[i]) != p_col; --i) 1322 { 1323 lidx[ll] = j; 1324 lval[ll] = rezi * p_work[j]; 1325 ++ll; 1326 } 1327 1328 lidx[ll] = p_col; 1329 1330 lval[ll] = 1 - rezi; 1331 ++ll; 1332 1333 for(--i; i >= 0; --i) 1334 { 1335 j = p_idx[i]; 1336 lidx[ll] = j; 1337 lval[ll] = x = rezi * p_work[j]; 1338 ++ll; 1339 1340 if(spxAbs(x) > maxabs) 1341 maxabs = spxAbs(x); 1342 } 1343 1344 stat = SLinSolverRational::OK; 1345 } 1346 1347 /*****************************************************************************/ 1348 /* 1349 * Temporary data structures. 1350 */ 1351 1352 /* 1353 * For the i=th column the situation might look like this: 1354 * 1355 * \begin{verbatim} 1356 * idx = ....................iiiIIIIII-----.............. 1357 * cbeg[i] = ^ 1358 * cact[i] = +----+ 1359 * clen[i] = +-------+ 1360 * cmax[i] = +------------+ 1361 * 1362 * Indices clen[i]-cbeg[i]: ^^^ 1363 * \end{verbatim} 1364 * belong to column i, but have already been pivotal and don't belong to 1365 * the active submatrix. 1366 */ 1367 1368 /****************************************************************************/ 1369 /* 1370 * Initialize row and column file of working matrix and 1371 * mark column singletons. 1372 */ 1373 inline void CLUFactorRational::initFactorMatrix(const SVectorRational** vec) 1374 { 1375 1376 Rational x; 1377 int m; 1378 int tot; 1379 Dring* rring, *lastrring; 1380 Dring* cring, *lastcring; 1381 const SVectorRational* psv; 1382 int* sing = temp.s_mark; 1383 1384 /* Initialize: 1385 * - column file thereby remembering column singletons in |sing|. 1386 * - nonzeros counts per row 1387 * - total number of nonzeros 1388 */ 1389 1390 for(int i = 0; i < thedim; i++) 1391 u.row.max[i] = u.row.len[i] = 0; 1392 1393 tot = 0; 1394 1395 for(int i = 0; i < thedim; i++) 1396 { 1397 int k; 1398 1399 psv = vec[i]; 1400 k = psv->size(); 1401 1402 if(k > 1) 1403 { 1404 tot += k; 1405 1406 for(int j = 0; j < k; ++j) 1407 u.row.max[psv->index(j)]++; 1408 } 1409 else if(k == 0) 1410 { 1411 stat = SLinSolverRational::SINGULAR; 1412 return; 1413 } 1414 } 1415 1416 /* Resize nonzero memory if necessary 1417 */ 1418 minRowMem(int(rowMemMult * tot)); 1419 1420 minColMem(int(colMemMult * tot)); 1421 1422 minLMem(int(lMemMult * tot)); 1423 1424 1425 /* Initialize: 1426 * - row ring lists 1427 * - row vectors in file 1428 * - column ring lists 1429 */ 1430 u.row.start[0] = 0; 1431 1432 rring = u.row.elem; 1433 1434 lastrring = &(u.row.list); 1435 1436 lastrring->idx = thedim; 1437 1438 lastrring->next = rring; 1439 1440 cring = u.col.elem; 1441 1442 lastcring = &(u.col.list); 1443 1444 lastcring->idx = thedim; 1445 1446 lastcring->next = cring; 1447 1448 m = 0; 1449 1450 for(int i = 0; i < thedim; i++) 1451 { 1452 u.row.start[i] = m; 1453 m += u.row.max[i]; 1454 1455 rring->idx = i; 1456 rring->prev = lastrring; 1457 lastrring->next = rring; 1458 lastrring = rring; 1459 ++rring; 1460 1461 cring->idx = i; 1462 cring->prev = lastcring; 1463 lastcring->next = cring; 1464 lastcring = cring; 1465 ++cring; 1466 } 1467 1468 u.row.start[thedim] = 0; 1469 1470 u.row.max[thedim] = 0; 1471 u.row.used = m; 1472 1473 lastrring->next = &(u.row.list); 1474 lastrring->next->prev = lastrring; 1475 1476 lastcring->next = &(u.col.list); 1477 lastcring->next->prev = lastcring; 1478 1479 /* Copy matrix to row and column file 1480 * excluding and marking column singletons! 1481 */ 1482 m = 0; 1483 temp.stage = 0; 1484 1485 initMaxabs = 0; 1486 1487 for(int i = 0; i < thedim; i++) 1488 { 1489 int nnonzeros; 1490 1491 psv = vec[i]; 1492 u.col.start[i] = m; 1493 1494 /* check whether number of nonzeros above tolerance is 0, 1 or >= 2 */ 1495 nnonzeros = 0; 1496 1497 for(int j = 0; j < psv->size() && nnonzeros <= 1; j++) 1498 { 1499 if(psv->value(j) != 0) 1500 nnonzeros++; 1501 } 1502 1503 /* basis is singular due to empty column */ 1504 if(nnonzeros == 0) 1505 { 1506 stat = SLinSolverRational::SINGULAR; 1507 return; 1508 } 1509 1510 /* exclude column singletons */ 1511 else if(nnonzeros == 1) 1512 { 1513 int j; 1514 1515 /* find nonzero */ 1516 1517 for(j = 0; psv->value(j) == 0; j++) 1518 ; 1519 1520 assert(j < psv->size()); 1521 1522 /* basis is singular due to two linearly dependent column singletons */ 1523 if(row.perm[psv->index(j)] >= 0) 1524 { 1525 stat = SLinSolverRational::SINGULAR; 1526 return; 1527 } 1528 1529 /* update maximum absolute nonzero value */ 1530 x = psv->value(j); 1531 1532 if(spxAbs(x) > initMaxabs) 1533 initMaxabs = spxAbs(x); 1534 1535 /* permute to front and mark as singleton */ 1536 setPivot(temp.stage, i, psv->index(j), x); 1537 1538 sing[temp.stage] = i; 1539 1540 temp.stage++; 1541 1542 /* set column length to zero */ 1543 temp.s_cact[i] = u.col.len[i] = u.col.max[i] = 0; 1544 } 1545 1546 /* add to active matrix if not a column singleton */ 1547 else 1548 { 1549 int end; 1550 int k; 1551 1552 /* go through all nonzeros in column */ 1553 assert(nnonzeros >= 2); 1554 nnonzeros = 0; 1555 1556 for(int j = 0; j < psv->size(); j++) 1557 { 1558 x = psv->value(j); 1559 1560 if(x != 0) 1561 { 1562 /* add to column array */ 1563 k = psv->index(j); 1564 u.col.idx[m] = k; 1565 m++; 1566 1567 /* add to row array */ 1568 end = u.row.start[k] + u.row.len[k]; 1569 u.row.idx[end] = i; 1570 u.row.val[end] = x; 1571 u.row.len[k]++; 1572 1573 /* update maximum absolute nonzero value */ 1574 1575 if(spxAbs(x) > initMaxabs) 1576 initMaxabs = spxAbs(x); 1577 1578 nnonzeros++; 1579 } 1580 } 1581 1582 assert(nnonzeros >= 2); 1583 1584 /* set column length */ 1585 temp.s_cact[i] = u.col.len[i] = u.col.max[i] = nnonzeros; 1586 } 1587 } 1588 1589 u.col.used = m; 1590 } 1591 1592 1593 1594 /*****************************************************************************/ 1595 /* 1596 * Remove column singletons 1597 */ 1598 1599 inline void CLUFactorRational::colSingletons() 1600 { 1601 int i, j, k, n; 1602 int len; 1603 int p_col, p_row, newrow; 1604 int* idx; 1605 int* rorig = row.orig; 1606 int* rperm = row.perm; 1607 int* sing = temp.s_mark; 1608 1609 1610 /* Iteratively update column counts due to removed column singletons 1611 * thereby removing new arising columns singletons 1612 * and computing the index of the first row singleton (-1) 1613 * until no more can be found. 1614 */ 1615 1616 for(i = 0; i < temp.stage; ++i) 1617 { 1618 p_row = rorig[i]; 1619 assert(p_row >= 0); 1620 idx = &(u.row.idx[u.row.start[p_row]]); 1621 len = u.row.len[p_row]; 1622 1623 for(j = 0; j < len; ++j) 1624 { 1625 /* Move pivotal nonzeros to front of column. 1626 */ 1627 p_col = idx[j]; 1628 assert(temp.s_cact[p_col] > 0); 1629 1630 n = u.col.start[p_col] + u.col.len[p_col] - temp.s_cact[p_col]; 1631 1632 for(k = n; u.col.idx[k] != p_row; ++k) 1633 ; 1634 1635 assert(k < u.col.start[p_col] + u.col.len[p_col]); 1636 1637 u.col.idx[k] = u.col.idx[n]; 1638 1639 u.col.idx[n] = p_row; 1640 1641 n = --(temp.s_cact[p_col]); /* column nonzeros of ACTIVE matrix */ 1642 1643 if(n == 1) /* Here is another singleton */ 1644 { 1645 newrow = u.col.idx[--u.col.len[p_col] + u.col.start[p_col]]; 1646 1647 /* Ensure, matrix not singular 1648 */ 1649 1650 if(rperm[newrow] >= 0) 1651 { 1652 stat = SLinSolverRational::SINGULAR; 1653 return; 1654 } 1655 1656 /* Find singleton in row. 1657 */ 1658 n = u.row.start[newrow] + (--(u.row.len[newrow])); 1659 1660 for(k = n; u.row.idx[k] != p_col; --k) 1661 ; 1662 1663 /* Remove singleton from column. 1664 */ 1665 setPivot(temp.stage, p_col, newrow, u.row.val[k]); 1666 1667 sing[temp.stage++] = p_col; 1668 1669 /* Move pivot element to diag. 1670 */ 1671 u.row.val[k] = u.row.val[n]; 1672 1673 u.row.idx[k] = u.row.idx[n]; 1674 } 1675 else if(n == 0) 1676 { 1677 stat = SLinSolverRational::SINGULAR; 1678 return; 1679 } 1680 } 1681 } 1682 1683 assert(temp.stage <= thedim); 1684 } 1685 1686 1687 /*****************************************************************************/ 1688 /* 1689 * Remove row singletons 1690 */ 1691 inline void CLUFactorRational::rowSingletons() 1692 { 1693 Rational pval; 1694 int i, j, k, ll, r; 1695 int p_row, p_col, len, rs, lk; 1696 int* idx; 1697 int* rperm = row.perm; 1698 int* sing = temp.s_mark; 1699 1700 /* Mark row singletons 1701 */ 1702 rs = temp.stage; 1703 1704 for(i = 0; i < thedim; ++i) 1705 { 1706 if(rperm[i] < 0 && u.row.len[i] == 1) 1707 sing[temp.stage++] = i; 1708 } 1709 1710 /* Eliminate row singletons 1711 * thereby marking newly arising ones 1712 * until no more can be found. 1713 */ 1714 for(; rs < temp.stage; ++rs) 1715 { 1716 /* Move pivot element from row file to diag 1717 */ 1718 p_row = sing[rs]; 1719 j = u.row.start[p_row]; 1720 p_col = u.row.idx[j]; 1721 pval = u.row.val[j]; 1722 setPivot(rs, p_col, p_row, pval); 1723 u.row.len[p_row] = 0; 1724 1725 /* Remove pivot column form workingmatrix 1726 * thereby building up L vector. 1727 */ 1728 idx = &(u.col.idx[u.col.start[p_col]]); 1729 i = temp.s_cact[p_col]; /* nr. nonzeros of new L vector */ 1730 lk = makeLvec(i - 1, p_row); 1731 len = u.col.len[p_col]; 1732 i = (u.col.len[p_col] -= i); /* remove pivot column from U */ 1733 1734 for(; i < len; ++i) 1735 { 1736 r = idx[i]; 1737 1738 if(r != p_row) 1739 { 1740 /* Find pivot column in row. 1741 */ 1742 ll = --(u.row.len[r]); 1743 k = u.row.start[r] + ll; 1744 1745 for(j = k; u.row.idx[j] != p_col; --j) 1746 ; 1747 1748 assert(k >= u.row.start[r]); 1749 1750 /* Initialize L vector 1751 */ 1752 l.idx[lk] = r; 1753 1754 l.val[lk] = u.row.val[j] / pval; 1755 1756 ++lk; 1757 1758 /* Remove pivot column from row. 1759 */ 1760 u.row.idx[j] = u.row.idx[k]; 1761 1762 u.row.val[j] = u.row.val[k]; 1763 1764 /* Check new row length. 1765 */ 1766 if(ll == 1) 1767 sing[temp.stage++] = r; 1768 else if(ll == 0) 1769 { 1770 stat = SLinSolverRational::SINGULAR; 1771 return; 1772 } 1773 } 1774 } 1775 } 1776 } 1777 1778 1779 /*****************************************************************************/ 1780 /* 1781 * Init nonzero number Ring lists 1782 * and required entries of arrays max and mark 1783 */ 1784 1785 inline void CLUFactorRational::initFactorRings() 1786 { 1787 int i; 1788 int* rperm = row.perm; 1789 int* cperm = col.perm; 1790 CLUFactorRational::Pring* ring; 1791 1792 assert(thedim >= 0); 1793 spx_alloc(temp.pivot_col, thedim + 1); 1794 spx_alloc(temp.pivot_colNZ, thedim + 1); 1795 spx_alloc(temp.pivot_row, thedim + 1); 1796 spx_alloc(temp.pivot_rowNZ, thedim + 1); 1797 1798 for(i = thedim - temp.stage; i >= 0; --i) 1799 { 1800 initDR(temp.pivot_colNZ[i]); 1801 initDR(temp.pivot_rowNZ[i]); 1802 } 1803 1804 for(i = 0; i < thedim; ++i) 1805 { 1806 if(rperm[i] < 0) 1807 { 1808 if(u.row.len[i] <= 0) 1809 { 1810 stat = SLinSolverRational::SINGULAR; 1811 return; 1812 } 1813 1814 ring = &(temp.pivot_rowNZ[u.row.len[i]]); 1815 1816 init2DR(temp.pivot_row[i], *ring); 1817 temp.pivot_row[i].idx = i; 1818 temp.s_max[i] = -1; 1819 } 1820 1821 if(cperm[i] < 0) 1822 { 1823 if(temp.s_cact[i] <= 0) 1824 { 1825 stat = SLinSolverRational::SINGULAR; 1826 return; 1827 } 1828 1829 ring = &(temp.pivot_colNZ[temp.s_cact[i]]); 1830 1831 init2DR(temp.pivot_col[i], *ring); 1832 temp.pivot_col[i].idx = i; 1833 temp.s_mark[i] = 0; 1834 } 1835 } 1836 } 1837 1838 inline void CLUFactorRational::freeFactorRings(void) 1839 { 1840 1841 if(temp.pivot_col) 1842 spx_free(temp.pivot_col); 1843 1844 if(temp.pivot_colNZ) 1845 spx_free(temp.pivot_colNZ); 1846 1847 if(temp.pivot_row) 1848 spx_free(temp.pivot_row); 1849 1850 if(temp.pivot_rowNZ) 1851 spx_free(temp.pivot_rowNZ); 1852 } 1853 1854 1855 /* 1856 * Eliminate all row singletons from nucleus. 1857 * A row singleton may well be column singleton at the same time! 1858 */ 1859 inline void CLUFactorRational::eliminateRowSingletons() 1860 { 1861 int i, j, k, ll, r; 1862 int len, lk; 1863 int pcol, prow; 1864 Rational pval; 1865 int* idx; 1866 CLUFactorRational::Pring* sing; 1867 1868 for(sing = temp.pivot_rowNZ[1].prev; sing != &(temp.pivot_rowNZ[1]); sing = sing->prev) 1869 { 1870 prow = sing->idx; 1871 i = u.row.start[prow]; 1872 pcol = u.row.idx[i]; 1873 pval = u.row.val[i]; 1874 setPivot(temp.stage++, pcol, prow, pval); 1875 u.row.len[prow] = 0; 1876 removeDR(temp.pivot_col[pcol]); 1877 1878 /* Eliminate pivot column and build L vector. 1879 */ 1880 i = temp.s_cact[pcol]; 1881 1882 if(i > 1) 1883 { 1884 idx = &(u.col.idx[u.col.start[pcol]]); 1885 len = u.col.len[pcol]; 1886 lk = makeLvec(i - 1, prow); 1887 i = u.col.len[pcol] -= i; 1888 1889 for(; (r = idx[i]) != prow; ++i) 1890 { 1891 /* Find pivot column in row. 1892 */ 1893 ll = --(u.row.len[r]); 1894 k = u.row.start[r] + ll; 1895 1896 for(j = k; u.row.idx[j] != pcol; --j) 1897 ; 1898 1899 assert(j >= u.row.start[r]); 1900 1901 /* Initialize L vector 1902 */ 1903 l.idx[lk] = r; 1904 1905 l.val[lk] = u.row.val[j] / pval; 1906 1907 ++lk; 1908 1909 /* Remove pivot column from row. 1910 */ 1911 u.row.idx[j] = u.row.idx[k]; 1912 1913 u.row.val[j] = u.row.val[k]; 1914 1915 /* Move column to appropriate nonzero ring. 1916 */ 1917 removeDR(temp.pivot_row[r]); 1918 1919 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]); 1920 1921 assert(row.perm[r] < 0); 1922 1923 temp.s_max[r] = -1; 1924 } 1925 1926 /* skip pivot element */ 1927 assert(i < len && "ERROR: pivot column does not contain pivot row"); 1928 1929 for(++i; i < len; ++i) 1930 { 1931 /* Find pivot column in row. 1932 */ 1933 r = idx[i]; 1934 ll = --(u.row.len[r]); 1935 k = u.row.start[r] + ll; 1936 1937 for(j = k; u.row.idx[j] != pcol; --j) 1938 ; 1939 1940 assert(j >= u.row.start[r]); 1941 1942 /* Initialize L vector 1943 */ 1944 l.idx[lk] = r; 1945 1946 l.val[lk] = u.row.val[j] / pval; 1947 1948 ++lk; 1949 1950 /* Remove pivot column from row. 1951 */ 1952 u.row.idx[j] = u.row.idx[k]; 1953 1954 u.row.val[j] = u.row.val[k]; 1955 1956 /* Move column to appropriate nonzero ring. 1957 */ 1958 removeDR(temp.pivot_row[r]); 1959 1960 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[ll]); 1961 1962 assert(row.perm[r] < 0); 1963 1964 temp.s_max[r] = -1; 1965 } 1966 } 1967 else 1968 u.col.len[pcol] -= i; 1969 } 1970 1971 initDR(temp.pivot_rowNZ[1]); /* Remove all row singletons from list */ 1972 } 1973 1974 1975 1976 /* 1977 * Eliminate all column singletons from nucleus. 1978 * A column singleton must not be row singleton at the same time! 1979 */ 1980 inline void CLUFactorRational::eliminateColSingletons() 1981 { 1982 int i, j, k, m, c; 1983 int pcol, prow; 1984 CLUFactorRational::Pring* sing; 1985 1986 for(sing = temp.pivot_colNZ[1].prev; 1987 sing != &(temp.pivot_colNZ[1]); 1988 sing = sing->prev) 1989 { 1990 /* Find pivot value 1991 */ 1992 pcol = sing->idx; 1993 j = --(u.col.len[pcol]) + u.col.start[pcol]; /* remove pivot column */ 1994 prow = u.col.idx[j]; 1995 removeDR(temp.pivot_row[prow]); 1996 1997 j = --(u.row.len[prow]) + u.row.start[prow]; 1998 1999 for(i = j; (c = u.row.idx[i]) != pcol; --i) 2000 { 2001 m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--; 2002 2003 for(k = m; u.col.idx[k] != prow; ++k) 2004 ; 2005 2006 u.col.idx[k] = u.col.idx[m]; 2007 2008 u.col.idx[m] = prow; 2009 2010 m = temp.s_cact[c]; 2011 2012 removeDR(temp.pivot_col[c]); 2013 2014 init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]); 2015 2016 assert(col.perm[c] < 0); 2017 } 2018 2019 /* remove pivot element from pivot row 2020 */ 2021 setPivot(temp.stage++, pcol, prow, u.row.val[i]); 2022 2023 u.row.idx[i] = u.row.idx[j]; 2024 2025 u.row.val[i] = u.row.val[j]; 2026 2027 j = u.row.start[prow]; 2028 2029 for(--i; i >= j; --i) 2030 { 2031 c = u.row.idx[i]; 2032 m = u.col.len[c] + u.col.start[c] - (temp.s_cact[c])--; 2033 2034 for(k = m; u.col.idx[k] != prow; ++k) 2035 ; 2036 2037 u.col.idx[k] = u.col.idx[m]; 2038 2039 u.col.idx[m] = prow; 2040 2041 m = temp.s_cact[c]; 2042 2043 removeDR(temp.pivot_col[c]); 2044 2045 init2DR(temp.pivot_col[c], temp.pivot_colNZ[m]); 2046 2047 assert(col.perm[c] < 0); 2048 } 2049 } 2050 2051 initDR(temp.pivot_colNZ[1]); /* Remove all column singletons from list */ 2052 } 2053 2054 /* 2055 * No singletons available: Select pivot elements. 2056 */ 2057 inline void CLUFactorRational::selectPivots(const Rational& threshold) 2058 { 2059 int ii; 2060 int i; 2061 int j; 2062 int k; 2063 int ll = -1; // This value should never be used. 2064 int kk; 2065 int m; 2066 int count; 2067 int num; 2068 int rw = -1; // This value should never be used. 2069 int cl = -1; // This value should never be used. 2070 int len; 2071 int beg; 2072 Rational l_maxabs; 2073 Rational x = 0; // This value should never be used. 2074 int mkwtz; 2075 int candidates; 2076 2077 candidates = thedim - temp.stage - 1; 2078 2079 if(candidates > 4) 2080 candidates = 4; 2081 2082 num = 0; 2083 2084 count = 2; 2085 2086 for(;;) 2087 { 2088 ii = -1; 2089 2090 if(temp.pivot_rowNZ[count].next != &(temp.pivot_rowNZ[count])) 2091 { 2092 rw = temp.pivot_rowNZ[count].next->idx; 2093 beg = u.row.start[rw]; 2094 len = u.row.len[rw] + beg - 1; 2095 2096 /* set l_maxabs to maximum absolute value in row 2097 * (compute it if necessary). 2098 */ 2099 2100 if(temp.s_max[rw] < 0) 2101 { 2102 l_maxabs = spxAbs(u.row.val[len]); 2103 2104 for(i = len - 1; i >= beg; --i) 2105 if(l_maxabs < spxAbs(u.row.val[i])) 2106 l_maxabs = spxAbs(u.row.val[i]); 2107 2108 temp.s_max[rw] = l_maxabs; /* ##### */ 2109 } 2110 else 2111 l_maxabs = temp.s_max[rw]; 2112 2113 l_maxabs *= threshold; 2114 2115 /* select pivot element with lowest markowitz number in row 2116 */ 2117 mkwtz = thedim + 1; 2118 2119 for(i = len; i >= beg; --i) 2120 { 2121 k = u.row.idx[i]; 2122 j = temp.s_cact[k]; 2123 x = u.row.val[i]; 2124 2125 if(j < mkwtz && spxAbs(x) > l_maxabs) 2126 { 2127 mkwtz = j; 2128 cl = k; 2129 ii = i; 2130 2131 if(j <= count) /* ##### */ 2132 break; 2133 } 2134 } 2135 } 2136 else if(temp.pivot_colNZ[count].next != &(temp.pivot_colNZ[count])) 2137 { 2138 cl = temp.pivot_colNZ[count].next->idx; 2139 beg = u.col.start[cl]; 2140 len = u.col.len[cl] + beg - 1; 2141 beg = len - temp.s_cact[cl] + 1; 2142 assert(count == temp.s_cact[cl]); 2143 2144 /* select pivot element with lowest markowitz number in column 2145 */ 2146 mkwtz = thedim + 1; 2147 2148 for(i = len; i >= beg; --i) 2149 { 2150 k = u.col.idx[i]; 2151 j = u.row.len[k]; 2152 2153 if(j < mkwtz) 2154 { 2155 /* ensure that element (cl,k) is stable. 2156 */ 2157 if(temp.s_max[k] > 0) 2158 { 2159 /* case 1: l_maxabs is known 2160 */ 2161 for(m = u.row.start[k], kk = m + u.row.len[k] - 1; 2162 kk >= m; --kk) 2163 { 2164 if(u.row.idx[kk] == cl) 2165 { 2166 x = u.row.val[kk]; 2167 ll = kk; 2168 break; 2169 } 2170 } 2171 2172 l_maxabs = temp.s_max[k]; 2173 } 2174 else 2175 { 2176 /* case 2: l_maxabs needs to be computed 2177 */ 2178 m = u.row.start[k]; 2179 l_maxabs = spxAbs(u.row.val[m]); 2180 2181 for(kk = m + u.row.len[k] - 1; kk >= m; --kk) 2182 { 2183 if(l_maxabs < spxAbs(u.row.val[kk])) 2184 l_maxabs = spxAbs(u.row.val[kk]); 2185 2186 if(u.row.idx[kk] == cl) 2187 { 2188 x = u.row.val[kk]; 2189 ll = kk; 2190 break; 2191 } 2192 } 2193 2194 for(--kk; kk > m; --kk) 2195 { 2196 if(l_maxabs < spxAbs(u.row.val[kk])) 2197 l_maxabs = spxAbs(u.row.val[kk]); 2198 } 2199 2200 temp.s_max[k] = l_maxabs; 2201 } 2202 2203 l_maxabs *= threshold; 2204 2205 if(spxAbs(x) > l_maxabs) 2206 { 2207 mkwtz = j; 2208 rw = k; 2209 ii = ll; 2210 2211 if(j <= count + 1) 2212 break; 2213 } 2214 } 2215 } 2216 } 2217 else 2218 { 2219 ++count; 2220 continue; 2221 } 2222 2223 assert(cl >= 0); 2224 2225 removeDR(temp.pivot_col[cl]); 2226 initDR(temp.pivot_col[cl]); 2227 2228 if(ii >= 0) 2229 { 2230 /* Initialize selected pivot element 2231 */ 2232 CLUFactorRational::Pring* pr; 2233 temp.pivot_row[rw].pos = ii - u.row.start[rw]; 2234 temp.pivot_row[rw].mkwtz = mkwtz = (mkwtz - 1) * (count - 1); 2235 // ??? mkwtz originally was long, 2236 // maybe to avoid an overflow in this instruction? 2237 2238 for(pr = temp.pivots.next; pr->idx >= 0; pr = pr->next) 2239 { 2240 if(pr->idx == rw || pr->mkwtz >= mkwtz) 2241 break; 2242 } 2243 2244 pr = pr->prev; 2245 2246 if(pr->idx != rw) 2247 { 2248 removeDR(temp.pivot_row[rw]); 2249 init2DR(temp.pivot_row[rw], *pr); 2250 } 2251 2252 num++; 2253 2254 if(num >= candidates) 2255 break; 2256 } 2257 } 2258 2259 /* 2260 * while(temp.temp.next->mkwtz < temp.temp.prev->mkwtz) 2261 * { 2262 * Pring *pr; 2263 * pr = temp.temp.prev; 2264 * removeDR(*pr); 2265 * init2DR (*pr, rowNZ[u.row.len[pr->idx]]); 2266 } 2267 */ 2268 2269 assert(row.perm[rw] < 0); 2270 2271 assert(col.perm[cl] < 0); 2272 } 2273 2274 2275 /* 2276 * Perform L and update loop for row r 2277 */ 2278 inline int CLUFactorRational::updateRow(int r, 2279 int lv, 2280 int prow, 2281 int pcol, 2282 const Rational& pval) 2283 { 2284 int fill; 2285 Rational x, lx; 2286 int c, i, j, k, ll, m, n; 2287 2288 n = u.row.start[r]; 2289 m = --(u.row.len[r]) + n; 2290 2291 /* compute L vector entry and 2292 * and remove pivot column form row file 2293 */ 2294 2295 for(j = m; u.row.idx[j] != pcol; --j) 2296 ; 2297 2298 lx = u.row.val[j] / pval; 2299 2300 l.val[lv] = lx; 2301 2302 l.idx[lv] = r; 2303 2304 ++lv; 2305 2306 u.row.idx[j] = u.row.idx[m]; 2307 2308 u.row.val[j] = u.row.val[m]; 2309 2310 2311 /* update loop (I) and 2312 * computing expected fill 2313 */ 2314 fill = u.row.len[prow]; 2315 2316 for(j = m - 1; j >= n; --j) 2317 { 2318 c = u.row.idx[j]; 2319 2320 if(temp.s_mark[c]) 2321 { 2322 /* count fill elements. 2323 */ 2324 temp.s_mark[c] = 0; 2325 --fill; 2326 2327 /* update row values 2328 */ 2329 x = u.row.val[j] -= work[c] * lx; 2330 2331 if(x == 0) 2332 { 2333 /* Eliminate zero from row r 2334 */ 2335 --u.row.len[r]; 2336 --m; 2337 u.row.val[j] = u.row.val[m]; 2338 u.row.idx[j] = u.row.idx[m]; 2339 2340 /* Eliminate zero from column c 2341 */ 2342 --(temp.s_cact[c]); 2343 k = --(u.col.len[c]) + u.col.start[c]; 2344 2345 for(i = k; u.col.idx[i] != r; --i) 2346 ; 2347 2348 u.col.idx[i] = u.col.idx[k]; 2349 } 2350 } 2351 } 2352 2353 2354 /* create space for fill in row file 2355 */ 2356 ll = u.row.len[r]; 2357 2358 if(ll + fill > u.row.max[r]) 2359 remaxRow(r, ll + fill); 2360 2361 ll += u.row.start[r]; 2362 2363 /* fill creating update loop (II) 2364 */ 2365 for(j = u.row.start[prow], m = j + u.row.len[prow]; j < m; ++j) 2366 { 2367 c = u.row.idx[j]; 2368 2369 if(temp.s_mark[c]) 2370 { 2371 x = - work[c] * lx; 2372 2373 if(x != 0) 2374 { 2375 /* produce fill element in row r 2376 */ 2377 u.row.val[ll] = x; 2378 u.row.idx[ll] = c; 2379 ll++; 2380 u.row.len[r]++; 2381 2382 /* produce fill element in column c 2383 */ 2384 2385 if(u.col.len[c] >= u.col.max[c]) 2386 remaxCol(c, u.col.len[c] + 1); 2387 2388 u.col.idx[u.col.start[c] + (u.col.len[c])++] = r; 2389 2390 temp.s_cact[c]++; 2391 } 2392 } 2393 else 2394 temp.s_mark[c] = 1; 2395 } 2396 2397 /* move row to appropriate list. 2398 */ 2399 removeDR(temp.pivot_row[r]); 2400 2401 init2DR(temp.pivot_row[r], temp.pivot_rowNZ[u.row.len[r]]); 2402 2403 assert(row.perm[r] < 0); 2404 2405 temp.s_max[r] = -1; 2406 2407 return lv; 2408 } 2409 2410 /* 2411 * Eliminate pivot element 2412 */ 2413 inline void CLUFactorRational::eliminatePivot(int prow, int pos) 2414 { 2415 int i, j, k, m = -1; 2416 int lv = -1; // This value should never be used. 2417 int pcol; 2418 Rational pval; 2419 int pbeg = u.row.start[prow]; 2420 int plen = --(u.row.len[prow]); 2421 int pend = pbeg + plen; 2422 2423 2424 /* extract pivot element */ 2425 i = pbeg + pos; 2426 pcol = u.row.idx[i]; 2427 pval = u.row.val[i]; 2428 removeDR(temp.pivot_col[pcol]); 2429 initDR(temp.pivot_col[pcol]); 2430 2431 /* remove pivot from pivot row */ 2432 u.row.idx[i] = u.row.idx[pend]; 2433 u.row.val[i] = u.row.val[pend]; 2434 2435 /* set pivot element and construct L vector */ 2436 setPivot(temp.stage++, pcol, prow, pval); 2437 2438 /**@todo If this test failes, lv has no value. I suppose that in this 2439 * case none of the loops below that uses lv is executed. 2440 * But this is unproven. 2441 */ 2442 2443 if(temp.s_cact[pcol] - 1 > 0) 2444 lv = makeLvec(temp.s_cact[pcol] - 1, prow); 2445 2446 /* init working vector, 2447 * remove pivot row from working matrix 2448 * and remove columns from list. 2449 */ 2450 for(i = pbeg; i < pend; ++i) 2451 { 2452 j = u.row.idx[i]; 2453 temp.s_mark[j] = 1; 2454 work[j] = u.row.val[i]; 2455 removeDR(temp.pivot_col[j]); 2456 m = u.col.start[j] + u.col.len[j] - temp.s_cact[j]; 2457 2458 for(k = m; u.col.idx[k] != prow; ++k) 2459 ; 2460 2461 u.col.idx[k] = u.col.idx[m]; 2462 2463 u.col.idx[m] = prow; 2464 2465 temp.s_cact[j]--; 2466 } 2467 2468 /* perform L and update loop 2469 */ 2470 for(i = u.col.len[pcol] - temp.s_cact[pcol]; 2471 (m = u.col.idx[u.col.start[pcol] + i]) != prow; 2472 ++i) 2473 { 2474 assert(row.perm[m] < 0); 2475 assert(lv >= 0); 2476 /* coverity[negative_returns] */ 2477 updateRow(m, lv++, prow, pcol, pval); 2478 } 2479 2480 /* skip pivot row */ 2481 2482 m = u.col.len[pcol]; 2483 2484 for(++i; i < m; ++i) 2485 { 2486 assert(lv >= 0); 2487 /* coverity[negative_returns] */ 2488 updateRow(u.col.idx[u.col.start[pcol] + i], lv++, prow, pcol, pval); 2489 } 2490 2491 /* remove pivot column from column file. 2492 */ 2493 u.col.len[pcol] -= temp.s_cact[pcol]; 2494 2495 /* clear working vector and reinsert columns to lists 2496 */ 2497 for(i = u.row.start[prow], pend = i + plen; i < pend; ++i) 2498 { 2499 j = u.row.idx[i]; 2500 work[j] = 0; 2501 temp.s_mark[j] = 0; 2502 init2DR(temp.pivot_col[j], temp.pivot_colNZ[temp.s_cact[j]]); 2503 assert(col.perm[j] < 0); 2504 } 2505 } 2506 2507 2508 /* 2509 * Factorize nucleus. 2510 */ 2511 inline void CLUFactorRational::eliminateNucleus(const Rational& threshold) 2512 { 2513 int r, c; 2514 CLUFactorRational::Pring* pivot; 2515 2516 if(stat == SLinSolverRational::SINGULAR) 2517 return; 2518 2519 temp.pivots.mkwtz = -1; 2520 2521 temp.pivots.idx = -1; 2522 2523 temp.pivots.pos = -1; 2524 2525 while(temp.stage < thedim - 1) 2526 { 2527 #ifndef NDEBUG 2528 int i; 2529 // CLUFactorRationalIsConsistent(fac); 2530 2531 for(i = 0; i < thedim; ++i) 2532 if(col.perm[i] < 0) 2533 assert(temp.s_mark[i] == 0); 2534 2535 #endif 2536 2537 if(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1])) 2538 { 2539 if(timeLimitReached()) 2540 return; 2541 2542 /* row singleton available */ 2543 eliminateRowSingletons(); 2544 } 2545 else if(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1])) 2546 { 2547 if(timeLimitReached()) 2548 return; 2549 2550 /* column singleton available */ 2551 eliminateColSingletons(); 2552 } 2553 else 2554 { 2555 initDR(temp.pivots); 2556 selectPivots(threshold); 2557 2558 assert(temp.pivots.next != &temp.pivots && 2559 "ERROR: no pivot element selected"); 2560 2561 for(pivot = temp.pivots.next; pivot != &temp.pivots; 2562 pivot = pivot->next) 2563 { 2564 if(timeLimitReached()) 2565 return; 2566 2567 eliminatePivot(pivot->idx, pivot->pos); 2568 } 2569 } 2570 2571 if(temp.pivot_rowNZ->next != temp.pivot_rowNZ || 2572 temp.pivot_colNZ->next != temp.pivot_colNZ) 2573 { 2574 stat = SLinSolverRational::SINGULAR; 2575 return; 2576 } 2577 } 2578 2579 if(temp.stage < thedim) 2580 { 2581 /* Eliminate remaining element. 2582 * Note, that this must be both, column and row singleton. 2583 */ 2584 assert(temp.pivot_rowNZ[1].next != &(temp.pivot_rowNZ[1]) && 2585 "ERROR: one row must be left"); 2586 assert(temp.pivot_colNZ[1].next != &(temp.pivot_colNZ[1]) && 2587 "ERROR: one col must be left"); 2588 r = temp.pivot_rowNZ[1].next->idx; 2589 c = temp.pivot_colNZ[1].next->idx; 2590 u.row.len[r] = 0; 2591 u.col.len[c]--; 2592 setPivot(temp.stage, c, r, u.row.val[u.row.start[r]]); 2593 } 2594 } 2595 2596 /*****************************************************************************/ 2597 2598 inline int CLUFactorRational::setupColVals() 2599 { 2600 int i; 2601 int n = thedim; 2602 2603 u.col.val.reDim(u.col.size); 2604 2605 for(i = 0; i < thedim; i++) 2606 u.col.len[i] = 0; 2607 2608 maxabs = 0; 2609 2610 for(i = 0; i < thedim; i++) 2611 { 2612 int k = u.row.start[i]; 2613 int* idx = &u.row.idx[k]; 2614 Rational* val = &u.row.val[k]; 2615 int len = u.row.len[i]; 2616 2617 n += len; 2618 2619 while(len-- > 0) 2620 { 2621 assert((*idx >= 0) && (*idx < thedim)); 2622 2623 k = u.col.start[*idx] + u.col.len[*idx]; 2624 2625 assert((k >= 0) && (k < u.col.size)); 2626 2627 u.col.len[*idx]++; 2628 2629 assert(u.col.len[*idx] <= u.col.max[*idx]); 2630 2631 u.col.idx[k] = i; 2632 u.col.val[k] = *val; 2633 2634 if(spxAbs(*val) > maxabs) 2635 maxabs = spxAbs(*val); 2636 2637 idx++; 2638 2639 val++; 2640 } 2641 } 2642 2643 return n; 2644 } 2645 2646 /*****************************************************************************/ 2647 2648 #ifdef SOPLEX_WITH_L_ROWS 2649 inline void CLUFactorRational::setupRowVals() 2650 { 2651 int i, j, k, m; 2652 int vecs, mem; 2653 int* l_row; 2654 int* idx; 2655 int* beg; 2656 int* l_ridx; 2657 int* l_rbeg; 2658 int* rorig; 2659 int* rrorig; 2660 int* rperm; 2661 int* rrperm; 2662 2663 vecs = l.firstUpdate; 2664 l_row = l.row; 2665 idx = l.idx; 2666 VectorRational& val = l.val; 2667 int validx = 0; 2668 beg = l.start; 2669 mem = beg[vecs]; 2670 2671 if(l.ridx) 2672 spx_free(l.ridx); 2673 2674 if(l.rbeg) 2675 spx_free(l.rbeg); 2676 2677 if(l.rorig) 2678 spx_free(l.rorig); 2679 2680 if(l.rperm) 2681 spx_free(l.rperm); 2682 2683 l.rval.reDim(mem + 1); 2684 2685 spx_alloc(l.ridx, mem); 2686 2687 spx_alloc(l.rbeg, thedim + 1); 2688 2689 spx_alloc(l.rorig, thedim); 2690 2691 spx_alloc(l.rperm, thedim); 2692 2693 l_ridx = l.ridx; 2694 2695 VectorRational& l_rval = l.rval; 2696 2697 l_rbeg = l.rbeg; 2698 2699 rorig = l.rorig; 2700 2701 rrorig = row.orig; 2702 2703 rperm = l.rperm; 2704 2705 rrperm = row.perm; 2706 2707 for(i = thedim; i--; *l_rbeg++ = 0) 2708 { 2709 *rorig++ = *rrorig++; 2710 *rperm++ = *rrperm++; 2711 } 2712 2713 *l_rbeg = 0; 2714 2715 l_rbeg = l.rbeg + 1; 2716 2717 for(i = mem; i--;) 2718 l_rbeg[*idx++]++; 2719 2720 idx = l.idx; 2721 2722 for(m = 0, i = thedim; i--; l_rbeg++) 2723 { 2724 j = *l_rbeg; 2725 *l_rbeg = m; 2726 m += j; 2727 } 2728 2729 assert(m == mem); 2730 2731 l_rbeg = l.rbeg + 1; 2732 2733 for(i = j = 0; i < vecs; ++i) 2734 { 2735 m = l_row[i]; 2736 assert(idx == &l.idx[l.start[i]]); 2737 2738 for(; j < beg[i + 1]; j++) 2739 { 2740 k = l_rbeg[*idx++]++; 2741 assert(k < mem); 2742 l_ridx[k] = m; 2743 l_rval[k] = val[validx]; 2744 validx++; // was l_rval[k] = *val++; with Rational* val 2745 } 2746 } 2747 2748 assert(l.rbeg[thedim] == mem); 2749 2750 assert(l.rbeg[0] == 0); 2751 } 2752 2753 #endif 2754 2755 /*****************************************************************************/ 2756 2757 inline void CLUFactorRational::factor( 2758 const SVectorRational** vec, ///< Array of column vector pointers 2759 const Rational& threshold ///< pivoting threshold 2760 ) 2761 { 2762 SPxOut::debug(this, "CLUFactorRational::factor()\n"); 2763 2764 factorTime->start(); 2765 2766 stat = SLinSolverRational::OK; 2767 2768 l.start[0] = 0; 2769 l.firstUpdate = 0; 2770 l.firstUnused = 0; 2771 2772 temp.init(thedim); 2773 initPerm(); 2774 2775 initFactorMatrix(vec); 2776 2777 if(stat) 2778 goto TERMINATE; 2779 2780 if(timeLimitReached()) 2781 goto TERMINATE; 2782 2783 // initMaxabs = initMaxabs; 2784 2785 colSingletons(); 2786 2787 if(stat != SLinSolverRational::OK) 2788 goto TERMINATE; 2789 2790 if(timeLimitReached()) 2791 goto TERMINATE; 2792 2793 rowSingletons(); 2794 2795 if(stat != SLinSolverRational::OK) 2796 goto TERMINATE; 2797 2798 if(temp.stage < thedim) 2799 { 2800 if(timeLimitReached()) 2801 goto TERMINATE; 2802 2803 initFactorRings(); 2804 eliminateNucleus(threshold); 2805 freeFactorRings(); 2806 } 2807 2808 TERMINATE: 2809 2810 l.firstUpdate = l.firstUnused; 2811 2812 if(stat == SLinSolverRational::OK) 2813 { 2814 #ifdef SOPLEX_WITH_L_ROWS 2815 setupRowVals(); 2816 #endif 2817 nzCnt = setupColVals(); 2818 } 2819 2820 factorTime->stop(); 2821 2822 factorCount++; 2823 } 2824 2825 inline void CLUFactorRational::dump() const 2826 { 2827 int i, j, k; 2828 2829 // Dump regardless of the verbosity level if this method is called; 2830 2831 /* Dump U: 2832 */ 2833 2834 for(i = 0; i < thedim; ++i) 2835 { 2836 if(row.perm[i] >= 0) 2837 std::cout << "DCLUFA01 diag[" << i << "]: [" << col.orig[row.perm[i]] 2838 << "] = " << diag[i] << std::endl; 2839 2840 for(j = 0; j < u.row.len[i]; ++j) 2841 std::cout << "DCLUFA02 u[" << i << "]: [" 2842 << u.row.idx[u.row.start[i] + j] << "] = " 2843 << u.row.val[u.row.start[i] + j] << std::endl; 2844 } 2845 2846 /* Dump L: 2847 */ 2848 for(i = 0; i < thedim; ++i) 2849 { 2850 for(j = 0; j < l.firstUnused; ++j) 2851 if(col.orig[row.perm[l.row[j]]] == i) 2852 { 2853 std::cout << "DCLUFA03 l[" << i << "]" << std::endl; 2854 2855 for(k = l.start[j]; k < l.start[j + 1]; ++k) 2856 std::cout << "DCLUFA04 l[" << k - l.start[j] 2857 << "]: [" << l.idx[k] 2858 << "] = " << l.val[k] << std::endl; 2859 2860 break; 2861 } 2862 } 2863 2864 return; 2865 } 2866 2867 /*****************************************************************************/ 2868 /* 2869 * Ensure that row memory is at least size. 2870 */ 2871 inline void CLUFactorRational::minRowMem(int size) 2872 { 2873 if(u.row.val.dim() < size) 2874 { 2875 u.row.val.reDim(size); 2876 spx_realloc(u.row.idx, size); 2877 } 2878 2879 assert(u.row.val.dim() >= size); 2880 } 2881 2882 /*****************************************************************************/ 2883 /* 2884 * Ensure that column memory is at least size. 2885 */ 2886 inline void CLUFactorRational::minColMem(int size) 2887 { 2888 2889 if(u.col.size < size) 2890 { 2891 u.col.size = size; 2892 spx_realloc(u.col.idx, size); 2893 } 2894 } 2895 2896 inline void CLUFactorRational::forestMinColMem(int size) 2897 { 2898 2899 if(u.col.size < size) 2900 { 2901 u.col.size = size; 2902 spx_realloc(u.col.idx, size); 2903 u.col.val.reDim(size); 2904 } 2905 } 2906 2907 inline void CLUFactorRational::minLMem(int size) 2908 { 2909 2910 if(size > l.val.dim()) 2911 { 2912 int newsize = int(0.2 * l.val.dim() + size); 2913 l.val.reDim(newsize); 2914 spx_realloc(l.idx, l.val.dim()); 2915 } 2916 } 2917 2918 2919 inline int CLUFactorRational::makeLvec(int p_len, int p_row) 2920 { 2921 2922 if(l.firstUnused >= l.startSize) 2923 { 2924 l.startSize += 100; 2925 spx_realloc(l.start, l.startSize); 2926 } 2927 2928 int* p_lrow = l.row; 2929 2930 int* p_lbeg = l.start; 2931 int first = p_lbeg[l.firstUnused]; 2932 2933 assert(p_len > 0 && "ERROR: no empty columns allowed in L vectors"); 2934 2935 minLMem(first + p_len); 2936 p_lrow[l.firstUnused] = p_row; 2937 l.start[++(l.firstUnused)] = first + p_len; 2938 2939 assert(l.start[l.firstUnused] <= l.val.dim()); 2940 assert(l.firstUnused <= l.startSize); 2941 return first; 2942 } 2943 2944 2945 /*****************************************************************************/ 2946 2947 inline bool CLUFactorRational::isConsistent() const 2948 { 2949 #ifdef ENABLE_CONSISTENCY_CHECKS 2950 int i, j, k, ll; 2951 Dring* ring; 2952 CLUFactorRational::Pring* pring; 2953 2954 /* Consistency only relevant for real factorizations 2955 */ 2956 2957 if(stat) 2958 return true; 2959 2960 /* Test column ring list consistency. 2961 */ 2962 i = 0; 2963 2964 for(ring = u.col.list.next; ring != &(u.col.list); ring = ring->next) 2965 { 2966 assert(ring->idx >= 0); 2967 assert(ring->idx < thedim); 2968 assert(ring->prev->next == ring); 2969 2970 if(ring != u.col.list.next) 2971 { 2972 assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx] 2973 == u.col.start[ring->idx]); 2974 } 2975 2976 ++i; 2977 } 2978 2979 assert(i == thedim); 2980 2981 assert(u.col.start[ring->prev->idx] + u.col.max[ring->prev->idx] 2982 == u.col.used); 2983 2984 2985 /* Test row ring list consistency. 2986 */ 2987 i = 0; 2988 2989 for(ring = u.row.list.next; ring != &(u.row.list); ring = ring->next) 2990 { 2991 assert(ring->idx >= 0); 2992 assert(ring->idx < thedim); 2993 assert(ring->prev->next == ring); 2994 assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx] 2995 == u.row.start[ring->idx]); 2996 ++i; 2997 } 2998 2999 assert(i == thedim); 3000 3001 assert(u.row.start[ring->prev->idx] + u.row.max[ring->prev->idx] 3002 == u.row.used); 3003 3004 3005 /* Test consistency of individual svectors. 3006 */ 3007 3008 for(i = 0; i < thedim; ++i) 3009 { 3010 assert(u.row.max[i] >= u.row.len[i]); 3011 assert(u.col.max[i] >= u.col.len[i]); 3012 } 3013 3014 3015 /* Test consistency of column file to row file of U 3016 */ 3017 for(i = 0; i < thedim; ++i) 3018 { 3019 for(j = u.row.start[i] + u.row.len[i] - 1; j >= u.row.start[i]; j--) 3020 { 3021 k = u.row.idx[j]; 3022 3023 for(ll = u.col.start[k] + u.col.len[k] - 1; ll >= u.col.start[k]; ll--) 3024 { 3025 if(u.col.idx[ll] == i) 3026 break; 3027 } 3028 3029 assert(u.col.idx[ll] == i); 3030 3031 if(row.perm[i] < 0) 3032 { 3033 assert(col.perm[k] < 0); 3034 } 3035 else 3036 { 3037 assert(col.perm[k] < 0 || col.perm[k] > row.perm[i]); 3038 } 3039 } 3040 } 3041 3042 /* Test consistency of row file to column file of U 3043 */ 3044 for(i = 0; i < thedim; ++i) 3045 { 3046 for(j = u.col.start[i] + u.col.len[i] - 1; j >= u.col.start[i]; j--) 3047 { 3048 k = u.col.idx[j]; 3049 3050 for(ll = u.row.start[k] + u.row.len[k] - 1; ll >= u.row.start[k]; ll--) 3051 { 3052 if(u.row.idx[ll] == i) 3053 break; 3054 } 3055 3056 assert(u.row.idx[ll] == i); 3057 3058 assert(col.perm[i] < 0 || row.perm[k] < col.perm[i]); 3059 } 3060 } 3061 3062 /* Test consistency of nonzero count lists 3063 */ 3064 if(temp.pivot_colNZ && temp.pivot_rowNZ) 3065 { 3066 for(i = 0; i < thedim - temp.stage; ++i) 3067 { 3068 for(pring = temp.pivot_rowNZ[i].next; pring != &(temp.pivot_rowNZ[i]); pring = pring->next) 3069 { 3070 assert(row.perm[pring->idx] < 0); 3071 } 3072 3073 for(pring = temp.pivot_colNZ[i].next; pring != &(temp.pivot_colNZ[i]); pring = pring->next) 3074 { 3075 assert(col.perm[pring->idx] < 0); 3076 } 3077 } 3078 } 3079 3080 #endif // CONSISTENCY_CHECKS 3081 3082 return true; 3083 } 3084 3085 inline void CLUFactorRational::solveUright(Rational* wrk, Rational* vec) 3086 { 3087 3088 for(int i = thedim - 1; i >= 0; i--) 3089 { 3090 int r = row.orig[i]; 3091 int c = col.orig[i]; 3092 Rational x = wrk[c] = diag[r] * vec[r]; 3093 3094 vec[r] = 0; 3095 3096 if(x != 0) 3097 //if (isNotZero(x)) 3098 { 3099 if(timeLimitReached()) 3100 return; 3101 3102 for(int j = u.col.start[c]; j < u.col.start[c] + u.col.len[c]; j++) 3103 vec[u.col.idx[j]] -= x * u.col.val[j]; 3104 } 3105 } 3106 } 3107 3108 inline int CLUFactorRational::solveUrightEps(Rational* vec, int* nonz, Rational* rhs) 3109 { 3110 int i, j, r, c, n; 3111 int* rorig, *corig; 3112 int* cidx, *clen, *cbeg; 3113 Rational x; 3114 3115 int* idx; 3116 Rational* val; 3117 3118 rorig = row.orig; 3119 corig = col.orig; 3120 3121 cidx = u.col.idx; 3122 VectorRational& cval = u.col.val; 3123 clen = u.col.len; 3124 cbeg = u.col.start; 3125 3126 n = 0; 3127 3128 for(i = thedim - 1; i >= 0; --i) 3129 { 3130 r = rorig[i]; 3131 x = diag[r] * rhs[r]; 3132 3133 if(x != 0) 3134 { 3135 c = corig[i]; 3136 vec[c] = x; 3137 nonz[n++] = c; 3138 val = &cval[cbeg[c]]; 3139 idx = &cidx[cbeg[c]]; 3140 j = clen[c]; 3141 3142 while(j-- > 0) 3143 rhs[*idx++] -= x * (*val++); 3144 } 3145 } 3146 3147 return n; 3148 } 3149 3150 inline void CLUFactorRational::solveUright2(Rational* p_work1, Rational* vec1, Rational* p_work2, 3151 Rational* vec2) 3152 { 3153 int i, j, r, c; 3154 int* rorig, *corig; 3155 int* cidx, *clen, *cbeg; 3156 Rational x1, x2; 3157 3158 int* idx; 3159 Rational* val; 3160 3161 rorig = row.orig; 3162 corig = col.orig; 3163 3164 cidx = u.col.idx; 3165 VectorRational& cval = u.col.val; 3166 clen = u.col.len; 3167 cbeg = u.col.start; 3168 3169 for(i = thedim - 1; i >= 0; --i) 3170 { 3171 r = rorig[i]; 3172 c = corig[i]; 3173 p_work1[c] = x1 = diag[r] * vec1[r]; 3174 p_work2[c] = x2 = diag[r] * vec2[r]; 3175 vec1[r] = vec2[r] = 0; 3176 3177 if(x1 != 0 && x2 != 0) 3178 { 3179 val = &cval[cbeg[c]]; 3180 idx = &cidx[cbeg[c]]; 3181 j = clen[c]; 3182 3183 while(j-- > 0) 3184 { 3185 vec1[*idx] -= x1 * (*val); 3186 vec2[*idx++] -= x2 * (*val++); 3187 } 3188 } 3189 else if(x1 != 0) 3190 { 3191 val = &cval[cbeg[c]]; 3192 idx = &cidx[cbeg[c]]; 3193 j = clen[c]; 3194 3195 while(j-- > 0) 3196 vec1[*idx++] -= x1 * (*val++); 3197 } 3198 else if(x2 != 0) 3199 { 3200 val = &cval[cbeg[c]]; 3201 idx = &cidx[cbeg[c]]; 3202 j = clen[c]; 3203 3204 while(j-- > 0) 3205 vec2[*idx++] -= x2 * (*val++); 3206 } 3207 } 3208 } 3209 3210 inline int CLUFactorRational::solveUright2eps(Rational* p_work1, Rational* vec1, Rational* p_work2, 3211 Rational* vec2, int* nonz) 3212 { 3213 int i, j, r, c, n; 3214 int* rorig, *corig; 3215 int* cidx, *clen, *cbeg; 3216 bool notzero1, notzero2; 3217 Rational x1, x2; 3218 3219 int* idx; 3220 Rational* val; 3221 3222 rorig = row.orig; 3223 corig = col.orig; 3224 3225 cidx = u.col.idx; 3226 VectorRational& cval = u.col.val; 3227 clen = u.col.len; 3228 cbeg = u.col.start; 3229 3230 n = 0; 3231 3232 for(i = thedim - 1; i >= 0; --i) 3233 { 3234 c = corig[i]; 3235 r = rorig[i]; 3236 p_work1[c] = x1 = diag[r] * vec1[r]; 3237 p_work2[c] = x2 = diag[r] * vec2[r]; 3238 vec1[r] = vec2[r] = 0; 3239 notzero1 = x1 != 0 ? 1 : 0; 3240 notzero2 = x2 != 0 ? 1 : 0; 3241 3242 if(notzero1 && notzero2) 3243 { 3244 *nonz++ = c; 3245 n++; 3246 val = &cval[cbeg[c]]; 3247 idx = &cidx[cbeg[c]]; 3248 j = clen[c]; 3249 3250 while(j-- > 0) 3251 { 3252 vec1[*idx] -= x1 * (*val); 3253 vec2[*idx++] -= x2 * (*val++); 3254 } 3255 } 3256 else if(notzero1) 3257 { 3258 p_work2[c] = 0; 3259 *nonz++ = c; 3260 n++; 3261 val = &cval[cbeg[c]]; 3262 idx = &cidx[cbeg[c]]; 3263 j = clen[c]; 3264 3265 while(j-- > 0) 3266 vec1[*idx++] -= x1 * (*val++); 3267 } 3268 else if(notzero2) 3269 { 3270 p_work1[c] = 0; 3271 val = &cval[cbeg[c]]; 3272 idx = &cidx[cbeg[c]]; 3273 j = clen[c]; 3274 3275 while(j-- > 0) 3276 vec2[*idx++] -= x2 * (*val++); 3277 } 3278 else 3279 { 3280 p_work1[c] = 0; 3281 p_work2[c] = 0; 3282 } 3283 } 3284 3285 return n; 3286 } 3287 3288 inline void CLUFactorRational::solveLright(Rational* vec) 3289 { 3290 int i, j, k; 3291 int end; 3292 Rational x; 3293 Rational* val; 3294 int* lrow, *lidx, *idx; 3295 int* lbeg; 3296 3297 VectorRational& lval = l.val; 3298 lidx = l.idx; 3299 lrow = l.row; 3300 lbeg = l.start; 3301 3302 end = l.firstUpdate; 3303 3304 for(i = 0; i < end; ++i) 3305 { 3306 if((x = vec[lrow[i]]) != 0) 3307 { 3308 if(timeLimitReached()) 3309 return; 3310 3311 SPxOut::debug(this, "y{}={}\n", lrow[i], vec[lrow[i]]); 3312 3313 k = lbeg[i]; 3314 idx = &(lidx[k]); 3315 val = &(lval[k]); 3316 3317 for(j = lbeg[i + 1]; j > k; --j) 3318 { 3319 SPxOut::debug(this, " -> y{} -= {} * {} = {} -> {}\n", *idx, x, *val, 3320 x * (*val), vec[*idx] - x * (*val)); 3321 vec[*idx++] -= x * (*val++); 3322 } 3323 } 3324 } 3325 3326 if(l.updateType) /* Forest-Tomlin Updates */ 3327 { 3328 SPxOut::debug(this, "performing FT updates...\n"); 3329 3330 end = l.firstUnused; 3331 3332 for(; i < end; ++i) 3333 { 3334 x = 0; 3335 k = lbeg[i]; 3336 idx = &(lidx[k]); 3337 val = &(lval[k]); 3338 3339 for(j = lbeg[i + 1]; j > k; --j) 3340 x += vec[*idx++] * (*val++); 3341 3342 vec[lrow[i]] -= x; 3343 3344 SPxOut::debug(this, "y{}={}\n", lrow[i], vec[lrow[i]]); 3345 } 3346 3347 SPxOut::debug(this, "finished FT updates.\n"); 3348 } 3349 } 3350 3351 inline void CLUFactorRational::solveLright2(Rational* vec1, Rational* vec2) 3352 { 3353 int i, j, k; 3354 int end; 3355 Rational x2; 3356 Rational x1; 3357 Rational* val; 3358 int* lrow, *lidx, *idx; 3359 int* lbeg; 3360 3361 VectorRational& lval = l.val; 3362 lidx = l.idx; 3363 lrow = l.row; 3364 lbeg = l.start; 3365 3366 end = l.firstUpdate; 3367 3368 for(i = 0; i < end; ++i) 3369 { 3370 x1 = vec1[lrow[i]]; 3371 x2 = vec2[lrow[i]]; 3372 3373 if(x1 != 0 && x2 != 0) 3374 { 3375 k = lbeg[i]; 3376 idx = &(lidx[k]); 3377 val = &(lval[k]); 3378 3379 for(j = lbeg[i + 1]; j > k; --j) 3380 { 3381 vec1[*idx] -= x1 * (*val); 3382 vec2[*idx++] -= x2 * (*val++); 3383 } 3384 } 3385 else if(x1 != 0) 3386 { 3387 k = lbeg[i]; 3388 idx = &(lidx[k]); 3389 val = &(lval[k]); 3390 3391 for(j = lbeg[i + 1]; j > k; --j) 3392 vec1[*idx++] -= x1 * (*val++); 3393 } 3394 else if(x2 != 0) 3395 { 3396 k = lbeg[i]; 3397 idx = &(lidx[k]); 3398 val = &(lval[k]); 3399 3400 for(j = lbeg[i + 1]; j > k; --j) 3401 vec2[*idx++] -= x2 * (*val++); 3402 } 3403 } 3404 3405 if(l.updateType) /* Forest-Tomlin Updates */ 3406 { 3407 end = l.firstUnused; 3408 3409 for(; i < end; ++i) 3410 { 3411 x1 = 0; 3412 x2 = 0; 3413 k = lbeg[i]; 3414 idx = &(lidx[k]); 3415 val = &(lval[k]); 3416 3417 for(j = lbeg[i + 1]; j > k; --j) 3418 { 3419 x1 += vec1[*idx] * (*val); 3420 x2 += vec2[*idx++] * (*val++); 3421 } 3422 3423 vec1[lrow[i]] -= x1; 3424 3425 vec2[lrow[i]] -= x2; 3426 } 3427 } 3428 } 3429 3430 inline void CLUFactorRational::solveUpdateRight(Rational* vec) 3431 { 3432 int i, j, k; 3433 int end; 3434 Rational x; 3435 Rational* val; 3436 int* lrow, *lidx, *idx; 3437 int* lbeg; 3438 3439 assert(!l.updateType); /* no Forest-Tomlin Updates */ 3440 3441 VectorRational& lval = l.val; 3442 lidx = l.idx; 3443 lrow = l.row; 3444 lbeg = l.start; 3445 3446 end = l.firstUnused; 3447 3448 for(i = l.firstUpdate; i < end; ++i) 3449 { 3450 if((x = vec[lrow[i]]) != 0) 3451 { 3452 k = lbeg[i]; 3453 idx = &(lidx[k]); 3454 val = &(lval[k]); 3455 3456 for(j = lbeg[i + 1]; j > k; --j) 3457 vec[*idx++] -= x * (*val++); 3458 } 3459 } 3460 } 3461 3462 inline void CLUFactorRational::solveUpdateRight2(Rational* vec1, Rational* vec2) 3463 { 3464 int i, j, k; 3465 int end; 3466 Rational x1, x2; 3467 int* lrow, *lidx; 3468 int* lbeg; 3469 3470 int* idx; 3471 Rational* val; 3472 3473 assert(!l.updateType); /* no Forest-Tomlin Updates */ 3474 3475 VectorRational& lval = l.val; 3476 lidx = l.idx; 3477 lrow = l.row; 3478 lbeg = l.start; 3479 3480 end = l.firstUnused; 3481 3482 for(i = l.firstUpdate; i < end; ++i) 3483 { 3484 x1 = vec1[lrow[i]]; 3485 x2 = vec2[lrow[i]]; 3486 3487 if(x1 != 0 && x2 != 0) 3488 { 3489 k = lbeg[i]; 3490 idx = &(lidx[k]); 3491 val = &(lval[k]); 3492 3493 for(j = lbeg[i + 1]; j > k; --j) 3494 { 3495 vec1[*idx] -= x1 * (*val); 3496 vec2[*idx++] -= x2 * (*val++); 3497 } 3498 } 3499 else if(x1 != 0) 3500 { 3501 k = lbeg[i]; 3502 idx = &(lidx[k]); 3503 val = &(lval[k]); 3504 3505 for(j = lbeg[i + 1]; j > k; --j) 3506 vec1[*idx++] -= x1 * (*val++); 3507 } 3508 else if(x2 != 0) 3509 { 3510 k = lbeg[i]; 3511 idx = &(lidx[k]); 3512 val = &(lval[k]); 3513 3514 for(j = lbeg[i + 1]; j > k; --j) 3515 vec2[*idx++] -= x2 * (*val++); 3516 } 3517 } 3518 } 3519 3520 inline int CLUFactorRational::solveRight4update(Rational* vec, int* nonz, 3521 Rational* rhs, Rational* forest, int* forestNum, int* forestIdx) 3522 { 3523 solveLright(rhs); 3524 3525 if(forest) 3526 { 3527 int n = 0; 3528 3529 for(int i = 0; i < thedim; i++) 3530 { 3531 forestIdx[n] = i; 3532 forest[i] = rhs[i]; 3533 n += rhs[i] != 0 ? 1 : 0; 3534 } 3535 3536 *forestNum = n; 3537 } 3538 3539 if(!l.updateType) /* no Forest-Tomlin Updates */ 3540 { 3541 solveUright(vec, rhs); 3542 solveUpdateRight(vec); 3543 return 0; 3544 } 3545 else 3546 return solveUrightEps(vec, nonz, rhs); 3547 } 3548 3549 inline void CLUFactorRational::solveRight(Rational* vec, Rational* rhs) 3550 { 3551 solveLright(rhs); 3552 solveUright(vec, rhs); 3553 3554 if(!l.updateType) /* no Forest-Tomlin Updates */ 3555 solveUpdateRight(vec); 3556 } 3557 3558 inline int CLUFactorRational::solveRight2update(Rational* vec1, 3559 Rational* vec2, 3560 Rational* rhs1, 3561 Rational* rhs2, 3562 int* nonz, 3563 Rational* forest, 3564 int* forestNum, 3565 int* forestIdx) 3566 { 3567 solveLright2(rhs1, rhs2); 3568 3569 if(forest) 3570 { 3571 int n = 0; 3572 3573 for(int i = 0; i < thedim; i++) 3574 { 3575 forestIdx[n] = i; 3576 forest[i] = rhs1[i]; 3577 n += rhs1[i] != 0 ? 1 : 0; 3578 } 3579 3580 *forestNum = n; 3581 } 3582 3583 if(!l.updateType) /* no Forest-Tomlin Updates */ 3584 { 3585 solveUright2(vec1, rhs1, vec2, rhs2); 3586 solveUpdateRight2(vec1, vec2); 3587 return 0; 3588 } 3589 else 3590 return solveUright2eps(vec1, rhs1, vec2, rhs2, nonz); 3591 } 3592 3593 inline void CLUFactorRational::solveRight2( 3594 Rational* vec1, 3595 Rational* vec2, 3596 Rational* rhs1, 3597 Rational* rhs2) 3598 { 3599 solveLright2(rhs1, rhs2); 3600 3601 if(l.updateType) /* Forest-Tomlin Updates */ 3602 solveUright2(vec1, rhs1, vec2, rhs2); 3603 else 3604 { 3605 solveUright2(vec1, rhs1, vec2, rhs2); 3606 solveUpdateRight2(vec1, vec2); 3607 } 3608 } 3609 3610 /*****************************************************************************/ 3611 inline void CLUFactorRational::solveUleft(Rational* p_work, Rational* vec) 3612 { 3613 for(int i = 0; i < thedim; ++i) 3614 { 3615 int c = col.orig[i]; 3616 int r = row.orig[i]; 3617 3618 assert(c >= 0); // Inna/Tobi: otherwise, vec[c] would be strange... 3619 assert(r >= 0); // Inna/Tobi: otherwise, diag[r] would be strange... 3620 3621 Rational x = vec[c]; 3622 3623 vec[c] = 0; 3624 3625 if(x != 0) 3626 { 3627 if(timeLimitReached()) 3628 return; 3629 3630 x *= diag[r]; 3631 p_work[r] = x; 3632 3633 int end = u.row.start[r] + u.row.len[r]; 3634 3635 for(int m = u.row.start[r]; m < end; m++) 3636 { 3637 vec[u.row.idx[m]] -= x * u.row.val[m]; 3638 } 3639 } 3640 } 3641 } 3642 3643 3644 3645 inline void CLUFactorRational::solveUleft2(Rational* p_work1, Rational* vec1, Rational* p_work2, 3646 Rational* vec2) 3647 { 3648 Rational x1; 3649 Rational x2; 3650 int i, k, r, c; 3651 int* rorig, *corig; 3652 int* ridx, *rlen, *rbeg, *idx; 3653 Rational* val; 3654 3655 rorig = row.orig; 3656 corig = col.orig; 3657 3658 ridx = u.row.idx; 3659 VectorRational& rval = u.row.val; 3660 rlen = u.row.len; 3661 rbeg = u.row.start; 3662 3663 for(i = 0; i < thedim; ++i) 3664 { 3665 c = corig[i]; 3666 r = rorig[i]; 3667 x1 = vec1[c]; 3668 x2 = vec2[c]; 3669 3670 if((x1 != 0) && (x2 != 0)) 3671 { 3672 x1 *= diag[r]; 3673 x2 *= diag[r]; 3674 p_work1[r] = x1; 3675 p_work2[r] = x2; 3676 k = rbeg[r]; 3677 idx = &ridx[k]; 3678 val = &rval[k]; 3679 3680 for(int m = rlen[r]; m != 0; --m) 3681 { 3682 vec1[*idx] -= x1 * (*val); 3683 vec2[*idx] -= x2 * (*val++); 3684 idx++; 3685 } 3686 } 3687 else if(x1 != 0) 3688 { 3689 x1 *= diag[r]; 3690 p_work1[r] = x1; 3691 k = rbeg[r]; 3692 idx = &ridx[k]; 3693 val = &rval[k]; 3694 3695 for(int m = rlen[r]; m != 0; --m) 3696 vec1[*idx++] -= x1 * (*val++); 3697 } 3698 else if(x2 != 0) 3699 { 3700 x2 *= diag[r]; 3701 p_work2[r] = x2; 3702 k = rbeg[r]; 3703 idx = &ridx[k]; 3704 val = &rval[k]; 3705 3706 for(int m = rlen[r]; m != 0; --m) 3707 vec2[*idx++] -= x2 * (*val++); 3708 } 3709 } 3710 } 3711 3712 inline int CLUFactorRational::solveLleft2forest(Rational* vec1, int* /* nonz */, Rational* vec2) 3713 { 3714 int i; 3715 int j; 3716 int k; 3717 int end; 3718 Rational x1, x2; 3719 Rational* val; 3720 int* lidx, *idx, *lrow; 3721 int* lbeg; 3722 3723 VectorRational& lval = l.val; 3724 lidx = l.idx; 3725 lrow = l.row; 3726 lbeg = l.start; 3727 3728 end = l.firstUpdate; 3729 3730 for(i = l.firstUnused - 1; i >= end; --i) 3731 { 3732 j = lrow[i]; 3733 x1 = vec1[j]; 3734 x2 = vec2[j]; 3735 3736 if(x1 != 0) 3737 { 3738 if(x2 != 0) 3739 { 3740 k = lbeg[i]; 3741 val = &lval[k]; 3742 idx = &lidx[k]; 3743 3744 for(j = lbeg[i + 1]; j > k; --j) 3745 { 3746 vec1[*idx] -= x1 * (*val); 3747 vec2[*idx++] -= x2 * (*val++); 3748 } 3749 } 3750 else 3751 { 3752 k = lbeg[i]; 3753 val = &lval[k]; 3754 idx = &lidx[k]; 3755 3756 for(j = lbeg[i + 1]; j > k; --j) 3757 vec1[*idx++] -= x1 * (*val++); 3758 } 3759 } 3760 else if(x2 != 0) 3761 { 3762 k = lbeg[i]; 3763 val = &lval[k]; 3764 idx = &lidx[k]; 3765 3766 for(j = lbeg[i + 1]; j > k; --j) 3767 vec2[*idx++] -= x2 * (*val++); 3768 } 3769 } 3770 3771 return 0; 3772 } 3773 3774 inline void CLUFactorRational::solveLleft2(Rational* vec1, int* /* nonz */, Rational* vec2) 3775 { 3776 int i, j, k, r; 3777 int x1not0, x2not0; 3778 Rational x1, x2; 3779 3780 Rational* val; 3781 int* ridx, *idx; 3782 int* rbeg; 3783 int* rorig; 3784 3785 ridx = l.ridx; 3786 VectorRational& rval = l.rval; 3787 rbeg = l.rbeg; 3788 rorig = l.rorig; 3789 3790 #ifndef SOPLEX_WITH_L_ROWS 3791 VectorRational& lval = l.val; 3792 int* lidx = l.idx; 3793 int* lrow = l.row; 3794 int* lbeg = l.start; 3795 3796 i = l.firstUpdate - 1; 3797 3798 for(; i >= 0; --i) 3799 { 3800 k = lbeg[i]; 3801 val = &lval[k]; 3802 idx = &lidx[k]; 3803 x1 = 0; 3804 x2 = 0; 3805 3806 for(j = lbeg[i + 1]; j > k; --j) 3807 { 3808 x1 += vec1[*idx] * (*val); 3809 x2 += vec2[*idx++] * (*val++); 3810 } 3811 3812 vec1[lrow[i]] -= x1; 3813 3814 vec2[lrow[i]] -= x2; 3815 } 3816 3817 #else 3818 3819 for(i = thedim; i--;) 3820 { 3821 r = rorig[i]; 3822 x1 = vec1[r]; 3823 x2 = vec2[r]; 3824 x1not0 = (x1 != 0); 3825 x2not0 = (x2 != 0); 3826 3827 if(x1not0 && x2not0) 3828 { 3829 k = rbeg[r]; 3830 j = rbeg[r + 1] - k; 3831 val = &rval[k]; 3832 idx = &ridx[k]; 3833 3834 while(j-- > 0) 3835 { 3836 assert(row.perm[*idx] < i); 3837 vec1[*idx] -= x1 * *val; 3838 vec2[*idx++] -= x2 * *val++; 3839 } 3840 } 3841 else if(x1not0) 3842 { 3843 k = rbeg[r]; 3844 j = rbeg[r + 1] - k; 3845 val = &rval[k]; 3846 idx = &ridx[k]; 3847 3848 while(j-- > 0) 3849 { 3850 assert(row.perm[*idx] < i); 3851 vec1[*idx++] -= x1 * *val++; 3852 } 3853 } 3854 else if(x2not0) 3855 { 3856 k = rbeg[r]; 3857 j = rbeg[r + 1] - k; 3858 val = &rval[k]; 3859 idx = &ridx[k]; 3860 3861 while(j-- > 0) 3862 { 3863 assert(row.perm[*idx] < i); 3864 vec2[*idx++] -= x2 * *val++; 3865 } 3866 } 3867 } 3868 3869 #endif 3870 } 3871 3872 inline int CLUFactorRational::solveLleftForest(Rational* vec, int* /* nonz */) 3873 { 3874 int i, j, k, end; 3875 Rational x; 3876 Rational* val; 3877 int* idx, *lidx, *lrow, *lbeg; 3878 3879 VectorRational& lval = l.val; 3880 lidx = l.idx; 3881 lrow = l.row; 3882 lbeg = l.start; 3883 3884 end = l.firstUpdate; 3885 3886 for(i = l.firstUnused - 1; i >= end; --i) 3887 { 3888 if((x = vec[lrow[i]]) != 0) 3889 { 3890 if(timeLimitReached()) 3891 return 0; 3892 3893 k = lbeg[i]; 3894 val = &lval[k]; 3895 idx = &lidx[k]; 3896 3897 for(j = lbeg[i + 1]; j > k; --j) 3898 vec[*idx++] -= x * (*val++); 3899 } 3900 } 3901 3902 return 0; 3903 } 3904 3905 inline void CLUFactorRational::solveLleft(Rational* vec) 3906 { 3907 3908 #ifndef SOPLEX_WITH_L_ROWS 3909 int* idx; 3910 Rational* val; 3911 VectorRational& lval = l.val; 3912 int* lidx = l.idx; 3913 int* lrow = l.row; 3914 int* lbeg = l.start; 3915 3916 for(int i = l.firstUpdate - 1; i >= 0; --i) 3917 { 3918 if(timeLimitReached()) 3919 return; 3920 3921 int k = lbeg[i]; 3922 val = &lval[k]; 3923 idx = &lidx[k]; 3924 x = 0; 3925 3926 for(int j = lbeg[i + 1]; j > k; --j) 3927 x += vec[*idx++] * (*val++); 3928 3929 vec[lrow[i]] -= x; 3930 } 3931 3932 #else 3933 3934 for(int i = thedim - 1; i >= 0; --i) 3935 { 3936 int r = l.rorig[i]; 3937 Rational x = vec[r]; 3938 3939 if(x != 0) 3940 { 3941 if(timeLimitReached()) 3942 return; 3943 3944 for(int k = l.rbeg[r]; k < l.rbeg[r + 1]; k++) 3945 { 3946 int j = l.ridx[k]; 3947 3948 assert(l.rperm[j] < i); 3949 3950 vec[j] -= x * l.rval[k]; 3951 } 3952 } 3953 } 3954 3955 #endif // SOPLEX_WITH_L_ROWS 3956 } 3957 3958 inline int CLUFactorRational::solveLleftEps(Rational* vec, int* nonz) 3959 { 3960 int i, j, k, n; 3961 int r; 3962 Rational x; 3963 Rational* val; 3964 int* ridx, *idx; 3965 int* rbeg; 3966 int* rorig; 3967 3968 ridx = l.ridx; 3969 VectorRational& rval = l.rval; 3970 rbeg = l.rbeg; 3971 rorig = l.rorig; 3972 n = 0; 3973 #ifndef SOPLEX_WITH_L_ROWS 3974 VectorRational& lval = l.val; 3975 int* lidx = l.idx; 3976 int* lrow = l.row; 3977 int* lbeg = l.start; 3978 3979 for(i = l.firstUpdate - 1; i >= 0; --i) 3980 { 3981 k = lbeg[i]; 3982 val = &lval[k]; 3983 idx = &lidx[k]; 3984 x = 0; 3985 3986 for(j = lbeg[i + 1]; j > k; --j) 3987 x += vec[*idx++] * (*val++); 3988 3989 vec[lrow[i]] -= x; 3990 } 3991 3992 #else 3993 3994 for(i = thedim; i--;) 3995 { 3996 r = rorig[i]; 3997 x = vec[r]; 3998 3999 if(x != 0) 4000 { 4001 *nonz++ = r; 4002 n++; 4003 k = rbeg[r]; 4004 j = rbeg[r + 1] - k; 4005 val = &rval[k]; 4006 idx = &ridx[k]; 4007 4008 while(j-- > 0) 4009 { 4010 assert(row.perm[*idx] < i); 4011 vec[*idx++] -= x * *val++; 4012 } 4013 } 4014 else 4015 vec[r] = 0; 4016 } 4017 4018 #endif 4019 4020 return n; 4021 } 4022 4023 inline void CLUFactorRational::solveUpdateLeft(Rational* vec) 4024 { 4025 int i, j, k, end; 4026 Rational x; 4027 Rational* val; 4028 int* lrow, *lidx, *idx; 4029 int* lbeg; 4030 4031 VectorRational& lval = l.val; 4032 lidx = l.idx; 4033 lrow = l.row; 4034 lbeg = l.start; 4035 4036 assert(!l.updateType); /* Forest-Tomlin Updates */ 4037 4038 end = l.firstUpdate; 4039 4040 for(i = l.firstUnused - 1; i >= end; --i) 4041 { 4042 k = lbeg[i]; 4043 val = &lval[k]; 4044 idx = &lidx[k]; 4045 x = 0; 4046 4047 for(j = lbeg[i + 1]; j > k; --j) 4048 x += vec[*idx++] * (*val++); 4049 4050 vec[lrow[i]] -= x; 4051 } 4052 } 4053 4054 inline void CLUFactorRational::solveUpdateLeft2(Rational* vec1, Rational* vec2) 4055 { 4056 int i, j, k, end; 4057 Rational x1, x2; 4058 Rational* val; 4059 int* lrow, *lidx, *idx; 4060 int* lbeg; 4061 4062 VectorRational& lval = l.val; 4063 lidx = l.idx; 4064 lrow = l.row; 4065 lbeg = l.start; 4066 4067 assert(!l.updateType); /* Forest-Tomlin Updates */ 4068 4069 end = l.firstUpdate; 4070 4071 for(i = l.firstUnused - 1; i >= end; --i) 4072 { 4073 k = lbeg[i]; 4074 val = &lval[k]; 4075 idx = &lidx[k]; 4076 x1 = 0; 4077 x2 = 0; 4078 4079 for(j = lbeg[i + 1]; j > k; --j) 4080 { 4081 x1 += vec1[*idx] * (*val); 4082 x2 += vec2[*idx++] * (*val++); 4083 } 4084 4085 vec1[lrow[i]] -= x1; 4086 4087 vec2[lrow[i]] -= x2; 4088 } 4089 } 4090 4091 inline int CLUFactorRational::solveUpdateLeft(Rational* vec, int* nonz, int n) 4092 { 4093 int i, j, k, end; 4094 Rational x, y; 4095 Rational* val; 4096 int* lrow, *lidx, *idx; 4097 int* lbeg; 4098 4099 assert(!l.updateType); /* no Forest-Tomlin Updates! */ 4100 4101 VectorRational& lval = l.val; 4102 lidx = l.idx; 4103 lrow = l.row; 4104 lbeg = l.start; 4105 4106 end = l.firstUpdate; 4107 4108 for(i = l.firstUnused - 1; i >= end; --i) 4109 { 4110 k = lbeg[i]; 4111 assert(k >= 0 && k < l.val.dim()); 4112 val = &lval[k]; 4113 idx = &lidx[k]; 4114 x = 0; 4115 4116 for(j = lbeg[i + 1]; j > k; --j) 4117 { 4118 assert(*idx >= 0 && *idx < thedim); 4119 x += vec[*idx++] * (*val++); 4120 } 4121 4122 k = lrow[i]; 4123 4124 y = vec[k]; 4125 4126 if(y == 0) 4127 { 4128 y = -x; 4129 4130 if(y != 0) 4131 { 4132 nonz[n++] = k; 4133 vec[k] = y; 4134 } 4135 } 4136 else 4137 { 4138 y -= x; 4139 //vec[k] = ( y != 0 ) ? y : SOPLEX_MARKER; 4140 vec[k] = y; 4141 } 4142 } 4143 4144 return n; 4145 } 4146 4147 inline void CLUFactorRational::solveLeft(Rational* vec, Rational* rhs) 4148 { 4149 4150 if(!l.updateType) /* no Forest-Tomlin Updates */ 4151 { 4152 solveUpdateLeft(rhs); 4153 solveUleft(vec, rhs); 4154 solveLleft(vec); 4155 } 4156 else 4157 { 4158 solveUleft(vec, rhs); 4159 solveLleftForest(vec, 0); 4160 solveLleft(vec); 4161 } 4162 } 4163 4164 inline int CLUFactorRational::solveLeftEps(Rational* vec, Rational* rhs, int* nonz) 4165 { 4166 4167 if(!l.updateType) /* no Forest-Tomlin Updates */ 4168 { 4169 solveUpdateLeft(rhs); 4170 solveUleft(vec, rhs); 4171 return solveLleftEps(vec, nonz); 4172 } 4173 else 4174 { 4175 solveUleft(vec, rhs); 4176 solveLleftForest(vec, nonz); 4177 return solveLleftEps(vec, nonz); 4178 } 4179 } 4180 4181 inline int CLUFactorRational::solveLeft2( 4182 Rational* vec1, 4183 int* nonz, 4184 Rational* vec2, 4185 Rational* rhs1, 4186 Rational* rhs2) 4187 { 4188 4189 if(!l.updateType) /* no Forest-Tomlin Updates */ 4190 { 4191 solveUpdateLeft2(rhs1, rhs2); 4192 solveUleft2(vec1, rhs1, vec2, rhs2); 4193 solveLleft2(vec1, nonz, vec2); 4194 return 0; 4195 } 4196 else 4197 { 4198 solveUleft2(vec1, rhs1, vec2, rhs2); 4199 solveLleft2forest(vec1, nonz, vec2); 4200 solveLleft2(vec1, nonz, vec2); 4201 return 0; 4202 } 4203 } 4204 4205 inline int CLUFactorRational::solveUleft(Rational* vec, int* vecidx, 4206 Rational* rhs, int* rhsidx, int rhsn) 4207 { 4208 Rational x, y; 4209 int i, j, k, n, r, c; 4210 int* rorig, *corig, *cperm; 4211 int* ridx, *rlen, *rbeg, *idx; 4212 Rational* val; 4213 4214 rorig = row.orig; 4215 corig = col.orig; 4216 cperm = col.perm; 4217 4218 /* move rhsidx to a heap 4219 */ 4220 4221 for(i = 0; i < rhsn;) 4222 enQueueMinRat(rhsidx, &i, cperm[rhsidx[i]]); 4223 4224 ridx = u.row.idx; 4225 4226 VectorRational& rval = u.row.val; 4227 4228 rlen = u.row.len; 4229 4230 rbeg = u.row.start; 4231 4232 n = 0; 4233 4234 while(rhsn > 0) 4235 { 4236 i = deQueueMinRat(rhsidx, &rhsn); 4237 assert(i >= 0 && i < thedim); 4238 c = corig[i]; 4239 assert(c >= 0 && c < thedim); 4240 x = rhs[c]; 4241 rhs[c] = 0; 4242 4243 if(x != 0) 4244 { 4245 r = rorig[i]; 4246 assert(r >= 0 && r < thedim); 4247 vecidx[n++] = r; 4248 x *= diag[r]; 4249 vec[r] = x; 4250 k = rbeg[r]; 4251 assert(k >= 0 && k < u.row.val.dim()); 4252 idx = &ridx[k]; 4253 val = &rval[k]; 4254 4255 for(int m = rlen[r]; m; --m) 4256 { 4257 j = *idx++; 4258 assert(j >= 0 && j < thedim); 4259 y = rhs[j]; 4260 4261 if(y == 0) 4262 { 4263 y = -x * (*val++); 4264 4265 if(y != 0) 4266 { 4267 rhs[j] = y; 4268 enQueueMinRat(rhsidx, &rhsn, cperm[j]); 4269 } 4270 } 4271 else 4272 { 4273 y -= x * (*val++); 4274 // rhs[j] = ( y != 0 ) ? y : SOPLEX_MARKER; 4275 rhs[j] = y; 4276 } 4277 } 4278 } 4279 } 4280 4281 return n; 4282 } 4283 4284 4285 inline void CLUFactorRational::solveUleftNoNZ(Rational* vec, Rational* rhs, int* rhsidx, int rhsn) 4286 { 4287 Rational x, y; 4288 int i, j, k, r, c; 4289 int* rorig, *corig, *cperm; 4290 int* ridx, *rlen, *rbeg, *idx; 4291 Rational* val; 4292 4293 rorig = row.orig; 4294 corig = col.orig; 4295 cperm = col.perm; 4296 4297 /* move rhsidx to a heap 4298 */ 4299 4300 for(i = 0; i < rhsn;) 4301 enQueueMinRat(rhsidx, &i, cperm[rhsidx[i]]); 4302 4303 ridx = u.row.idx; 4304 4305 VectorRational& rval = u.row.val; 4306 4307 rlen = u.row.len; 4308 4309 rbeg = u.row.start; 4310 4311 while(rhsn > 0) 4312 { 4313 i = deQueueMinRat(rhsidx, &rhsn); 4314 assert(i >= 0 && i < thedim); 4315 c = corig[i]; 4316 assert(c >= 0 && c < thedim); 4317 x = rhs[c]; 4318 rhs[c] = 0; 4319 4320 if(x != 0) 4321 { 4322 r = rorig[i]; 4323 assert(r >= 0 && r < thedim); 4324 x *= diag[r]; 4325 vec[r] = x; 4326 k = rbeg[r]; 4327 assert(k >= 0 && k < u.row.val.dim()); 4328 idx = &ridx[k]; 4329 val = &rval[k]; 4330 4331 for(int m = rlen[r]; m; --m) 4332 { 4333 j = *idx++; 4334 assert(j >= 0 && j < thedim); 4335 y = rhs[j]; 4336 4337 if(y == 0) 4338 { 4339 y = -x * (*val++); 4340 4341 if(y != 0) 4342 { 4343 rhs[j] = y; 4344 enQueueMinRat(rhsidx, &rhsn, cperm[j]); 4345 } 4346 } 4347 else 4348 { 4349 y -= x * (*val++); 4350 // rhs[j] = ( y != 0 ) ? y : SOPLEX_MARKER; 4351 rhs[j] = y; 4352 } 4353 } 4354 } 4355 } 4356 } 4357 4358 4359 inline int CLUFactorRational::solveLleftForest(Rational* vec, int* nonz, int n) 4360 { 4361 int i, j, k, end; 4362 Rational x, y; 4363 Rational* val; 4364 int* idx, *lidx, *lrow, *lbeg; 4365 4366 VectorRational& lval = l.val; 4367 lidx = l.idx; 4368 lrow = l.row; 4369 lbeg = l.start; 4370 end = l.firstUpdate; 4371 4372 for(i = l.firstUnused - 1; i >= end; --i) 4373 { 4374 assert(i >= 0 && i < l.val.dim()); 4375 4376 if((x = vec[lrow[i]]) != 0) 4377 { 4378 k = lbeg[i]; 4379 assert(k >= 0 && k < l.val.dim()); 4380 val = &lval[k]; 4381 idx = &lidx[k]; 4382 4383 for(j = lbeg[i + 1]; j > k; --j) 4384 { 4385 int m = *idx++; 4386 assert(m >= 0 && m < thedim); 4387 y = vec[m]; 4388 4389 if(y == 0) 4390 { 4391 y = -x * (*val++); 4392 4393 if(y != 0) 4394 { 4395 vec[m] = y; 4396 nonz[n++] = m; 4397 } 4398 } 4399 else 4400 { 4401 y -= x * (*val++); 4402 4403 } 4404 } 4405 } 4406 } 4407 4408 return n; 4409 } 4410 4411 4412 inline void CLUFactorRational::solveLleftForestNoNZ(Rational* vec) 4413 { 4414 int i, j, k, end; 4415 Rational x; 4416 Rational* val; 4417 int* idx, *lidx, *lrow, *lbeg; 4418 4419 VectorRational& lval = l.val; 4420 lidx = l.idx; 4421 lrow = l.row; 4422 lbeg = l.start; 4423 end = l.firstUpdate; 4424 4425 for(i = l.firstUnused - 1; i >= end; --i) 4426 { 4427 if((x = vec[lrow[i]]) != 0) 4428 { 4429 assert(i >= 0 && i < l.val.dim()); 4430 k = lbeg[i]; 4431 assert(k >= 0 && k < l.val.dim()); 4432 val = &lval[k]; 4433 idx = &lidx[k]; 4434 4435 for(j = lbeg[i + 1]; j > k; --j) 4436 { 4437 assert(*idx >= 0 && *idx < thedim); 4438 vec[*idx++] -= x * (*val++); 4439 } 4440 } 4441 } 4442 } 4443 4444 4445 inline int CLUFactorRational::solveLleft(Rational* vec, int* nonz, int rn) 4446 { 4447 int i, j, k, n; 4448 int r; 4449 Rational x, y; 4450 Rational* val; 4451 int* ridx, *idx; 4452 int* rbeg; 4453 int* rorig, *rperm; 4454 int* last; 4455 4456 ridx = l.ridx; 4457 VectorRational& rval = l.rval; 4458 rbeg = l.rbeg; 4459 rorig = l.rorig; 4460 rperm = l.rperm; 4461 n = 0; 4462 4463 i = l.firstUpdate - 1; 4464 #ifndef SOPLEX_WITH_L_ROWS 4465 #pragma warn "Not yet implemented, define SOPLEX_WITH_L_ROWS" 4466 VectorRational& lval = l.val; 4467 int* lidx = l.idx; 4468 int* lrow = l.row; 4469 int* lbeg = l.start; 4470 4471 for(; i >= 0; --i) 4472 { 4473 k = lbeg[i]; 4474 val = &lval[k]; 4475 idx = &lidx[k]; 4476 x = 0; 4477 4478 for(j = lbeg[i + 1]; j > k; --j) 4479 x += vec[*idx++] * (*val++); 4480 4481 vec[lrow[i]] -= x; 4482 } 4483 4484 #else 4485 4486 /* move rhsidx to a heap 4487 */ 4488 for(i = 0; i < rn;) 4489 enQueueMaxRat(nonz, &i, rperm[nonz[i]]); 4490 4491 last = nonz + thedim; 4492 4493 while(rn > 0) 4494 { 4495 i = deQueueMaxRat(nonz, &rn); 4496 r = rorig[i]; 4497 x = vec[r]; 4498 4499 if(x != 0) 4500 { 4501 *(--last) = r; 4502 n++; 4503 k = rbeg[r]; 4504 j = rbeg[r + 1] - k; 4505 val = &rval[k]; 4506 idx = &ridx[k]; 4507 4508 while(j-- > 0) 4509 { 4510 assert(l.rperm[*idx] < i); 4511 int m = *idx++; 4512 y = vec[m]; 4513 4514 if(y == 0) 4515 { 4516 y = -x * *val++; 4517 4518 if(y != 0) 4519 { 4520 vec[m] = y; 4521 enQueueMaxRat(nonz, &rn, rperm[m]); 4522 } 4523 } 4524 else 4525 { 4526 y -= x * *val++; 4527 // vec[m] = ( y != 0 ) ? y : SOPLEX_MARKER; 4528 vec[m] = y; 4529 } 4530 } 4531 } 4532 else 4533 vec[r] = 0; 4534 } 4535 4536 for(i = 0; i < n; ++i) 4537 *nonz++ = *last++; 4538 4539 #endif 4540 4541 return n; 4542 } 4543 4544 4545 inline void CLUFactorRational::solveLleftNoNZ(Rational* vec) 4546 { 4547 int i, j, k; 4548 int r; 4549 Rational x; 4550 Rational* val; 4551 int* ridx, *idx; 4552 int* rbeg; 4553 int* rorig; 4554 4555 ridx = l.ridx; 4556 VectorRational& rval = l.rval; 4557 rbeg = l.rbeg; 4558 rorig = l.rorig; 4559 4560 #ifndef SOPLEX_WITH_L_ROWS 4561 VectorRational& lval = l.val; 4562 int* lidx = l.idx; 4563 int* lrow = l.row; 4564 int* lbeg = l.start; 4565 4566 i = l.firstUpdate - 1; 4567 assert(i < thedim); 4568 4569 for(; i >= 0; --i) 4570 { 4571 k = lbeg[i]; 4572 assert(k >= 0 && k < l.val.dim()); 4573 val = &lval[k]; 4574 idx = &lidx[k]; 4575 x = 0; 4576 4577 for(j = lbeg[i + 1]; j > k; --j) 4578 { 4579 assert(*idx >= 0 && *idx < thedim); 4580 x += vec[*idx++] * (*val++); 4581 } 4582 4583 vec[lrow[i]] -= x; 4584 } 4585 4586 #else 4587 4588 for(i = thedim; i--;) 4589 { 4590 r = rorig[i]; 4591 x = vec[r]; 4592 4593 if(x != 0) 4594 { 4595 k = rbeg[r]; 4596 j = rbeg[r + 1] - k; 4597 val = &rval[k]; 4598 idx = &ridx[k]; 4599 4600 while(j-- > 0) 4601 { 4602 assert(l.rperm[*idx] < i); 4603 vec[*idx++] -= x * *val++; 4604 } 4605 } 4606 } 4607 4608 #endif 4609 } 4610 4611 inline int CLUFactorRational::vSolveLright(Rational* vec, int* ridx, int rn) 4612 { 4613 int i, j, k, n; 4614 int end; 4615 Rational x; 4616 Rational* val; 4617 int* lrow, *lidx, *idx; 4618 int* lbeg; 4619 4620 VectorRational& lval = l.val; 4621 lidx = l.idx; 4622 lrow = l.row; 4623 lbeg = l.start; 4624 4625 end = l.firstUpdate; 4626 4627 for(i = 0; i < end; ++i) 4628 { 4629 x = vec[lrow[i]]; 4630 4631 if(x != 0) 4632 { 4633 k = lbeg[i]; 4634 idx = &(lidx[k]); 4635 val = &(lval[k]); 4636 4637 for(j = lbeg[i + 1]; j > k; --j) 4638 { 4639 assert(*idx >= 0 && *idx < thedim); 4640 ridx[rn] = n = *idx++; 4641 rn += (vec[n] == 0) ? 1 : 0; 4642 vec[n] -= x * (*val++); 4643 // vec[n] += ( vec[n] == 0 ) ? SOPLEX_MARKER : 0; 4644 } 4645 } 4646 } 4647 4648 if(l.updateType) /* Forest-Tomlin Updates */ 4649 { 4650 end = l.firstUnused; 4651 4652 for(; i < end; ++i) 4653 { 4654 x = 0; 4655 k = lbeg[i]; 4656 idx = &(lidx[k]); 4657 val = &(lval[k]); 4658 4659 for(j = lbeg[i + 1]; j > k; --j) 4660 { 4661 assert(*idx >= 0 && *idx < thedim); 4662 x += vec[*idx++] * (*val++); 4663 } 4664 4665 ridx[rn] = j = lrow[i]; 4666 4667 rn += (vec[j] == 0) ? 1 : 0; 4668 vec[j] -= x; 4669 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0; 4670 } 4671 } 4672 4673 return rn; 4674 } 4675 4676 inline void CLUFactorRational::vSolveLright2(Rational* vec, int* ridx, int* rnptr, 4677 Rational* vec2, int* ridx2, int* rn2ptr) 4678 { 4679 int i, j, k, n; 4680 int end; 4681 Rational x, y; 4682 Rational x2, y2; 4683 Rational* val; 4684 int* lrow, *lidx, *idx; 4685 int* lbeg; 4686 4687 int rn = *rnptr; 4688 int rn2 = *rn2ptr; 4689 4690 VectorRational& lval = l.val; 4691 lidx = l.idx; 4692 lrow = l.row; 4693 lbeg = l.start; 4694 4695 end = l.firstUpdate; 4696 4697 for(i = 0; i < end; ++i) 4698 { 4699 j = lrow[i]; 4700 x2 = vec2[j]; 4701 x = vec[j]; 4702 4703 if(x != 0) 4704 { 4705 k = lbeg[i]; 4706 idx = &(lidx[k]); 4707 val = &(lval[k]); 4708 4709 if(x2 != 0) 4710 { 4711 for(j = lbeg[i + 1]; j > k; --j) 4712 { 4713 assert(*idx >= 0 && *idx < thedim); 4714 ridx[rn] = ridx2[rn2] = n = *idx++; 4715 y = vec[n]; 4716 y2 = vec2[n]; 4717 rn += (y == 0) ? 1 : 0; 4718 rn2 += (y2 == 0) ? 1 : 0; 4719 y -= x * (*val); 4720 y2 -= x2 * (*val++); 4721 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4722 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4723 vec[n] = y; 4724 vec2[n] = y2; 4725 } 4726 } 4727 else 4728 { 4729 for(j = lbeg[i + 1]; j > k; --j) 4730 { 4731 assert(*idx >= 0 && *idx < thedim); 4732 ridx[rn] = n = *idx++; 4733 y = vec[n]; 4734 rn += (y == 0) ? 1 : 0; 4735 y -= x * (*val++); 4736 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4737 vec[n] = y; 4738 } 4739 } 4740 } 4741 else if(x2 != 0) 4742 { 4743 k = lbeg[i]; 4744 idx = &(lidx[k]); 4745 val = &(lval[k]); 4746 4747 for(j = lbeg[i + 1]; j > k; --j) 4748 { 4749 assert(*idx >= 0 && *idx < thedim); 4750 ridx2[rn2] = n = *idx++; 4751 y2 = vec2[n]; 4752 rn2 += (y2 == 0) ? 1 : 0; 4753 y2 -= x2 * (*val++); 4754 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4755 vec2[n] = y2; 4756 } 4757 } 4758 } 4759 4760 if(l.updateType) /* Forest-Tomlin Updates */ 4761 { 4762 end = l.firstUnused; 4763 4764 for(; i < end; ++i) 4765 { 4766 x = x2 = 0; 4767 k = lbeg[i]; 4768 idx = &(lidx[k]); 4769 val = &(lval[k]); 4770 4771 for(j = lbeg[i + 1]; j > k; --j) 4772 { 4773 assert(*idx >= 0 && *idx < thedim); 4774 x += vec[*idx] * (*val); 4775 x2 += vec2[*idx++] * (*val++); 4776 } 4777 4778 ridx[rn] = ridx2[rn2] = j = lrow[i]; 4779 4780 rn += (vec[j] == 0) ? 1 : 0; 4781 rn2 += (vec2[j] == 0) ? 1 : 0; 4782 vec[j] -= x; 4783 vec2[j] -= x2; 4784 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0; 4785 // vec2[j] += ( vec2[j] == 0 ) ? SOPLEX_MARKER : 0; 4786 } 4787 } 4788 4789 *rnptr = rn; 4790 4791 *rn2ptr = rn2; 4792 } 4793 4794 inline void CLUFactorRational::vSolveLright3(Rational* vec, int* ridx, int* rnptr, 4795 Rational* vec2, int* ridx2, int* rn2ptr, 4796 Rational* vec3, int* ridx3, int* rn3ptr) 4797 { 4798 int i, j, k, n; 4799 int end; 4800 Rational x, y; 4801 Rational x2, y2; 4802 Rational x3, y3; 4803 Rational* val; 4804 int* lrow, *lidx, *idx; 4805 int* lbeg; 4806 4807 int rn = *rnptr; 4808 int rn2 = *rn2ptr; 4809 int rn3 = *rn3ptr; 4810 4811 VectorRational& lval = l.val; 4812 lidx = l.idx; 4813 lrow = l.row; 4814 lbeg = l.start; 4815 4816 end = l.firstUpdate; 4817 4818 for(i = 0; i < end; ++i) 4819 { 4820 j = lrow[i]; 4821 x3 = vec3[j]; 4822 x2 = vec2[j]; 4823 x = vec[j]; 4824 4825 if(x != 0) 4826 { 4827 k = lbeg[i]; 4828 idx = &(lidx[k]); 4829 val = &(lval[k]); 4830 4831 if(x2 != 0) 4832 { 4833 if(x3 != 0) 4834 { 4835 // case 1: all three vectors are nonzero at j 4836 for(j = lbeg[i + 1]; j > k; --j) 4837 { 4838 assert(*idx >= 0 && *idx < thedim); 4839 ridx[rn] = ridx2[rn2] = ridx3[rn3] = n = *idx++; 4840 y = vec[n]; 4841 y2 = vec2[n]; 4842 y3 = vec3[n]; 4843 rn += (y == 0) ? 1 : 0; 4844 rn2 += (y2 == 0) ? 1 : 0; 4845 rn3 += (y3 == 0) ? 1 : 0; 4846 y -= x * (*val); 4847 y2 -= x2 * (*val); 4848 y3 -= x3 * (*val++); 4849 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4850 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4851 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 ); 4852 vec[n] = y; 4853 vec2[n] = y2; 4854 vec3[n] = y3; 4855 } 4856 } 4857 else 4858 { 4859 // case 2: 1 and 2 are nonzero at j 4860 for(j = lbeg[i + 1]; j > k; --j) 4861 { 4862 assert(*idx >= 0 && *idx < thedim); 4863 ridx[rn] = ridx2[rn2] = n = *idx++; 4864 y = vec[n]; 4865 y2 = vec2[n]; 4866 rn += (y == 0) ? 1 : 0; 4867 rn2 += (y2 == 0) ? 1 : 0; 4868 y -= x * (*val); 4869 y2 -= x2 * (*val++); 4870 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4871 //vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4872 vec[n] = y; 4873 vec2[n] = y2; 4874 } 4875 } 4876 } 4877 else if(x3 != 0) 4878 { 4879 // case 3: 1 and 3 are nonzero at j 4880 for(j = lbeg[i + 1]; j > k; --j) 4881 { 4882 assert(*idx >= 0 && *idx < thedim); 4883 ridx[rn] = ridx3[rn3] = n = *idx++; 4884 y = vec[n]; 4885 y3 = vec3[n]; 4886 rn += (y == 0) ? 1 : 0; 4887 rn3 += (y3 == 0) ? 1 : 0; 4888 y -= x * (*val); 4889 y3 -= x3 * (*val++); 4890 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4891 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 ); 4892 vec[n] = y; 4893 vec3[n] = y3; 4894 } 4895 } 4896 else 4897 { 4898 // case 4: only 1 is nonzero at j 4899 for(j = lbeg[i + 1]; j > k; --j) 4900 { 4901 assert(*idx >= 0 && *idx < thedim); 4902 ridx[rn] = n = *idx++; 4903 y = vec[n]; 4904 rn += (y == 0) ? 1 : 0; 4905 y -= x * (*val++); 4906 // vec[n] = y + ( y == 0 ? SOPLEX_MARKER : 0 ); 4907 vec[n] = y; 4908 } 4909 } 4910 } 4911 else if(x2 != 0) 4912 { 4913 k = lbeg[i]; 4914 idx = &(lidx[k]); 4915 val = &(lval[k]); 4916 4917 if(x3 != 0) 4918 { 4919 // case 5: 2 and 3 are nonzero at j 4920 for(j = lbeg[i + 1]; j > k; --j) 4921 { 4922 assert(*idx >= 0 && *idx < thedim); 4923 ridx2[rn2] = ridx3[rn3] = n = *idx++; 4924 y2 = vec2[n]; 4925 y3 = vec3[n]; 4926 rn2 += (y2 == 0) ? 1 : 0; 4927 rn3 += (y3 == 0) ? 1 : 0; 4928 y2 -= x2 * (*val); 4929 y3 -= x3 * (*val++); 4930 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4931 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 ); 4932 vec2[n] = y2; 4933 vec3[n] = y3; 4934 } 4935 } 4936 else 4937 { 4938 // case 6: only 2 is nonzero at j 4939 for(j = lbeg[i + 1]; j > k; --j) 4940 { 4941 assert(*idx >= 0 && *idx < thedim); 4942 ridx2[rn2] = n = *idx++; 4943 y2 = vec2[n]; 4944 rn2 += (y2 == 0) ? 1 : 0; 4945 y2 -= x2 * (*val++); 4946 // vec2[n] = y2 + ( y2 == 0 ? SOPLEX_MARKER : 0 ); 4947 vec2[n] = y2; 4948 } 4949 } 4950 } 4951 else if(x3 != 0) 4952 { 4953 // case 7: only 3 is nonzero at j 4954 k = lbeg[i]; 4955 idx = &(lidx[k]); 4956 val = &(lval[k]); 4957 4958 for(j = lbeg[i + 1]; j > k; --j) 4959 { 4960 assert(*idx >= 0 && *idx < thedim); 4961 ridx3[rn3] = n = *idx++; 4962 y3 = vec3[n]; 4963 rn3 += (y3 == 0) ? 1 : 0; 4964 y3 -= x3 * (*val++); 4965 // vec3[n] = y3 + ( y3 == 0 ? SOPLEX_MARKER : 0 ); 4966 vec3[n] = y3; 4967 } 4968 } 4969 } 4970 4971 if(l.updateType) /* Forest-Tomlin Updates */ 4972 { 4973 end = l.firstUnused; 4974 4975 for(; i < end; ++i) 4976 { 4977 x = x2 = x3 = 0; 4978 k = lbeg[i]; 4979 idx = &(lidx[k]); 4980 val = &(lval[k]); 4981 4982 for(j = lbeg[i + 1]; j > k; --j) 4983 { 4984 assert(*idx >= 0 && *idx < thedim); 4985 x += vec[*idx] * (*val); 4986 x2 += vec2[*idx] * (*val); 4987 x3 += vec3[*idx++] * (*val++); 4988 } 4989 4990 ridx[rn] = ridx2[rn2] = ridx3[rn3] = j = lrow[i]; 4991 4992 rn += (vec[j] == 0) ? 1 : 0; 4993 rn2 += (vec2[j] == 0) ? 1 : 0; 4994 rn3 += (vec3[j] == 0) ? 1 : 0; 4995 vec[j] -= x; 4996 vec2[j] -= x2; 4997 vec3[j] -= x3; 4998 // vec[j] += ( vec[j] == 0 ) ? SOPLEX_MARKER : 0; 4999 // vec2[j] += ( vec2[j] == 0 ) ? SOPLEX_MARKER : 0; 5000 // vec3[j] += ( vec3[j] == 0 ) ? SOPLEX_MARKER : 0; 5001 } 5002 } 5003 5004 *rnptr = rn; 5005 5006 *rn2ptr = rn2; 5007 *rn3ptr = rn3; 5008 } 5009 5010 inline int CLUFactorRational::vSolveUright(Rational* vec, int* vidx, 5011 Rational* rhs, int* ridx, int rn) 5012 { 5013 int i, j, k, r, c, n; 5014 int* rorig, *corig; 5015 int* rperm; 5016 int* cidx, *clen, *cbeg; 5017 Rational x, y; 5018 5019 int* idx; 5020 Rational* val; 5021 5022 rorig = row.orig; 5023 corig = col.orig; 5024 rperm = row.perm; 5025 5026 cidx = u.col.idx; 5027 VectorRational& cval = u.col.val; 5028 clen = u.col.len; 5029 cbeg = u.col.start; 5030 5031 n = 0; 5032 5033 while(rn > 0) 5034 { 5035 /* Find nonzero with highest permuted row index and setup i and r 5036 */ 5037 i = deQueueMaxRat(ridx, &rn); 5038 assert(i >= 0 && i < thedim); 5039 r = rorig[i]; 5040 assert(r >= 0 && r < thedim); 5041 5042 x = diag[r] * rhs[r]; 5043 rhs[r] = 0; 5044 5045 if(x != 0) 5046 { 5047 c = corig[i]; 5048 assert(c >= 0 && c < thedim); 5049 vidx[n++] = c; 5050 vec[c] = x; 5051 val = &cval[cbeg[c]]; 5052 idx = &cidx[cbeg[c]]; 5053 j = clen[c]; 5054 5055 while(j-- > 0) 5056 { 5057 assert(*idx >= 0 && *idx < thedim); 5058 k = *idx++; 5059 assert(k >= 0 && k < thedim); 5060 y = rhs[k]; 5061 5062 if(y == 0) 5063 { 5064 y = -x * (*val++); 5065 5066 if(y != 0) 5067 { 5068 rhs[k] = y; 5069 enQueueMaxRat(ridx, &rn, rperm[k]); 5070 } 5071 } 5072 else 5073 { 5074 y -= x * (*val++); 5075 // y += ( y == 0 ) ? SOPLEX_MARKER : 0; 5076 rhs[k] = y; 5077 } 5078 } 5079 5080 if(rn > i * verySparseFactor4rightRat) 5081 { 5082 /* continue with dense case */ 5083 for(i = *ridx; i >= 0; --i) 5084 { 5085 r = rorig[i]; 5086 assert(r >= 0 && r < thedim); 5087 x = diag[r] * rhs[r]; 5088 rhs[r] = 0; 5089 5090 if(x != 0) 5091 { 5092 c = corig[i]; 5093 assert(c >= 0 && c < thedim); 5094 vidx[n++] = c; 5095 vec[c] = x; 5096 val = &cval[cbeg[c]]; 5097 idx = &cidx[cbeg[c]]; 5098 j = clen[c]; 5099 5100 while(j-- > 0) 5101 { 5102 assert(*idx >= 0 && *idx < thedim); 5103 rhs[*idx++] -= x * (*val++); 5104 } 5105 } 5106 } 5107 5108 break; 5109 } 5110 } 5111 } 5112 5113 return n; 5114 } 5115 5116 5117 inline void CLUFactorRational::vSolveUrightNoNZ(Rational* vec, Rational* rhs, int* ridx, int rn) 5118 { 5119 int i, j, k, r, c; 5120 int* rorig, *corig; 5121 int* rperm; 5122 int* cidx, *clen, *cbeg; 5123 Rational x, y; 5124 5125 int* idx; 5126 Rational* val; 5127 5128 rorig = row.orig; 5129 corig = col.orig; 5130 rperm = row.perm; 5131 5132 cidx = u.col.idx; 5133 VectorRational& cval = u.col.val; 5134 clen = u.col.len; 5135 cbeg = u.col.start; 5136 5137 while(rn > 0) 5138 { 5139 if(rn > *ridx * verySparseFactor4rightRat) 5140 { 5141 /* continue with dense case */ 5142 for(i = *ridx; i >= 0; --i) 5143 { 5144 assert(i >= 0 && i < thedim); 5145 r = rorig[i]; 5146 assert(r >= 0 && r < thedim); 5147 x = diag[r] * rhs[r]; 5148 rhs[r] = 0; 5149 5150 if(x != 0) 5151 { 5152 c = corig[i]; 5153 vec[c] = x; 5154 val = &cval[cbeg[c]]; 5155 idx = &cidx[cbeg[c]]; 5156 j = clen[c]; 5157 5158 while(j-- > 0) 5159 { 5160 assert(*idx >= 0 && *idx < thedim); 5161 rhs[*idx++] -= x * (*val++); 5162 } 5163 } 5164 } 5165 5166 break; 5167 } 5168 5169 /* Find nonzero with highest permuted row index and setup i and r 5170 */ 5171 i = deQueueMaxRat(ridx, &rn); 5172 5173 assert(i >= 0 && i < thedim); 5174 5175 r = rorig[i]; 5176 5177 assert(r >= 0 && r < thedim); 5178 5179 x = diag[r] * rhs[r]; 5180 5181 rhs[r] = 0; 5182 5183 if(x != 0) 5184 { 5185 c = corig[i]; 5186 vec[c] = x; 5187 val = &cval[cbeg[c]]; 5188 idx = &cidx[cbeg[c]]; 5189 j = clen[c]; 5190 5191 while(j-- > 0) 5192 { 5193 k = *idx++; 5194 assert(k >= 0 && k < thedim); 5195 y = rhs[k]; 5196 5197 if(y == 0) 5198 { 5199 y = -x * (*val++); 5200 5201 if(y != 0) 5202 { 5203 rhs[k] = y; 5204 enQueueMaxRat(ridx, &rn, rperm[k]); 5205 } 5206 } 5207 else 5208 { 5209 y -= x * (*val++); 5210 // y += ( y == 0 ) ? SOPLEX_MARKER : 0; 5211 rhs[k] = y; 5212 } 5213 } 5214 } 5215 } 5216 } 5217 5218 5219 inline int CLUFactorRational::vSolveUright2(Rational* vec, int* vidx, Rational* rhs, int* ridx, 5220 int rn, 5221 Rational* vec2, Rational* rhs2, int* ridx2, int rn2) 5222 { 5223 int i, j, k, r, c, n; 5224 int* rorig, *corig; 5225 int* rperm; 5226 int* cidx, *clen, *cbeg; 5227 Rational x, y; 5228 Rational x2, y2; 5229 5230 int* idx; 5231 Rational* val; 5232 5233 rorig = row.orig; 5234 corig = col.orig; 5235 rperm = row.perm; 5236 5237 cidx = u.col.idx; 5238 VectorRational& cval = u.col.val; 5239 clen = u.col.len; 5240 cbeg = u.col.start; 5241 5242 n = 0; 5243 5244 while(rn + rn2 > 0) 5245 { 5246 /* Find nonzero with highest permuted row index and setup i and r 5247 */ 5248 if(rn <= 0) 5249 i = deQueueMaxRat(ridx2, &rn2); 5250 else if(rn2 <= 0) 5251 i = deQueueMaxRat(ridx, &rn); 5252 else if(*ridx2 > *ridx) 5253 i = deQueueMaxRat(ridx2, &rn2); 5254 else if(*ridx2 < *ridx) 5255 i = deQueueMaxRat(ridx, &rn); 5256 else 5257 { 5258 (void) deQueueMaxRat(ridx, &rn); 5259 i = deQueueMaxRat(ridx2, &rn2); 5260 } 5261 5262 assert(i >= 0 && i < thedim); 5263 5264 r = rorig[i]; 5265 assert(r >= 0 && r < thedim); 5266 5267 x = diag[r] * rhs[r]; 5268 x2 = diag[r] * rhs2[r]; 5269 rhs[r] = 0; 5270 rhs2[r] = 0; 5271 5272 if(x != 0) 5273 { 5274 c = corig[i]; 5275 vidx[n++] = c; 5276 vec[c] = x; 5277 vec2[c] = x2; 5278 val = &cval[cbeg[c]]; 5279 idx = &cidx[cbeg[c]]; 5280 j = clen[c]; 5281 5282 if(x2 != 0) 5283 { 5284 while(j-- > 0) 5285 { 5286 k = *idx++; 5287 assert(k >= 0 && k < thedim); 5288 y2 = rhs2[k]; 5289 5290 if(y2 == 0) 5291 { 5292 y2 = -x2 * (*val); 5293 5294 if(y2 != 0) 5295 { 5296 rhs2[k] = y2; 5297 enQueueMaxRat(ridx2, &rn2, rperm[k]); 5298 } 5299 } 5300 else 5301 { 5302 y2 -= x2 * (*val); 5303 // rhs2[k] = ( y2 != 0 ) ? y2 : SOPLEX_MARKER; 5304 rhs2[k] = y2; 5305 } 5306 5307 y = rhs[k]; 5308 5309 if(y == 0) 5310 { 5311 y = -x * (*val++); 5312 5313 if(y != 0) 5314 { 5315 rhs[k] = y; 5316 enQueueMaxRat(ridx, &rn, rperm[k]); 5317 } 5318 } 5319 else 5320 { 5321 y -= x * (*val++); 5322 // y += ( y == 0 ) ? SOPLEX_MARKER : 0; 5323 rhs[k] = y; 5324 } 5325 } 5326 } 5327 else 5328 { 5329 while(j-- > 0) 5330 { 5331 k = *idx++; 5332 assert(k >= 0 && k < thedim); 5333 y = rhs[k]; 5334 5335 if(y == 0) 5336 { 5337 y = -x * (*val++); 5338 5339 if(y != 0) 5340 { 5341 rhs[k] = y; 5342 enQueueMaxRat(ridx, &rn, rperm[k]); 5343 } 5344 } 5345 else 5346 { 5347 y -= x * (*val++); 5348 // y += ( y == 0 ) ? SOPLEX_MARKER : 0; 5349 rhs[k] = y; 5350 } 5351 } 5352 } 5353 } 5354 else if(x2 != 0) 5355 { 5356 c = corig[i]; 5357 assert(c >= 0 && c < thedim); 5358 vec2[c] = x2; 5359 val = &cval[cbeg[c]]; 5360 idx = &cidx[cbeg[c]]; 5361 j = clen[c]; 5362 5363 while(j-- > 0) 5364 { 5365 k = *idx++; 5366 assert(k >= 0 && k < thedim); 5367 y2 = rhs2[k]; 5368 5369 if(y2 == 0) 5370 { 5371 y2 = -x2 * (*val++); 5372 5373 if(y2 != 0) 5374 { 5375 rhs2[k] = y2; 5376 enQueueMaxRat(ridx2, &rn2, rperm[k]); 5377 } 5378 } 5379 else 5380 { 5381 y2 -= x2 * (*val++); 5382 // rhs2[k] = ( y2 != 0 ) ? y2 : SOPLEX_MARKER; 5383 rhs2[k] = y2; 5384 } 5385 } 5386 } 5387 5388 if(rn + rn2 > i * verySparseFactor4rightRat) 5389 { 5390 /* continue with dense case */ 5391 if(*ridx > *ridx2) 5392 i = *ridx; 5393 else 5394 i = *ridx2; 5395 5396 for(; i >= 0; --i) 5397 { 5398 assert(i < thedim); 5399 r = rorig[i]; 5400 assert(r >= 0 && r < thedim); 5401 x = diag[r] * rhs[r]; 5402 x2 = diag[r] * rhs2[r]; 5403 rhs[r] = 0; 5404 rhs2[r] = 0; 5405 5406 if(x2 != 0) 5407 { 5408 c = corig[i]; 5409 assert(c >= 0 && c < thedim); 5410 vec2[c] = x2; 5411 val = &cval[cbeg[c]]; 5412 idx = &cidx[cbeg[c]]; 5413 j = clen[c]; 5414 5415 if(x != 0) 5416 { 5417 vidx[n++] = c; 5418 vec[c] = x; 5419 5420 while(j-- > 0) 5421 { 5422 assert(*idx >= 0 && *idx < thedim); 5423 rhs[*idx] -= x * (*val); 5424 rhs2[*idx++] -= x2 * (*val++); 5425 } 5426 } 5427 else 5428 { 5429 while(j-- > 0) 5430 { 5431 assert(*idx >= 0 && *idx < thedim); 5432 rhs2[*idx++] -= x2 * (*val++); 5433 } 5434 } 5435 } 5436 else if(x != 0) 5437 { 5438 c = corig[i]; 5439 assert(c >= 0 && c < thedim); 5440 vidx[n++] = c; 5441 vec[c] = x; 5442 val = &cval[cbeg[c]]; 5443 idx = &cidx[cbeg[c]]; 5444 j = clen[c]; 5445 5446 while(j-- > 0) 5447 { 5448 assert(*idx >= 0 && *idx < thedim); 5449 rhs[*idx++] -= x * (*val++); 5450 } 5451 } 5452 } 5453 5454 break; 5455 } 5456 } 5457 5458 return n; 5459 } 5460 5461 inline int CLUFactorRational::vSolveUpdateRight(Rational* vec, int* ridx, int n) 5462 { 5463 int i, j, k; 5464 int end; 5465 Rational x, y; 5466 Rational* val; 5467 int* lrow, *lidx, *idx; 5468 int* lbeg; 5469 5470 assert(!l.updateType); /* no Forest-Tomlin Updates */ 5471 5472 VectorRational& lval = l.val; 5473 lidx = l.idx; 5474 lrow = l.row; 5475 lbeg = l.start; 5476 end = l.firstUnused; 5477 5478 for(i = l.firstUpdate; i < end; ++i) 5479 { 5480 assert(i >= 0 && i < thedim); 5481 x = vec[lrow[i]]; 5482 5483 if(x != 0) 5484 { 5485 k = lbeg[i]; 5486 assert(k >= 0 && k < l.val.dim()); 5487 idx = &(lidx[k]); 5488 val = &(lval[k]); 5489 5490 for(j = lbeg[i + 1]; j > k; --j) 5491 { 5492 int m = ridx[n] = *idx++; 5493 assert(m >= 0 && m < thedim); 5494 y = vec[m]; 5495 n += (y == 0) ? 1 : 0; 5496 y = y - x * (*val++); 5497 // vec[m] = ( y != 0 ) ? y : SOPLEX_MARKER; 5498 vec[m] = y; 5499 } 5500 } 5501 } 5502 5503 return n; 5504 } 5505 5506 inline void CLUFactorRational::vSolveUpdateRightNoNZ(Rational* vec) 5507 { 5508 int i, j, k; 5509 int end; 5510 Rational x; 5511 Rational* val; 5512 int* lrow, *lidx, *idx; 5513 int* lbeg; 5514 5515 assert(!l.updateType); /* no Forest-Tomlin Updates */ 5516 5517 VectorRational& lval = l.val; 5518 lidx = l.idx; 5519 lrow = l.row; 5520 lbeg = l.start; 5521 end = l.firstUnused; 5522 5523 for(i = l.firstUpdate; i < end; ++i) 5524 { 5525 assert(i >= 0 && i < thedim); 5526 5527 if((x = vec[lrow[i]]) != 0) 5528 { 5529 k = lbeg[i]; 5530 assert(k >= 0 && k < l.val.dim()); 5531 idx = &(lidx[k]); 5532 val = &(lval[k]); 5533 5534 for(j = lbeg[i + 1]; j > k; --j) 5535 { 5536 assert(*idx >= 0 && *idx < thedim); 5537 vec[*idx++] -= x * (*val++); 5538 } 5539 } 5540 } 5541 } 5542 5543 5544 inline int CLUFactorRational::vSolveRight4update(Rational* vec, 5545 int* idx, /* result */ 5546 Rational* rhs, int* ridx, int rn, /* rhs */ 5547 Rational* forest, int* forestNum, int* forestIdx) 5548 { 5549 rn = vSolveLright(rhs, ridx, rn); 5550 5551 /* turn index list into a heap 5552 */ 5553 5554 if(forest) 5555 { 5556 Rational x; 5557 int i, j, k; 5558 int* rperm; 5559 int* it = forestIdx; 5560 5561 rperm = row.perm; 5562 5563 for(i = j = 0; i < rn; ++i) 5564 { 5565 k = ridx[i]; 5566 assert(k >= 0 && k < thedim); 5567 x = rhs[k]; 5568 5569 if(x != 0) 5570 { 5571 enQueueMaxRat(ridx, &j, rperm[*it++ = k]); 5572 forest[k] = x; 5573 } 5574 else 5575 rhs[k] = 0; 5576 } 5577 5578 *forestNum = rn = j; 5579 } 5580 else 5581 { 5582 Rational x; 5583 int i, j, k; 5584 int* rperm; 5585 5586 rperm = row.perm; 5587 5588 for(i = j = 0; i < rn; ++i) 5589 { 5590 k = ridx[i]; 5591 assert(k >= 0 && k < thedim); 5592 x = rhs[k]; 5593 5594 if(x != 0) 5595 enQueueMaxRat(ridx, &j, rperm[k]); 5596 else 5597 rhs[k] = 0; 5598 } 5599 5600 rn = j; 5601 } 5602 5603 rn = vSolveUright(vec, idx, rhs, ridx, rn); 5604 5605 if(!l.updateType) /* no Forest-Tomlin Updates */ 5606 rn = vSolveUpdateRight(vec, idx, rn); 5607 5608 return rn; 5609 } 5610 5611 inline int CLUFactorRational::vSolveRight4update2(Rational* vec, 5612 int* idx, /* result1 */ 5613 Rational* rhs, int* ridx, int rn, /* rhs1 */ 5614 Rational* vec2, /* result2 */ 5615 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 5616 Rational* forest, int* forestNum, int* forestIdx) 5617 { 5618 vSolveLright2(rhs, ridx, &rn, rhs2, ridx2, &rn2); 5619 5620 /* turn index list into a heap 5621 */ 5622 5623 if(forest) 5624 { 5625 Rational x; 5626 int i, j, k; 5627 int* rperm; 5628 int* it = forestIdx; 5629 5630 rperm = row.perm; 5631 5632 for(i = j = 0; i < rn; ++i) 5633 { 5634 k = ridx[i]; 5635 assert(k >= 0 && k < thedim); 5636 x = rhs[k]; 5637 5638 if(x != 0) 5639 { 5640 enQueueMaxRat(ridx, &j, rperm[*it++ = k]); 5641 forest[k] = x; 5642 } 5643 else 5644 rhs[k] = 0; 5645 } 5646 5647 *forestNum = rn = j; 5648 } 5649 else 5650 { 5651 Rational x; 5652 int i, j, k; 5653 int* rperm; 5654 5655 rperm = row.perm; 5656 5657 for(i = j = 0; i < rn; ++i) 5658 { 5659 k = ridx[i]; 5660 assert(k >= 0 && k < thedim); 5661 x = rhs[k]; 5662 5663 if(x != 0) 5664 enQueueMaxRat(ridx, &j, rperm[k]); 5665 else 5666 rhs[k] = 0; 5667 } 5668 5669 rn = j; 5670 } 5671 5672 if(rn2 > thedim * verySparseFactor4rightRat) 5673 { 5674 ridx2[0] = thedim - 1; 5675 /* ridx2[1] = thedim - 2; */ 5676 } 5677 else 5678 { 5679 Rational x; 5680 /* Rational maxabs; */ 5681 int i, j, k; 5682 int* rperm; 5683 5684 /* maxabs = 1; */ 5685 rperm = row.perm; 5686 5687 for(i = j = 0; i < rn2; ++i) 5688 { 5689 k = ridx2[i]; 5690 assert(k >= 0 && k < thedim); 5691 x = rhs2[k]; 5692 5693 if(x == 0) 5694 { 5695 /* maxabs = (maxabs < -x) ? -x : maxabs; */ 5696 enQueueMaxRat(ridx2, &j, rperm[k]); 5697 } 5698 else 5699 rhs2[k] = 0; 5700 } 5701 5702 rn2 = j; 5703 5704 } 5705 5706 rn = vSolveUright(vec, idx, rhs, ridx, rn); 5707 5708 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2); 5709 5710 /* 5711 * rn = vSolveUright2(vec, idx, rhs, ridx, rn, vec2, rhs2, ridx2, rn2 ); 5712 */ 5713 5714 if(!l.updateType) /* no Forest-Tomlin Updates */ 5715 { 5716 rn = vSolveUpdateRight(vec, idx, rn); 5717 vSolveUpdateRightNoNZ(vec2); 5718 } 5719 5720 return rn; 5721 } 5722 5723 inline int CLUFactorRational::vSolveRight4update3(Rational* vec, 5724 int* idx, /* result1 */ 5725 Rational* rhs, int* ridx, int rn, /* rhs1 */ 5726 Rational* vec2, /* result2 */ 5727 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 5728 Rational* vec3, /* result3 */ 5729 Rational* rhs3, int* ridx3, int rn3, /* rhs3 */ 5730 Rational* forest, int* forestNum, int* forestIdx) 5731 { 5732 5733 vSolveLright3(rhs, ridx, &rn, rhs2, ridx2, &rn2, rhs3, ridx3, &rn3); 5734 assert(rn >= 0 && rn <= thedim); 5735 assert(rn2 >= 0 && rn2 <= thedim); 5736 assert(rn3 >= 0 && rn3 <= thedim); 5737 5738 /* turn index list into a heap 5739 */ 5740 5741 if(forest) 5742 { 5743 Rational x; 5744 int i, j, k; 5745 int* rperm; 5746 int* it = forestIdx; 5747 5748 rperm = row.perm; 5749 5750 for(i = j = 0; i < rn; ++i) 5751 { 5752 k = ridx[i]; 5753 assert(k >= 0 && k < thedim); 5754 x = rhs[k]; 5755 5756 if(x != 0) 5757 { 5758 enQueueMaxRat(ridx, &j, rperm[*it++ = k]); 5759 forest[k] = x; 5760 } 5761 else 5762 rhs[k] = 0; 5763 } 5764 5765 *forestNum = rn = j; 5766 } 5767 else 5768 { 5769 Rational x; 5770 int i, j, k; 5771 int* rperm; 5772 5773 rperm = row.perm; 5774 5775 for(i = j = 0; i < rn; ++i) 5776 { 5777 k = ridx[i]; 5778 assert(k >= 0 && k < thedim); 5779 x = rhs[k]; 5780 5781 if(x != 0) 5782 enQueueMaxRat(ridx, &j, rperm[k]); 5783 else 5784 rhs[k] = 0; 5785 } 5786 5787 rn = j; 5788 } 5789 5790 if(rn2 > thedim * verySparseFactor4rightRat) 5791 { 5792 ridx2[0] = thedim - 1; 5793 } 5794 else 5795 { 5796 Rational x; 5797 int i, j, k; 5798 int* rperm; 5799 5800 rperm = row.perm; 5801 5802 for(i = j = 0; i < rn2; ++i) 5803 { 5804 k = ridx2[i]; 5805 assert(k >= 0 && k < thedim); 5806 x = rhs2[k]; 5807 5808 if(x == 0) 5809 { 5810 enQueueMaxRat(ridx2, &j, rperm[k]); 5811 } 5812 else 5813 rhs2[k] = 0; 5814 } 5815 5816 rn2 = j; 5817 } 5818 5819 if(rn3 > thedim * verySparseFactor4rightRat) 5820 { 5821 ridx3[0] = thedim - 1; 5822 } 5823 else 5824 { 5825 Rational x; 5826 int i, j, k; 5827 int* rperm; 5828 5829 rperm = row.perm; 5830 5831 for(i = j = 0; i < rn3; ++i) 5832 { 5833 k = ridx3[i]; 5834 assert(k >= 0 && k < thedim); 5835 x = rhs3[k]; 5836 5837 if(x == 0) 5838 { 5839 enQueueMaxRat(ridx3, &j, rperm[k]); 5840 } 5841 else 5842 rhs3[k] = 0; 5843 } 5844 5845 rn3 = j; 5846 } 5847 5848 rn = vSolveUright(vec, idx, rhs, ridx, rn); 5849 5850 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2); 5851 vSolveUrightNoNZ(vec3, rhs3, ridx3, rn3); 5852 5853 if(!l.updateType) /* no Forest-Tomlin Updates */ 5854 { 5855 rn = vSolveUpdateRight(vec, idx, rn); 5856 vSolveUpdateRightNoNZ(vec2); 5857 vSolveUpdateRightNoNZ(vec3); 5858 } 5859 5860 return rn; 5861 } 5862 5863 inline void CLUFactorRational::vSolveRightNoNZ(Rational* 5864 vec2, /* result2 */ 5865 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */ 5866 { 5867 rn2 = vSolveLright(rhs2, ridx2, rn2); 5868 assert(rn2 >= 0 && rn2 <= thedim); 5869 5870 if(rn2 > thedim * verySparseFactor4rightRat) 5871 { 5872 *ridx2 = thedim - 1; 5873 } 5874 else 5875 { 5876 Rational x; 5877 /* Rational maxabs; */ 5878 int i, j, k; 5879 int* rperm; 5880 5881 /* maxabs = 1; */ 5882 rperm = row.perm; 5883 5884 for(i = j = 0; i < rn2; ++i) 5885 { 5886 k = ridx2[i]; 5887 assert(k >= 0 && k < thedim); 5888 x = rhs2[k]; 5889 5890 if(x == 0) 5891 { 5892 /* maxabs = (maxabs < -x) ? -x : maxabs; */ 5893 enQueueMaxRat(ridx2, &j, rperm[k]); 5894 } 5895 else 5896 rhs2[k] = 0; 5897 } 5898 5899 rn2 = j; 5900 } 5901 5902 vSolveUrightNoNZ(vec2, rhs2, ridx2, rn2); 5903 5904 if(!l.updateType) /* no Forest-Tomlin Updates */ 5905 vSolveUpdateRightNoNZ(vec2); 5906 } 5907 5908 inline int CLUFactorRational::vSolveLeft(Rational* vec, int* idx, /* result */ 5909 Rational* rhs, int* ridx, int rn) /* rhs */ 5910 { 5911 5912 if(!l.updateType) /* no Forest-Tomlin Updates */ 5913 { 5914 rn = solveUpdateLeft(rhs, ridx, rn); 5915 rn = solveUleft(vec, idx, rhs, ridx, rn); 5916 } 5917 else 5918 { 5919 rn = solveUleft(vec, idx, rhs, ridx, rn); 5920 rn = solveLleftForest(vec, idx, rn); 5921 } 5922 5923 return solveLleft(vec, idx, rn); 5924 } 5925 5926 inline int CLUFactorRational::vSolveLeft2(Rational* vec, 5927 int* idx, /* result */ 5928 Rational* rhs, int* ridx, int rn, /* rhs */ 5929 Rational* vec2, /* result2 */ 5930 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */ 5931 { 5932 5933 if(!l.updateType) /* no Forest-Tomlin Updates */ 5934 { 5935 rn = solveUpdateLeft(rhs, ridx, rn); 5936 rn = solveUleft(vec, idx, rhs, ridx, rn); 5937 rn2 = solveUpdateLeft(rhs2, ridx2, rn2); 5938 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 5939 } 5940 else 5941 { 5942 rn = solveUleft(vec, idx, rhs, ridx, rn); 5943 rn = solveLleftForest(vec, idx, rn); 5944 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 5945 solveLleftForestNoNZ(vec2); 5946 } 5947 5948 rn = solveLleft(vec, idx, rn); 5949 5950 solveLleftNoNZ(vec2); 5951 5952 return rn; 5953 } 5954 5955 inline int CLUFactorRational::vSolveLeft3(Rational* vec, 5956 int* idx, /* result */ 5957 Rational* rhs, int* ridx, int rn, /* rhs */ 5958 Rational* vec2, /* result2 */ 5959 Rational* rhs2, int* ridx2, int rn2, /* rhs2 */ 5960 Rational* vec3, /* result3 */ 5961 Rational* rhs3, int* ridx3, int rn3) /* rhs3 */ 5962 { 5963 5964 if(!l.updateType) /* no Forest-Tomlin Updates */ 5965 { 5966 rn = solveUpdateLeft(rhs, ridx, rn); 5967 rn = solveUleft(vec, idx, rhs, ridx, rn); 5968 rn2 = solveUpdateLeft(rhs2, ridx2, rn2); 5969 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 5970 rn3 = solveUpdateLeft(rhs3, ridx3, rn3); 5971 solveUleftNoNZ(vec3, rhs3, ridx3, rn3); 5972 } 5973 else 5974 { 5975 rn = solveUleft(vec, idx, rhs, ridx, rn); 5976 rn = solveLleftForest(vec, idx, rn); 5977 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 5978 solveLleftForestNoNZ(vec2); 5979 solveUleftNoNZ(vec3, rhs3, ridx3, rn3); 5980 solveLleftForestNoNZ(vec3); 5981 } 5982 5983 rn = solveLleft(vec, idx, rn); 5984 5985 solveLleftNoNZ(vec2); 5986 solveLleftNoNZ(vec3); 5987 5988 return rn; 5989 } 5990 5991 inline void CLUFactorRational::vSolveLeftNoNZ(Rational* 5992 vec2, /* result2 */ 5993 Rational* rhs2, int* ridx2, int rn2) /* rhs2 */ 5994 { 5995 5996 if(!l.updateType) /* no Forest-Tomlin Updates */ 5997 { 5998 rn2 = solveUpdateLeft(rhs2, ridx2, rn2); 5999 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 6000 } 6001 else 6002 { 6003 solveUleftNoNZ(vec2, rhs2, ridx2, rn2); 6004 solveLleftForestNoNZ(vec2); 6005 } 6006 6007 solveLleftNoNZ(vec2); 6008 } 6009 } // namespace soplex 6010