1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 2 /* */ 3 /* This file is part of the class library */ 4 /* SoPlex --- the Sequential object-oriented simPlex. */ 5 /* */ 6 /* Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB) */ 7 /* */ 8 /* Licensed under the Apache License, Version 2.0 (the "License"); */ 9 /* you may not use this file except in compliance with the License. */ 10 /* You may obtain a copy of the License at */ 11 /* */ 12 /* http://www.apache.org/licenses/LICENSE-2.0 */ 13 /* */ 14 /* Unless required by applicable law or agreed to in writing, software */ 15 /* distributed under the License is distributed on an "AS IS" BASIS, */ 16 /* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */ 17 /* See the License for the specific language governing permissions and */ 18 /* limitations under the License. */ 19 /* */ 20 /* You should have received a copy of the Apache-2.0 license */ 21 /* along with SoPlex; see the file LICENSE. If not email to soplex@zib.de. */ 22 /* */ 23 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ 24 25 /**@file basevectors.h 26 * @brief Collection of dense, sparse, and semi-sparse vectors. 27 */ 28 #ifndef _BASEVECTORS_H_ 29 #define _BASEVECTORS_H_ 30 31 /* undefine SOPLEX_DEBUG flag from including files; if SOPLEX_DEBUG should be defined in this file, do so below */ 32 #ifdef SOPLEX_DEBUG 33 #define SOPLEX_DEBUG_BASEVECTORS 34 #undef SOPLEX_DEBUG 35 #endif 36 37 #include "soplex/spxdefines.h" 38 #include "soplex/rational.h" 39 #include "soplex/vectorbase.h" 40 #include "soplex/ssvectorbase.h" 41 #include "soplex/svectorbase.h" 42 #include "soplex/dsvectorbase.h" 43 #include "soplex/unitvectorbase.h" 44 #include "soplex/svsetbase.h" 45 #include "soplex/timer.h" 46 47 #define SOPLEX_VECTOR_MARKER 1e-100 48 49 namespace soplex 50 { 51 52 // --------------------------------------------------------------------------------------------------------------------- 53 // Methods of VectorBase 54 // --------------------------------------------------------------------------------------------------------------------- 55 56 /// Assignment operator. 57 /** Assigning an SVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p vec. 58 * This is different in method assign(). 59 */ 60 61 62 63 template < class R > 64 template < class S > 65 inline 66 VectorBase<R>& VectorBase<R>::operator=(const SVectorBase<S>& vec) 67 { 68 clear(); 69 70 for(int i = 0; i < vec.size(); ++i) 71 { 72 assert(vec.index(i) < dim()); 73 val[vec.index(i)] = vec.value(i); 74 } 75 76 return *this; 77 } 78 79 80 81 /// Assign values of \p vec. 82 /** Assigns all nonzeros of \p vec to the vector. All other values remain unchanged. */ 83 template < class R > 84 template < class S > 85 inline 86 VectorBase<R>& VectorBase<R>::assign(const SVectorBase<S>& vec) 87 { 88 for(int i = vec.size() - 1; i >= 0; --i) 89 { 90 assert(vec.index(i) < dim()); 91 val[vec.index(i)] = vec.value(i); 92 } 93 94 return *this; 95 } 96 97 98 99 /// Assignment operator. 100 /** Assigning an SSVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p vec. 101 * This is different in method assign(). 102 */ 103 template < class R > 104 template < class S > 105 inline 106 VectorBase<R>& VectorBase<R>::operator=(const SSVectorBase<S>& vec) 107 { 108 if(vec.isSetup()) 109 { 110 clear(); 111 assign(vec); 112 } 113 else 114 operator=(static_cast<const VectorBase<R>&>(vec)); 115 116 return *this; 117 } 118 119 120 121 /// Assign values of \p vec. 122 /** Assigns all nonzeros of \p vec to the vector. All other values remain unchanged. */ 123 template < class R > 124 template < class S > 125 inline 126 VectorBase<R>& VectorBase<R>::assign(const SSVectorBase<S>& vec) 127 { 128 assert(vec.dim() <= dim()); 129 130 if(vec.isSetup()) 131 { 132 const int* idx = vec.indexMem(); 133 134 for(int i = vec.size() - 1; i >= 0; --i) 135 { 136 val[*idx] = vec.val[*idx]; 137 idx++; 138 } 139 } 140 else 141 operator=(static_cast<const VectorBase<R>&>(vec)); 142 143 return *this; 144 } 145 146 147 148 /// Addition. 149 template < class R > 150 template < class S > 151 inline 152 VectorBase<R>& VectorBase<R>::operator+=(const SVectorBase<S>& vec) 153 { 154 for(int i = vec.size() - 1; i >= 0; --i) 155 { 156 assert(vec.index(i) >= 0); 157 assert(vec.index(i) < dim()); 158 val[vec.index(i)] += vec.value(i); 159 } 160 161 return *this; 162 } 163 164 165 166 /// Addition. 167 template < class R > 168 template < class S > 169 inline 170 VectorBase<R>& VectorBase<R>::operator+=(const SSVectorBase<S>& vec) 171 { 172 assert(dim() == vec.dim()); 173 174 if(vec.isSetup()) 175 { 176 for(int i = vec.size() - 1; i >= 0 ; --i) 177 val[vec.index(i)] += vec.value(i); 178 } 179 else 180 { 181 for(int i = dim() - 1; i >= 0; --i) 182 val[i] += vec[i]; 183 } 184 185 return *this; 186 } 187 188 189 190 /// Subtraction. 191 template < class R > 192 template < class S > 193 inline 194 VectorBase<R>& VectorBase<R>::operator-=(const SVectorBase<S>& vec) 195 { 196 for(int i = vec.size() - 1; i >= 0; --i) 197 { 198 assert(vec.index(i) >= 0); 199 assert(vec.index(i) < dim()); 200 val[vec.index(i)] -= vec.value(i); 201 } 202 203 return *this; 204 } 205 206 207 208 /// Subtraction. 209 template < class R > 210 template < class S > 211 inline 212 VectorBase<R>& VectorBase<R>::operator-=(const SSVectorBase<S>& vec) 213 { 214 assert(dim() == vec.dim()); 215 216 if(vec.isSetup()) 217 { 218 for(int i = vec.size() - 1; i >= 0; --i) 219 val[vec.index(i)] -= vec.value(i); 220 } 221 else 222 { 223 for(int i = dim() - 1; i >= 0; --i) 224 val[i] -= vec[i]; 225 } 226 227 return *this; 228 } 229 230 231 232 /// Inner product. 233 template < class R > 234 inline 235 R VectorBase<R>::operator*(const SVectorBase<R>& vec) const 236 { 237 assert(dim() >= vec.dim()); 238 239 StableSum<R> x; 240 241 for(int i = vec.size() - 1; i >= 0; --i) 242 x += val[vec.index(i)] * vec.value(i); 243 244 return x; 245 } 246 247 248 249 /// Inner product. 250 template < class R > 251 inline 252 R VectorBase<R>::operator*(const SSVectorBase<R>& vec) const 253 { 254 assert(dim() == vec.dim()); 255 256 if(vec.isSetup()) 257 { 258 const int* idx = vec.indexMem(); 259 260 StableSum<R> x; 261 262 for(int i = vec.size() - 1; i >= 0; --i) 263 { 264 x += val[*idx] * vec.val[*idx]; 265 idx++; 266 } 267 268 return x; 269 } 270 else 271 return operator*(static_cast<const VectorBase<R>&>(vec)); 272 } 273 274 275 276 /// Addition of scaled vector. 277 template < class R > 278 template < class S, class T > 279 inline 280 VectorBase<R>& VectorBase<R>::multAdd(const S& x, const SVectorBase<T>& vec) 281 { 282 for(int i = vec.size() - 1; i >= 0; --i) 283 { 284 assert(vec.index(i) < dim()); 285 val[vec.index(i)] += x * vec.value(i); 286 } 287 288 return *this; 289 } 290 291 292 293 /// Subtraction of scaled vector. 294 template < class R > 295 template < class S, class T > 296 inline 297 VectorBase<R>& VectorBase<R>::multSub(const S& x, const SVectorBase<T>& vec) 298 { 299 for(int i = vec.size() - 1; i >= 0; --i) 300 { 301 assert(vec.index(i) < dim()); 302 val[vec.index(i)] -= x * vec.value(i); 303 } 304 305 return *this; 306 } 307 308 309 310 /// Addition of scaled vector. 311 template < class R > 312 template < class S, class T > 313 inline 314 VectorBase<R>& VectorBase<R>::multAdd(const S& x, const SSVectorBase<T>& vec) 315 { 316 assert(vec.dim() <= dim()); 317 318 if(vec.isSetup()) 319 { 320 const int* idx = vec.indexMem(); 321 322 for(int i = vec.size() - 1; i >= 0; --i) 323 val[idx[i]] += x * vec[idx[i]]; 324 } 325 else 326 { 327 assert(vec.dim() == dim()); 328 329 for(int i = dim() - 1; i >= 0; --i) 330 val[i] += x * vec.val[i]; 331 } 332 333 return *this; 334 } 335 336 337 338 339 340 // --------------------------------------------------------------------------------------------------------------------- 341 // Methods of SSVectorBase 342 // --------------------------------------------------------------------------------------------------------------------- 343 344 345 346 /// Addition. 347 template < class R > 348 template < class S > 349 inline 350 SSVectorBase<R>& SSVectorBase<R>::operator+=(const SVectorBase<S>& vec) 351 { 352 VectorBase<R>::operator+=(vec); 353 354 if(isSetup()) 355 { 356 setupStatus = false; 357 setup(); 358 } 359 360 return *this; 361 } 362 363 364 365 /// Subtraction. 366 template < class R > 367 template < class S > 368 inline 369 SSVectorBase<R>& SSVectorBase<R>::operator-=(const SVectorBase<S>& vec) 370 { 371 VectorBase<R>::operator-=(vec); 372 373 if(isSetup()) 374 { 375 setupStatus = false; 376 setup(); 377 } 378 379 return *this; 380 } 381 382 383 384 /// Addition of a scaled vector. 385 ///@todo SSVectorBase::multAdd() should be rewritten without pointer arithmetic. 386 template < class R > 387 template < class S, class T > 388 inline 389 SSVectorBase<R>& SSVectorBase<R>::multAdd(S xx, const SVectorBase<T>& vec) 390 { 391 if(isSetup()) 392 { 393 R* v = VectorBase<R>::val.data(); 394 R x; 395 bool adjust = false; 396 int j; 397 398 for(int i = vec.size() - 1; i >= 0; --i) 399 { 400 j = vec.index(i); 401 402 if(v[j] != 0) 403 { 404 x = v[j] + xx * vec.value(i); 405 406 if(isNotZero(x, this->tolerances()->epsilon())) 407 v[j] = x; 408 else 409 { 410 adjust = true; 411 v[j] = SOPLEX_VECTOR_MARKER; 412 } 413 } 414 else 415 { 416 x = xx * vec.value(i); 417 418 if(isNotZero(x, this->tolerances()->epsilon())) 419 { 420 v[j] = x; 421 addIdx(j); 422 } 423 } 424 } 425 426 if(adjust) 427 { 428 int* iptr = idx; 429 int* iiptr = idx; 430 int* endptr = idx + num; 431 432 for(; iptr < endptr; ++iptr) 433 { 434 x = v[*iptr]; 435 436 if(isNotZero(x, this->tolerances()->epsilon())) 437 *iiptr++ = *iptr; 438 else 439 v[*iptr] = 0; 440 } 441 442 num = int(iiptr - idx); 443 } 444 } 445 else 446 VectorBase<R>::multAdd(xx, vec); 447 448 assert(isConsistent()); 449 450 return *this; 451 } 452 453 454 /// Assigns pair wise vector product of setup x and setup y to SSVectorBase. 455 template < class R > 456 template < class S, class T > 457 inline 458 SSVectorBase<R>& SSVectorBase<R>::assignPWproduct4setup(const SSVectorBase<S>& x, 459 const SSVectorBase<T>& y) 460 { 461 assert(dim() == x.dim()); 462 assert(x.dim() == y.dim()); 463 assert(x.isSetup()); 464 assert(y.isSetup()); 465 466 clear(); 467 setupStatus = false; 468 469 int i = 0; 470 int j = 0; 471 int n = x.size() - 1; 472 int m = y.size() - 1; 473 474 /* both x and y non-zero vectors? */ 475 if(m >= 0 && n >= 0) 476 { 477 int xi = x.index(i); 478 int yj = y.index(j); 479 480 while(i < n && j < m) 481 { 482 if(xi == yj) 483 { 484 VectorBase<R>::val[xi] = R(x.val[xi]) * R(y.val[xi]); 485 xi = x.index(++i); 486 yj = y.index(++j); 487 } 488 else if(xi < yj) 489 xi = x.index(++i); 490 else 491 yj = y.index(++j); 492 } 493 494 /* check (possible) remaining indices */ 495 496 while(i < n && xi != yj) 497 xi = x.index(++i); 498 499 while(j < m && xi != yj) 500 yj = y.index(++j); 501 502 if(xi == yj) 503 VectorBase<R>::val[xi] = R(x.val[xi]) * R(y.val[xi]); 504 } 505 506 setup(); 507 508 assert(isConsistent()); 509 510 return *this; 511 } 512 513 514 515 /// Assigns \f$x^T \cdot A\f$ to SSVectorBase. 516 template < class R > 517 template < class S, class T > 518 inline 519 SSVectorBase<R>& SSVectorBase<R>::assign2product(const SSVectorBase<S>& x, const SVSetBase<T>& A) 520 { 521 assert(A.num() == dim()); 522 523 R y; 524 525 clear(); 526 527 for(int i = dim() - 1; i >= 0; --i) 528 { 529 y = A[i] * x; 530 531 if(isNotZero(y, this->tolerances()->epsilon())) 532 { 533 VectorBase<R>::val[i] = y; 534 IdxSet::addIdx(i); 535 } 536 } 537 538 assert(isConsistent()); 539 540 return *this; 541 } 542 543 544 545 /// Assigns SSVectorBase to \f$A \cdot x\f$ for a setup \p x. 546 #define SOPLEX_SHORTPRODUCT_FACTOR 0.5 547 template < class R > 548 template < class S, class T > 549 inline 550 SSVectorBase<R>& SSVectorBase<R>::assign2product4setup(const SVSetBase<S>& A, 551 const SSVectorBase<T>& x, 552 Timer* timeSparse, Timer* timeFull, 553 int& nCallsSparse, int& nCallsFull 554 ) 555 { 556 assert(A.num() == x.dim()); 557 assert(x.isSetup()); 558 clear(); 559 560 if(x.size() == 1) 561 { 562 if(timeSparse != 0) 563 timeSparse->start(); 564 565 assign2product1(A, x); 566 setupStatus = true; 567 568 if(timeSparse != 0) 569 timeSparse->stop(); 570 571 ++nCallsSparse; 572 } 573 else if(isSetup() 574 && (double(x.size()) * A.memSize() <= SOPLEX_SHORTPRODUCT_FACTOR * dim() * A.num())) 575 { 576 if(timeSparse != 0) 577 timeSparse->start(); 578 579 assign2productShort(A, x); 580 setupStatus = true; 581 582 if(timeSparse != 0) 583 timeSparse->stop(); 584 585 ++nCallsSparse; 586 } 587 else 588 { 589 if(timeFull != 0) 590 timeFull->start(); 591 592 assign2productFull(A, x); 593 setupStatus = false; 594 595 if(timeFull != 0) 596 timeFull->stop(); 597 598 ++nCallsFull; 599 } 600 601 assert(isConsistent()); 602 603 return *this; 604 } 605 606 607 608 /// Assignment helper. 609 template < class R > 610 template < class S, class T > 611 inline 612 SSVectorBase<R>& SSVectorBase<R>::assign2product1(const SVSetBase<S>& A, const SSVectorBase<T>& x) 613 { 614 assert(x.isSetup()); 615 assert(x.size() == 1); 616 617 // get the nonzero value of x and the corresponding vector in A: 618 const int nzidx = x.idx[0]; 619 const T nzval = x.val[nzidx]; 620 const SVectorBase<S>& Ai = A[nzidx]; 621 622 // compute A[nzidx] * nzval: 623 if(isZero(nzval, this->tolerances()->epsilon()) || Ai.size() == 0) 624 clear(); // this := zero vector 625 else 626 { 627 num = Ai.size(); 628 629 for(int j = num - 1; j >= 0; --j) 630 { 631 const Nonzero<S>& Aij = Ai.element(j); 632 idx[j] = Aij.idx; 633 VectorBase<R>::val[Aij.idx] = nzval * Aij.val; 634 } 635 } 636 637 assert(isConsistent()); 638 639 return *this; 640 } 641 642 643 644 /// Assignment helper. 645 template < class R > 646 template < class S, class T > 647 inline 648 SSVectorBase<R>& SSVectorBase<R>::assign2productShort(const SVSetBase<S>& A, 649 const SSVectorBase<T>& x) 650 { 651 assert(x.isSetup()); 652 653 if(x.size() == 0) // x can be setup but have size 0 => this := zero vector 654 { 655 clear(); 656 return *this; 657 } 658 659 // compute x[0] * A[0] 660 int curidx = x.idx[0]; 661 const T x0 = x.val[curidx]; 662 const SVectorBase<S>& A0 = A[curidx]; 663 int nonzero_idx = 0; 664 int xsize = x.size(); 665 int Aisize; 666 667 num = A0.size(); 668 669 if(isZero(x0, this->tolerances()->epsilon()) || num == 0) 670 { 671 // A[0] == 0 or x[0] == 0 => this := zero vector 672 clear(); 673 } 674 else 675 { 676 for(int j = 0; j < num; ++j) 677 { 678 const Nonzero<S>& elt = A0.element(j); 679 const R product = x0 * elt.val; 680 681 // store the value in any case 682 idx[nonzero_idx] = elt.idx; 683 VectorBase<R>::val[elt.idx] = product; 684 685 // count only non-zero values; not 'isNotZero(product, epsilon)' 686 if(product != 0) 687 ++nonzero_idx; 688 } 689 } 690 691 // Compute the other x[i] * A[i] and add them to the existing vector. 692 for(int i = 1; i < xsize; ++i) 693 { 694 curidx = x.idx[i]; 695 const T xi = x.val[curidx]; 696 const SVectorBase<S>& Ai = A[curidx]; 697 698 // If A[i] == 0 or x[i] == 0, do nothing. 699 Aisize = Ai.size(); 700 701 if(isNotZero(xi, this->tolerances()->epsilon()) || Aisize == 0) 702 { 703 // Compute x[i] * A[i] and add it to the existing vector. 704 for(int j = 0; j < Aisize; ++j) 705 { 706 const Nonzero<S>& elt = Ai.element(j); 707 idx[nonzero_idx] = elt.idx; 708 R oldval = VectorBase<R>::val[elt.idx]; 709 710 // An old value of exactly 0 means the position is still unused. 711 // It will be used now (either by a new nonzero or by a SOPLEX_VECTOR_MARKER), 712 // so increase the counter. If oldval != 0, we just 713 // change an existing NZ-element, so don't increase the counter. 714 if(oldval == 0) 715 ++nonzero_idx; 716 717 // Add the current product x[i] * A[i][j]; if oldval was 718 // SOPLEX_VECTOR_MARKER before, it does not hurt because SOPLEX_VECTOR_MARKER is really small. 719 oldval += xi * elt.val; 720 721 // If the new value is exactly 0, mark the index as used 722 // by setting a value which is nearly 0; otherwise, store 723 // the value. Values below epsilon will be removed later. 724 if(oldval == 0) 725 VectorBase<R>::val[elt.idx] = SOPLEX_VECTOR_MARKER; 726 else 727 VectorBase<R>::val[elt.idx] = oldval; 728 } 729 } 730 } 731 732 // Clean up by shifting all nonzeros (w.r.t. epsilon) to the front of idx, 733 // zeroing all values which are nearly 0, and setting #num# appropriately. 734 int nz_counter = 0; 735 736 for(int i = 0; i < nonzero_idx; ++i) 737 { 738 curidx = idx[i]; 739 740 if(isZero(VectorBase<R>::val[curidx], this->tolerances()->epsilon())) 741 VectorBase<R>::val[curidx] = 0; 742 else 743 { 744 idx[nz_counter] = curidx; 745 ++nz_counter; 746 } 747 748 num = nz_counter; 749 } 750 751 assert(isConsistent()); 752 753 return *this; 754 } 755 756 757 758 /// Assignment helper. 759 template < class R > 760 template < class S, class T > 761 inline 762 SSVectorBase<R>& SSVectorBase<R>::assign2productFull(const SVSetBase<S>& A, 763 const SSVectorBase<T>& x) 764 { 765 assert(x.isSetup()); 766 767 if(x.size() == 0) // x can be setup but have size 0 => this := zero vector 768 { 769 clear(); 770 return *this; 771 } 772 773 bool A_is_zero = true; 774 int xsize = x.size(); 775 int Aisize; 776 777 for(int i = 0; i < xsize; ++i) 778 { 779 const int curidx = x.idx[i]; 780 const T xi = x.val[curidx]; 781 const SVectorBase<S>& Ai = A[curidx]; 782 Aisize = Ai.size(); 783 784 if(A_is_zero && Aisize > 0) 785 A_is_zero = false; 786 787 for(int j = 0; j < Aisize; ++j) 788 { 789 const Nonzero<S>& elt = Ai.element(j); 790 VectorBase<R>::val[elt.idx] += xi * elt.val; 791 } 792 } 793 794 if(A_is_zero) 795 clear(); // case x != 0 but A == 0 796 797 return *this; 798 } 799 800 801 802 /// Assigns SSVectorBase to \f$A \cdot x\f$ thereby setting up \p x. 803 template < class R > 804 template < class S, class T > 805 inline 806 SSVectorBase<R>& SSVectorBase<R>::assign2productAndSetup(const SVSetBase<S>& A, SSVectorBase<T>& x) 807 { 808 assert(!x.isSetup()); 809 810 if(x.dim() == 0) 811 { 812 // x == 0 => this := zero vector 813 clear(); 814 x.num = 0; 815 } 816 else 817 { 818 // x is not setup, so walk through its value vector 819 int nzcount = 0; 820 int end = x.dim(); 821 822 for(int i = 0; i < end; ++i) 823 { 824 // advance to the next element != 0 825 T& xval = x.val[i]; 826 827 if(xval != 0) 828 { 829 // If x[i] is really nonzero, compute A[i] * x[i] and adapt x.idx, 830 // otherwise set x[i] to 0. 831 if(isNotZero(xval, this->tolerances()->epsilon())) 832 { 833 const SVectorBase<S>& Ai = A[i]; 834 x.idx[ nzcount++ ] = i; 835 836 for(int j = Ai.size() - 1; j >= 0; --j) 837 { 838 const Nonzero<S>& elt = Ai.element(j); 839 VectorBase<R>::val[elt.idx] += xval * elt.val; 840 } 841 } 842 else 843 xval = 0; 844 } 845 } 846 847 x.num = nzcount; 848 setupStatus = false; 849 } 850 851 x.setupStatus = true; 852 853 assert(isConsistent()); 854 855 return *this; 856 } 857 858 859 860 /// Assigns only the elements of \p rhs. 861 template < class R > 862 template < class S > 863 inline 864 SSVectorBase<R>& SSVectorBase<R>::assign(const SVectorBase<S>& rhs) 865 { 866 assert(rhs.dim() <= VectorBase<R>::dim()); 867 868 int s = rhs.size(); 869 num = 0; 870 871 for(int i = 0; i < s; ++i) 872 { 873 int k = rhs.index(i); 874 S v = rhs.value(i); 875 876 if(isZero(v, this->tolerances()->epsilon())) 877 VectorBase<R>::val[k] = 0; 878 else 879 { 880 VectorBase<R>::val[k] = v; 881 idx[num++] = k; 882 } 883 } 884 885 setupStatus = true; 886 887 assert(isConsistent()); 888 889 return *this; 890 } 891 892 893 894 /// Assigns only the elements of \p rhs. 895 template < > 896 template < > 897 inline 898 SSVectorBase<Rational>& SSVectorBase<Rational>::assign(const SVectorBase<Rational>& rhs) 899 { 900 assert(rhs.dim() <= VectorBase<Rational>::dim()); 901 902 int s = rhs.size(); 903 num = 0; 904 905 for(int i = 0; i < s; ++i) 906 { 907 int k = rhs.index(i); 908 const Rational& v = rhs.value(i); 909 910 if(v == 0) 911 VectorBase<Rational>::val[k] = 0; 912 else 913 { 914 VectorBase<Rational>::val[k] = v; 915 idx[num++] = k; 916 } 917 } 918 919 setupStatus = true; 920 921 assert(isConsistent()); 922 923 return *this; 924 } 925 926 927 928 /// Assignment operator. 929 template < class R > 930 template < class S > 931 inline 932 SSVectorBase<R>& SSVectorBase<R>::operator=(const SVectorBase<S>& rhs) 933 { 934 clear(); 935 936 return assign(rhs); 937 } 938 939 940 941 // --------------------------------------------------------------------------------------------------------------------- 942 // Methods of SVectorBase 943 // --------------------------------------------------------------------------------------------------------------------- 944 945 946 947 /// Assignment operator. 948 template < class R > 949 template < class S > 950 inline 951 SVectorBase<R>& SVectorBase<R>::operator=(const VectorBase<S>& vec) 952 { 953 int n = 0; 954 Nonzero<R>* e = m_elem; 955 956 clear(); 957 958 for(int i = vec.dim() - 1; i >= 0; --i) 959 { 960 if(vec[i] != 0) 961 { 962 assert(n < max()); 963 964 e->idx = i; 965 e->val = vec[i]; 966 ++e; 967 ++n; 968 } 969 } 970 971 set_size(n); 972 973 return *this; 974 } 975 976 977 978 /// Assignment operator (specialization for Real). 979 template <> 980 template < class S > 981 inline 982 SVectorBase<Real>& SVectorBase<Real>::operator=(const VectorBase<S>& vec) 983 { 984 int n = 0; 985 Nonzero<Real>* e = m_elem; 986 987 clear(); 988 989 for(int i = vec.dim() - 1; i >= 0; --i) 990 { 991 if(vec[i] != 0) 992 { 993 assert(n < max()); 994 995 e->idx = i; 996 e->val = Real(vec[i]); 997 ++e; 998 ++n; 999 } 1000 } 1001 1002 set_size(n); 1003 1004 return *this; 1005 } 1006 1007 1008 1009 /// Assignment operator. 1010 template < class R > 1011 template < class S > 1012 inline 1013 SVectorBase<R>& SVectorBase<R>::operator=(const SSVectorBase<S>& sv) 1014 { 1015 assert(sv.isSetup()); 1016 assert(max() >= sv.size()); 1017 1018 int nnz = 0; 1019 int idx; 1020 1021 Nonzero<R>* e = m_elem; 1022 1023 for(int i = 0; i < nnz; ++i) 1024 { 1025 idx = sv.index(i); 1026 1027 if(sv.value(idx) != 0.0) 1028 { 1029 e->idx = idx; 1030 e->val = sv[idx]; 1031 ++e; 1032 ++nnz; 1033 } 1034 } 1035 1036 set_size(nnz); 1037 1038 return *this; 1039 } 1040 1041 1042 1043 /// Inner product. 1044 template < class R > 1045 inline 1046 R SVectorBase<R>::operator*(const VectorBase<R>& w) const 1047 { 1048 StableSum<R> x; 1049 Nonzero<R>* e = m_elem; 1050 1051 for(int i = size() - 1; i >= 0; --i) 1052 { 1053 x += e->val * w[e->idx]; 1054 e++; 1055 } 1056 1057 return x; 1058 } 1059 1060 1061 1062 // --------------------------------------------------------------------------------------------------------------------- 1063 // Methods of DSVectorBase 1064 // --------------------------------------------------------------------------------------------------------------------- 1065 1066 1067 1068 /// Copy constructor. 1069 template < class R > 1070 template < class S > 1071 inline 1072 DSVectorBase<R>::DSVectorBase(const VectorBase<S>& vec) 1073 : theelem(0) 1074 { 1075 allocMem((vec.dim() < 1) ? 2 : vec.dim()); 1076 *this = vec; 1077 1078 assert(isConsistent()); 1079 } 1080 1081 1082 1083 /// Copy constructor. 1084 template < class R > 1085 template < class S > 1086 inline 1087 DSVectorBase<R>::DSVectorBase(const SSVectorBase<S>& old) 1088 : theelem(0) 1089 { 1090 allocMem(old.size() < 1 ? 2 : old.size()); 1091 SVectorBase<R>::operator=(old); 1092 1093 assert(isConsistent()); 1094 } 1095 1096 1097 1098 /// Assignment operator. 1099 template < class R > 1100 template < class S > 1101 inline 1102 DSVectorBase<R>& DSVectorBase<R>::operator=(const VectorBase<S>& vec) 1103 { 1104 assert(this != (const DSVectorBase<R>*)(&vec)); 1105 1106 SVectorBase<R>::clear(); 1107 setMax(vec.dim()); 1108 SVectorBase<R>::operator=(vec); 1109 1110 assert(isConsistent()); 1111 1112 return *this; 1113 } 1114 1115 1116 1117 /// Assignment operator. 1118 template < class R > 1119 template < class S > 1120 inline 1121 DSVectorBase<R>& DSVectorBase<R>::operator=(const SSVectorBase<S>& vec) 1122 { 1123 assert(this != &vec); 1124 1125 SVectorBase<R>::clear(); 1126 makeMem(vec.size()); 1127 SVectorBase<R>::operator=(vec); 1128 1129 return *this; 1130 } 1131 1132 1133 1134 // --------------------------------------------------------------------------------------------------------------------- 1135 // Operators 1136 // --------------------------------------------------------------------------------------------------------------------- 1137 1138 1139 1140 /// Output operator. 1141 template < class R > 1142 inline 1143 std::ostream& operator<<(std::ostream& s, const VectorBase<R>& vec) 1144 { 1145 int i; 1146 1147 s << '('; 1148 1149 for(i = 0; i < vec.dim() - 1; ++i) 1150 s << vec[i] << ", "; 1151 1152 s << vec[i] << ')'; 1153 1154 return s; 1155 } 1156 1157 1158 1159 /// Subtraction. 1160 template < class R > 1161 inline 1162 VectorBase<R> operator-(const SVectorBase<R>& v, const VectorBase<R>& w) 1163 { 1164 VectorBase<R> res(w.dim()); 1165 1166 for(int i = 0; i < res.dim(); ++i) 1167 res[i] = -w[i]; 1168 1169 res += v; 1170 1171 return res; 1172 } 1173 1174 1175 1176 1177 /// Scaling. 1178 template < class R > 1179 inline 1180 DSVectorBase<R> operator*(const SVectorBase<R>& v, R x) 1181 { 1182 DSVectorBase<R> res(v.size()); 1183 1184 for(int i = 0; i < v.size(); ++i) 1185 res.add(v.index(i), v.value(i) * x); 1186 1187 return res; 1188 } 1189 1190 1191 1192 /// Scaling. 1193 template < class R > 1194 inline 1195 DSVectorBase<R> operator*(R x, const SVectorBase<R>& v) 1196 { 1197 return v * x; 1198 } 1199 1200 1201 1202 template < class R > 1203 inline 1204 std::istream& operator>>(std::istream& s, VectorBase<R>& vec) 1205 { 1206 char c; 1207 R val; 1208 int i = 0; 1209 1210 while(s.get(c).good()) 1211 { 1212 if(c != ' ' && c != '\t' && c != '\n') 1213 break; 1214 } 1215 1216 if(c != '(') 1217 s.putback(c); 1218 else 1219 { 1220 do 1221 { 1222 s >> val; 1223 1224 if(i >= vec.dim() - 1) 1225 vec.reDim(i + 16); 1226 1227 vec[i++] = val; 1228 1229 while(s.get(c).good()) 1230 { 1231 if(c != ' ' && c != '\t' && c != '\n') 1232 break; 1233 } 1234 1235 if(c != ',') 1236 { 1237 if(c != ')') 1238 s.putback(c); 1239 1240 break; 1241 } 1242 } 1243 while(s.good()); 1244 } 1245 1246 vec.reDim(i); 1247 1248 return s; 1249 } 1250 1251 1252 1253 /// Output operator. 1254 template < class R > 1255 inline 1256 std::ostream& operator<<(std::ostream& os, const SVectorBase<R>& v) 1257 { 1258 for(int i = 0, j = 0; i < v.size(); ++i) 1259 { 1260 if(j) 1261 { 1262 if(v.value(i) < 0) 1263 os << " - " << -v.value(i); 1264 else 1265 os << " + " << v.value(i); 1266 } 1267 else 1268 os << v.value(i); 1269 1270 os << " x" << v.index(i); 1271 j = 1; 1272 1273 if((i + 1) % 4 == 0) 1274 os << "\n\t"; 1275 } 1276 1277 return os; 1278 } 1279 1280 1281 1282 /// Output operator. 1283 template < class R > 1284 inline 1285 std::ostream& operator<<(std::ostream& os, const SVSetBase<R>& s) 1286 { 1287 for(int i = 0; i < s.num(); ++i) 1288 os << s[i] << "\n"; 1289 1290 return os; 1291 } 1292 } 1293 1294 /* reset the SOPLEX_DEBUG flag to its original value */ 1295 #undef SOPLEX_DEBUG 1296 #ifdef SOPLEX_DEBUG_BASEVECTORS 1297 #define SOPLEX_DEBUG 1298 #undef SOPLEX_DEBUG_BASEVECTORS 1299 #endif 1300 1301 #endif // _BASEVECTORS_H_ 1302