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 svectorbase.h 26 * @brief Sparse vectors. 27 */ 28 #ifndef _SVECTORBASE_H_ 29 #define _SVECTORBASE_H_ 30 31 #include <iostream> 32 #include <assert.h> 33 #include <math.h> 34 #include <cmath> 35 #include "soplex/stablesum.h" 36 37 namespace soplex 38 { 39 template < class R > class VectorBase; 40 template < class R > class SSVectorBase; 41 42 /// Sparse vector nonzero element. 43 /** SVectorBase keeps its nonzeros in an array of Nonzero%s providing members for saving the index and value. 44 */ 45 template < class R > 46 class Nonzero 47 { 48 public: 49 50 R val; ///< Value of nonzero element. 51 int idx; ///< Index of nonzero element. 52 53 template < class S > 54 Nonzero<R>& operator=(const Nonzero<S>& vec) 55 { 56 // todo: is the cast really necessary? Previous code worked without a cast 57 val = (R) vec.val; 58 idx = vec.idx; 59 return *this; 60 } 61 62 template < class S > 63 Nonzero(const Nonzero<S>& vec) 64 : val(vec.val) 65 , idx(vec.idx) 66 { 67 } 68 69 Nonzero() 70 : val() 71 , idx(0) 72 { 73 } 74 }; 75 76 77 78 // specialized assignment operator 79 template <> 80 template < class S > 81 Nonzero<Real>& Nonzero<Real>::operator=(const Nonzero<S>& vec) 82 { 83 val = Real(vec.val); 84 idx = vec.idx; 85 return *this; 86 } 87 88 89 /**@brief Sparse vectors. 90 * @ingroup Algebra 91 * 92 * Class SVectorBase provides packed sparse vectors. Such are a sparse vectors, with a storage scheme that keeps all 93 * data in one contiguous block of memory. This is best suited for using them for parallel computing on a distributed 94 * memory multiprocessor. 95 * 96 * SVectorBase does not provide any memory management (this will be done by class DSVectorBase). This means, that the 97 * constructor of SVectorBase expects memory where to save the nonzeros. Further, adding nonzeros to an SVectorBase may 98 * fail if no more memory is available for saving them (see also DSVectorBase). 99 * 100 * When nonzeros are added to an SVectorBase, they are appended to the set of nonzeros, i.e., they recieve numbers 101 * size(), size()+1 ... . An SVectorBase can hold atmost max() nonzeros, where max() is given in the constructor. When 102 * removing nonzeros, the remaining nonzeros are renumbered. However, only the numbers greater than the number of the 103 * first removed nonzero are affected. 104 * 105 * The following mathematical operations are provided by class SVectorBase (SVectorBase \p a, \p b, \p c; R \p x): 106 * 107 * <TABLE> 108 * <TR><TD>Operation</TD><TD>Description </TD><TD></TD> </TR> 109 * <TR><TD>\c -= </TD><TD>subtraction </TD><TD>\c a \c -= \c b </TD></TR> 110 * <TR><TD>\c += </TD><TD>addition </TD><TD>\c a \c += \c b </TD></TR> 111 * <TR><TD>\c * </TD><TD>skalar product</TD> 112 * <TD>\c x = \c a \c * \c b </TD></TR> 113 * <TR><TD>\c *= </TD><TD>scaling </TD><TD>\c a \c *= \c x </TD></TR> 114 * <TR><TD>maxAbs() </TD><TD>infinity norm </TD> 115 * <TD>\c a.maxAbs() == \f$\|a\|_{\infty}\f$ </TD></TR> 116 * <TR><TD>length() </TD><TD>eucledian norm</TD> 117 * <TD>\c a.length() == \f$\sqrt{a^2}\f$ </TD></TR> 118 * <TR><TD>length2()</TD><TD>square norm </TD> 119 * <TD>\c a.length2() == \f$a^2\f$ </TD></TR> 120 * </TABLE> 121 * 122 * Operators \c += and \c -= should be used with caution, since no efficient implementation is available. One should 123 * think of assigning the left handside vector to a dense VectorBase first and perform the addition on it. The same 124 * applies to the scalar product \c *. 125 * 126 * There are two numberings of the nonzeros of an SVectorBase. First, an SVectorBase is supposed to act like a linear 127 * algebra VectorBase. An \em index refers to this view of an SVectorBase: operator[]() is provided which returns the 128 * value at the given index of the vector, i.e., 0 for all indices which are not in the set of nonzeros. The other view 129 * of SVectorBase%s is that of a set of nonzeros. The nonzeros are numbered from 0 to size()-1. The methods index(int 130 * n) and value(int n) allow to access the index and value of the \p n 'th nonzero. \p n is referred to as the \em 131 * number of a nonzero. 132 * 133 * @todo SVectorBase should get a new implementation. There maybe a lot of memory lost due to padding the Nonzero 134 * structure. A better idea seems to be class SVectorBase { int size; int used; int* idx; R* val; }; which for 135 * several reasons could be faster or slower. If SVectorBase is changed, also DSVectorBase and SVSet have to be 136 * modified. 137 */ 138 template < class R > 139 class SVectorBase 140 { 141 template < class S > friend class SVectorBase; 142 143 private: 144 145 // ------------------------------------------------------------------------------------------------------------------ 146 /**@name Data */ 147 ///@{ 148 149 Nonzero<R>* m_elem; 150 int memsize; 151 int memused; 152 153 ///@} 154 155 public: 156 157 typedef Nonzero<R> Element; 158 159 // ------------------------------------------------------------------------------------------------------------------ 160 /**@name Access */ 161 ///@{ 162 163 /// Number of used indices. 164 int size() const 165 { 166 assert(m_elem != 0 || memused == 0); 167 return memused; 168 } 169 170 /// Maximal number of indices. 171 int max() const 172 { 173 assert(m_elem != 0 || memused == 0); 174 return memsize; 175 } 176 177 /// Dimension of the vector defined as maximal index + 1 178 int dim() const 179 { 180 const Nonzero<R>* e = m_elem; 181 int d = -1; 182 int n = size(); 183 184 while(n--) 185 { 186 d = (d > e->idx) ? d : e->idx; 187 e++; 188 } 189 190 return d + 1; 191 } 192 193 /// Position of index \p i. 194 /** @return Finds the position of the first index \p i in the index set. If no such index \p i is found, 195 * -1 is returned. Otherwise, index(pos(i)) == i holds. 196 */ 197 int pos(int i) const 198 { 199 if(m_elem != 0) 200 { 201 int n = size(); 202 203 for(int p = 0; p < n; ++p) 204 { 205 if(m_elem[p].idx == i) 206 { 207 assert(index(p) == i); 208 return p; 209 } 210 } 211 } 212 213 return -1; 214 } 215 216 /// Value to index \p i. 217 R operator[](int i) const 218 { 219 int n = pos(i); 220 221 if(n >= 0) 222 return m_elem[n].val; 223 224 return 0; 225 } 226 227 /// Reference to the \p n 'th nonzero element. 228 Nonzero<R>& element(int n) 229 { 230 assert(n >= 0); 231 assert(n < max()); 232 233 return m_elem[n]; 234 } 235 236 /// The \p n 'th nonzero element. 237 const Nonzero<R>& element(int n) const 238 { 239 assert(n >= 0); 240 assert(n < size()); 241 242 return m_elem[n]; 243 } 244 245 /// Reference to index of \p n 'th nonzero. 246 int& index(int n) 247 { 248 assert(n >= 0); 249 assert(n < size()); 250 251 return m_elem[n].idx; 252 } 253 254 /// Index of \p n 'th nonzero. 255 int index(int n) const 256 { 257 assert(n >= 0); 258 assert(n < size()); 259 260 return m_elem[n].idx; 261 } 262 263 /// Reference to value of \p n 'th nonzero. 264 R& value(int n) 265 { 266 assert(n >= 0); 267 assert(n < size()); 268 269 return m_elem[n].val; 270 } 271 272 /// Value of \p n 'th nonzero. 273 const R& value(int n) const 274 { 275 assert(n >= 0); 276 assert(n < size()); 277 278 return m_elem[n].val; 279 } 280 281 /// Append one nonzero \p (i,v). 282 void add(int i, const R& v) 283 { 284 assert(m_elem != 0); 285 assert(size() < max()); 286 287 if(v != 0.0) 288 { 289 int n = size(); 290 291 m_elem[n].idx = i; 292 m_elem[n].val = v; 293 set_size(n + 1); 294 295 assert(size() <= max()); 296 } 297 } 298 299 /// Append one uninitialized nonzero. 300 void add(int i) 301 { 302 assert(m_elem != 0); 303 assert(size() < max()); 304 305 int n = size(); 306 307 m_elem[n].idx = i; 308 set_size(n + 1); 309 310 assert(size() <= max()); 311 } 312 313 /// Append nonzeros of \p sv. 314 void add(const SVectorBase& sv) 315 { 316 add(sv.size(), sv.m_elem); 317 } 318 319 /// Append \p n nonzeros. 320 void add(int n, const int i[], const R v[]) 321 { 322 assert(n + size() <= max()); 323 324 if(n <= 0) 325 return; 326 327 int newnnz = 0; 328 329 Nonzero<R>* e = m_elem + size(); 330 331 while(n--) 332 { 333 if(*v != 0.0) 334 { 335 assert(e != nullptr); 336 e->idx = *i; 337 e->val = *v; 338 e++; 339 ++newnnz; 340 } 341 342 i++; 343 v++; 344 } 345 346 set_size(size() + newnnz); 347 } 348 349 /// Append \p n nonzeros. 350 template < class S > 351 void add(int n, const int i[], const S v[]) 352 { 353 assert(n + size() <= max()); 354 355 if(n <= 0) 356 return; 357 358 int newnnz = 0; 359 360 Nonzero<R>* e = m_elem + size(); 361 362 while(n--) 363 { 364 if(*v != R(0.0)) 365 { 366 e->idx = *i; 367 e->val = *v; 368 e++; 369 ++newnnz; 370 } 371 372 i++; 373 v++; 374 } 375 376 set_size(size() + newnnz); 377 } 378 379 /// Append \p n nonzeros. 380 void add(int n, const Nonzero<R> e[]) 381 { 382 assert(n + size() <= max()); 383 384 if(n <= 0) 385 return; 386 387 int newnnz = 0; 388 389 Nonzero<R>* ee = m_elem + size(); 390 391 while(n--) 392 { 393 if(e->val != 0.0) 394 { 395 *ee++ = *e; 396 ++newnnz; 397 } 398 399 e++; 400 } 401 402 set_size(size() + newnnz); 403 } 404 405 /// Remove nonzeros \p n thru \p m. 406 void remove(int n, int m) 407 { 408 assert(n <= m); 409 assert(m < size()); 410 assert(n >= 0); 411 412 ++m; 413 414 int cpy = m - n; 415 cpy = (size() - m >= cpy) ? cpy : size() - m; 416 417 Nonzero<R>* e = &m_elem[size() - 1]; 418 Nonzero<R>* r = &m_elem[n]; 419 420 set_size(size() - cpy); 421 422 do 423 { 424 *r++ = *e--; 425 } 426 while(--cpy); 427 } 428 429 /// Remove \p n 'th nonzero. 430 void remove(int n) 431 { 432 assert(n >= 0); 433 assert(n < size()); 434 435 int newsize = size() - 1; 436 set_size(newsize); 437 438 if(n < newsize) 439 m_elem[n] = m_elem[newsize]; 440 } 441 442 /// Remove all indices. 443 void clear() 444 { 445 set_size(0); 446 } 447 448 /// Sort nonzeros to increasing indices. 449 void sort() 450 { 451 if(m_elem != 0) 452 { 453 Nonzero<R> dummy; 454 Nonzero<R>* w; 455 Nonzero<R>* l; 456 Nonzero<R>* s = &(m_elem[0]); 457 Nonzero<R>* e = s + size(); 458 459 for(l = s, w = s + 1; w < e; l = w, ++w) 460 { 461 if(l->idx > w->idx) 462 { 463 dummy = *w; 464 465 do 466 { 467 l[1] = *l; 468 469 if(l-- == s) 470 break; 471 } 472 while(l->idx > dummy.idx); 473 474 l[1] = dummy; 475 } 476 } 477 } 478 } 479 480 ///@} 481 482 // ------------------------------------------------------------------------------------------------------------------ 483 /**@name Arithmetic operations */ 484 ///@{ 485 486 /// Maximum absolute value, i.e., infinity norm. 487 R maxAbs() const 488 { 489 R maxi = 0; 490 491 for(int i = size() - 1; i >= 0; --i) 492 { 493 if(spxAbs(m_elem[i].val) > maxi) 494 maxi = spxAbs(m_elem[i].val); 495 } 496 497 assert(maxi >= 0); 498 499 return maxi; 500 } 501 502 /// Minimum absolute value. 503 R minAbs() const 504 { 505 R mini = R(infinity); 506 507 for(int i = size() - 1; i >= 0; --i) 508 { 509 if(spxAbs(m_elem[i].val) < mini) 510 mini = spxAbs(m_elem[i].val); 511 } 512 513 assert(mini >= 0); 514 515 return mini; 516 } 517 518 /// Floating point approximation of euclidian norm (without any approximation guarantee). 519 R length() const 520 { 521 return std::sqrt(R(length2())); 522 } 523 524 /// Squared norm. 525 R length2() const 526 { 527 R x = 0; 528 int n = size(); 529 const Nonzero<R>* e = m_elem; 530 531 while(n--) 532 { 533 x += e->val * e->val; 534 e++; 535 } 536 537 return x; 538 } 539 540 /// Scaling. 541 SVectorBase<R>& operator*=(const R& x) 542 { 543 int n = size(); 544 Nonzero<R>* e = m_elem; 545 546 assert(x != 0); 547 548 while(n--) 549 { 550 e->val *= x; 551 e++; 552 } 553 554 return *this; 555 } 556 557 /// Inner product. 558 R operator*(const VectorBase<R>& w) const; 559 560 /// inner product for sparse vectors 561 template < class S > 562 R operator*(const SVectorBase<S>& w) const 563 { 564 StableSum<R> x; 565 int n = size(); 566 int m = w.size(); 567 568 if(n == 0 || m == 0) 569 return x; 570 571 int i = 0; 572 int j = 0; 573 Element* e = m_elem; 574 typename SVectorBase<S>::Element wj = w.element(j); 575 576 while(true) 577 { 578 if(e->idx == wj.idx) 579 { 580 x += e->val * wj.val; 581 i++; 582 j++; 583 584 if(i == n || j == m) 585 break; 586 587 e++; 588 wj = w.element(j); 589 } 590 else if(e->idx < wj.idx) 591 { 592 i++; 593 594 if(i == n) 595 break; 596 597 e++; 598 } 599 else 600 { 601 j++; 602 603 if(j == m) 604 break; 605 606 wj = w.element(j); 607 } 608 } 609 610 return x; 611 } 612 613 ///@} 614 615 // ------------------------------------------------------------------------------------------------------------------ 616 /**@name Constructions, destruction, and assignment */ 617 ///@{ 618 619 /// Default constructor. 620 /** The constructor expects one memory block where to store the nonzero elements. This must be passed to the 621 * constructor, where the \em number of Nonzero%s needs that fit into the memory must be given and a pointer to the 622 * beginning of the memory block. Once this memory has been passed, it shall not be modified until the SVectorBase 623 * is no longer used. 624 */ 625 explicit SVectorBase(int n = 0, Nonzero<R>* p_mem = 0) 626 { 627 setMem(n, p_mem); 628 } 629 630 SVectorBase(const SVectorBase<R>& sv) = default; 631 632 /// Assignment operator. 633 template < class S > 634 SVectorBase<R>& operator=(const VectorBase<S>& vec); 635 636 /// Assignment operator. 637 SVectorBase<R>& operator=(const SVectorBase<R>& sv) 638 { 639 if(this != &sv) 640 { 641 assert(max() >= sv.size()); 642 643 int i = sv.size(); 644 int nnz = 0; 645 Nonzero<R>* e = m_elem; 646 const Nonzero<R>* s = sv.m_elem; 647 648 while(i--) 649 { 650 assert(e != 0); 651 652 if(s->val != 0.0) 653 { 654 *e++ = *s; 655 ++nnz; 656 } 657 658 ++s; 659 } 660 661 set_size(nnz); 662 } 663 664 return *this; 665 } 666 667 /// move assignement operator. 668 SVectorBase<R>& operator=(const SVectorBase<R>&& sv) 669 { 670 if(this != &sv) 671 { 672 this->m_elem = sv.m_elem; 673 this->memsize = sv.memsize; 674 this->memused = sv.memused; 675 } 676 677 return *this; 678 } 679 680 /// Assignment operator. 681 template < class S > 682 SVectorBase<R>& operator=(const SVectorBase<S>& sv) 683 { 684 if(this != (const SVectorBase<R>*)(&sv)) 685 { 686 assert(max() >= sv.size()); 687 688 int i = sv.size(); 689 int nnz = 0; 690 Nonzero<R>* e = m_elem; 691 const Nonzero<S>* s = sv.m_elem; 692 693 while(i--) 694 { 695 assert(e != 0); 696 697 if(s->val != 0.0) 698 { 699 *e++ = *s; 700 ++nnz; 701 } 702 703 ++s; 704 } 705 706 set_size(nnz); 707 } 708 709 return *this; 710 } 711 712 /// scale and assign 713 SVectorBase<Real>& scaleAssign(int scaleExp, const SVectorBase<Real>& sv) 714 { 715 if(this != &sv) 716 { 717 assert(max() >= sv.size()); 718 719 for(int i = 0; i < sv.size(); ++i) 720 { 721 m_elem[i].val = spxLdexp(sv.value(i), scaleExp); 722 m_elem[i].idx = sv.index(i); 723 } 724 725 assert(isConsistent()); 726 } 727 728 return *this; 729 } 730 731 /// scale and assign 732 SVectorBase<Real>& scaleAssign(const int* scaleExp, const SVectorBase<Real>& sv, 733 bool negateExp = false) 734 { 735 if(this != &sv) 736 { 737 assert(max() >= sv.size()); 738 739 if(negateExp) 740 { 741 for(int i = 0; i < sv.size(); ++i) 742 { 743 m_elem[i].val = spxLdexp(sv.value(i), -scaleExp[sv.index(i)]); 744 m_elem[i].idx = sv.index(i); 745 } 746 } 747 else 748 { 749 for(int i = 0; i < sv.size(); ++i) 750 { 751 m_elem[i].val = spxLdexp(sv.value(i), scaleExp[sv.index(i)]); 752 m_elem[i].idx = sv.index(i); 753 } 754 } 755 756 assert(isConsistent()); 757 } 758 759 return *this; 760 } 761 762 763 /// Assignment operator. 764 template < class S > 765 SVectorBase<R>& assignArray(const S* rowValues, const int* rowIndices, int rowSize) 766 { 767 assert(max() >= rowSize); 768 769 int i; 770 771 for(i = 0; i < rowSize && i < max(); i++) 772 { 773 m_elem[i].val = rowValues[i]; 774 m_elem[i].idx = rowIndices[i]; 775 } 776 777 set_size(i); 778 779 return *this; 780 } 781 782 /// Assignment operator. 783 template < class S > 784 SVectorBase<R>& operator=(const SSVectorBase<S>& sv); 785 786 ///@} 787 788 // ------------------------------------------------------------------------------------------------------------------ 789 /**@name Memory */ 790 ///@{ 791 792 /// get pointer to internal memory. 793 Nonzero<R>* mem() const 794 { 795 return m_elem; 796 } 797 798 /// Set size of the vector. 799 void set_size(int s) 800 { 801 assert(m_elem != 0 || s == 0); 802 memused = s; 803 } 804 805 /// Set the maximum number of nonzeros in the vector. 806 void set_max(int m) 807 { 808 assert(m_elem != 0 || m == 0); 809 memsize = m; 810 } 811 812 /// Set the memory area where the nonzeros will be stored. 813 void setMem(int n, Nonzero<R>* elmem) 814 { 815 assert(n >= 0); 816 assert(n == 0 || elmem != 0); 817 818 m_elem = elmem; 819 set_size(0); 820 set_max(n); 821 } 822 823 ///@} 824 825 // ------------------------------------------------------------------------------------------------------------------ 826 /**@name Utilities */ 827 ///@{ 828 829 /// Consistency check. 830 bool isConsistent() const 831 { 832 #ifdef ENABLE_CONSISTENCY_CHECKS 833 834 if(m_elem != 0) 835 { 836 const int my_size = size(); 837 const int my_max = max(); 838 839 if(my_size < 0 || my_max < 0 || my_size > my_max) 840 return SPX_MSG_INCONSISTENT("SVectorBase"); 841 842 for(int i = 1; i < my_size; ++i) 843 { 844 for(int j = 0; j < i; ++j) 845 { 846 // allow trailing zeros 847 if(m_elem[i].idx == m_elem[j].idx && m_elem[i].val != 0) 848 return SPX_MSG_INCONSISTENT("SVectorBase"); 849 } 850 } 851 } 852 853 #endif 854 855 return true; 856 } 857 858 ///@} 859 }; 860 861 862 863 /// specialization for inner product for sparse vectors 864 template <> 865 template < class S > 866 Real SVectorBase<Real>::operator*(const SVectorBase<S>& w) const 867 { 868 StableSum<Real> x; 869 int n = size(); 870 int m = w.size(); 871 872 if(n == 0 || m == 0) 873 return Real(0); 874 875 int i = 0; 876 int j = 0; 877 SVectorBase<Real>::Element* e = m_elem; 878 typename SVectorBase<S>::Element wj = w.element(j); 879 880 while(true) 881 { 882 if(e->idx == wj.idx) 883 { 884 x += e->val * Real(wj.val); 885 i++; 886 j++; 887 888 if(i == n || j == m) 889 break; 890 891 e++; 892 wj = w.element(j); 893 } 894 else if(e->idx < wj.idx) 895 { 896 i++; 897 898 if(i == n) 899 break; 900 901 e++; 902 } 903 else 904 { 905 j++; 906 907 if(j == m) 908 break; 909 910 wj = w.element(j); 911 } 912 } 913 914 return x; 915 } 916 917 } // namespace soplex 918 #endif // _SVECTORBASE_H_ 919