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