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