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
17 /**@file ssvectorbase.h
18 * @brief Semi sparse vector.
19 */
20 #ifndef _SSVECTORBASE_H_
21 #define _SSVECTORBASE_H_
22
23 #include <assert.h>
24
25 #include "soplex/spxdefines.h"
26 #include "soplex/vectorbase.h"
27 #include "soplex/idxset.h"
28 #include "soplex/spxalloc.h"
29 #include "soplex/timer.h"
30 #include "soplex/stablesum.h"
31
32 namespace soplex
33 {
34 template < class R > class SVectorBase;
35 template < class R > class SVSetBase;
36
37 /**@brief Semi sparse vector.
38 * @ingroup Algebra
39 *
40 * This class implements semi-sparse vectors. Such are #VectorBase%s where the indices of its nonzero elements can be
41 * stored in an extra IdxSet. Only elements with absolute value > #epsilon are considered to be nonzero. Since really
42 * storing the nonzeros is not always convenient, an SSVectorBase provides two different stati: setup and not setup.
43 * An SSVectorBase being setup means that the nonzero indices are available, otherwise an SSVectorBase is just an
44 * ordinary VectorBase with an empty IdxSet. Note that due to arithmetic operation, zeros can slip in, i.e., it is
45 * only guaranteed that at least every non-zero is in the IdxSet.
46 */
47 template < class R >
48 class SSVectorBase : public VectorBase<R>, protected IdxSet
49 {
50 private:
51
52 friend class VectorBase<R>;
53 template < class S > friend class DSVectorBase;
54
55 // ------------------------------------------------------------------------------------------------------------------
56 /**@name Data */
57 ///@{
58
59 /// Is the SSVectorBase set up?
60 bool setupStatus;
61
62 /// A value x with |x| < epsilon is considered zero.
63 R epsilon;
64
65 /// Allocates enough space to accommodate \p newmax values.
66 void setMax(int newmax)
67 {
68 assert(idx != 0);
69 assert(newmax != 0);
70 assert(newmax >= IdxSet::size());
71
72 len = newmax;
73 spx_realloc(idx, len);
74 }
75
76 ///@}
77
78 public:
79
80 // ------------------------------------------------------------------------------------------------------------------
81 /**@name Status of an SSVectorBase
82 *
83 * An SSVectorBase can be set up or not. In case it is set up, its IdxSet correctly contains all indices of nonzero
84 * elements of the SSVectorBase. Otherwise, it does not contain any useful data. Whether or not an SSVectorBase is
85 * setup can be determined with the method \ref soplex::SSVectorBase::isSetup() "isSetup()".
86 *
87 * There are three methods for directly affecting the setup status of an SSVectorBase:
88 *
89 * - unSetup(): This method sets the status to ``not setup''.
90 *
91 * - setup(): This method initializes the IdxSet to the SSVectorBase's nonzero indices and sets the status to
92 * ``setup''.
93 *
94 * - forceSetup(): This method sets the status to ``setup'' without verifying that the IdxSet correctly contains all
95 * nonzero indices. It may be used when the nonzero indices have been computed externally.
96 */
97 ///@{
98
99 /// Only used in slufactor.hpp.
100 R* get_ptr()
101 {
102 return VectorBase<R>::get_ptr();
103 }
104 /// Returns the non-zero epsilon used.
105 R getEpsilon() const
106 {
107 return epsilon;
108 }
109
110 /// Changes the non-zero epsilon, invalidating the setup. */
111 void setEpsilon(R eps)
112 {
113 if(eps != epsilon)
114 {
115 epsilon = eps;
116 setupStatus = false;
117 }
118 }
119
120 /// Returns setup status.
121 bool isSetup() const
122 {
123 return setupStatus;
124 }
125
126 /// Makes SSVectorBase not setup.
127 void unSetup()
128 {
129 setupStatus = false;
130 }
131
132 /// Initializes nonzero indices for elements with absolute values above #epsilon and sets all other elements to 0.
133 void setup()
134 {
135 if(!isSetup())
136 {
137 IdxSet::clear();
138
139 int d = dim();
140 num = 0;
141
142 for(int i = 0; i < d; ++i)
143 {
144 if(VectorBase<R>::val[i] != R(0))
145 {
146 if(spxAbs(VectorBase<R>::val[i]) <= epsilon)
147 VectorBase<R>::val[i] = R(0);
148 else
149 {
150 idx[num] = i;
151 num++;
152 }
153 }
154 }
155
156 setupStatus = true;
157
158 assert(isConsistent());
159 }
160 }
161
162 /// Forces setup status.
163 void forceSetup()
164 {
165 setupStatus = true;
166 }
167
168 ///@}
169
170 // ------------------------------------------------------------------------------------------------------------------
171 /**@name Methods for setup SSVectorBases */
172 ///@{
173
174 /// Returns index of the \p n 'th nonzero element.
175 int index(int n) const
176 {
177 assert(isSetup());
178
179 return IdxSet::index(n);
180 }
181
182 /// Returns value of the \p n 'th nonzero element.
183 R value(int n) const
184 {
185 assert(isSetup());
186 assert(n >= 0 && n < size());
187
188 return VectorBase<R>::val[idx[n]];
189 }
190
191 /// Finds the position of index \p i in the #IdxSet, or -1 if \p i doesn't exist.
192 int pos(int i) const
193 {
194 assert(isSetup());
195
196 return IdxSet::pos(i);
197 }
198
199 /// Returns the number of nonzeros.
200 int size() const
201 {
202 assert(isSetup());
203
204 return IdxSet::size();
205 }
206
207 /// Adds nonzero (\p i, \p x) to SSVectorBase.
208 /** No nonzero with index \p i must exist in the SSVectorBase. */
209 void add(int i, R x)
210 {
211 assert(VectorBase<R>::val[i] == R(0));
212 assert(pos(i) < 0);
213
214 addIdx(i);
215 VectorBase<R>::val[i] = x;
216 }
217
218 /// Sets \p i 'th element to \p x.
219 void setValue(int i, R x)
220 {
221 assert(i >= 0);
222 assert(i < VectorBase<R>::dim());
223
224 if(isSetup())
225 {
226 int n = pos(i);
227
228 if(n < 0)
229 {
230 if(spxAbs(x) > epsilon)
231 IdxSet::add(1, &i);
232 }
233 else if(x == R(0))
234 clearNum(n);
235 }
236
237 VectorBase<R>::val[i] = x;
238
239 assert(isConsistent());
240 }
241
242 /// Scale \p i 'th element by a
243 void scaleValue(int i, int scaleExp)
244 {
245 assert(i >= 0);
246 assert(i < VectorBase<R>::dim());
247
248 VectorBase<R>::val[i] = spxLdexp(VectorBase<R>::val[i], scaleExp);
249
250 assert(isConsistent());
251 }
252
253 /// Clears element \p i.
254 void clearIdx(int i)
255 {
256 if(isSetup())
257 {
258 int n = pos(i);
259
260 if(n >= 0)
261 remove(n);
262 }
263
264 VectorBase<R>::val[i] = 0;
265
266 assert(isConsistent());
267 }
268
269 /// Sets \p n 'th nonzero element to 0 (index \p n must exist).
270 void clearNum(int n)
271 {
272 assert(isSetup());
273 assert(index(n) >= 0);
274
275 VectorBase<R>::val[index(n)] = 0;
276 remove(n);
277
278 assert(isConsistent());
279 }
280
281 ///@}
282
283 // ------------------------------------------------------------------------------------------------------------------
284 /**@name Methods independent of the Status */
285 ///@{
286
287 /// Returns \p i 'th value.
288 R operator[](int i) const
289 {
290 return VectorBase<R>::val[i];
291 }
292
293 /// Returns array indices.
294 const int* indexMem() const
295 {
296 return idx;
297 }
298
299 /// Returns array values.
300 const R* values() const
301 {
302 return VectorBase<R>::val.data();
303 }
304
305 /// Returns indices.
306 const IdxSet& indices() const
307 {
308 return *this;
309 }
310
311 /// Returns array indices.
312 int* altIndexMem()
313 {
314 unSetup();
315 return idx;
316 }
317
318 /// Returns array values.
319 R* altValues()
320 {
321 unSetup();
322 return VectorBase<R>::val.data();
323 }
324
325 /// Returns indices.
326 IdxSet& altIndices()
327 {
328 unSetup();
329 return *this;
330 }
331
332 ///@}
333
334 // ------------------------------------------------------------------------------------------------------------------
335 /**@name Arithmetic operations */
336 ///@{
337
338 /// Addition.
339 template < class S >
340 SSVectorBase<R>& operator+=(const VectorBase<S>& vec)
341 {
342 VectorBase<S>::operator+=(vec);
343
344 if(isSetup())
345 {
346 setupStatus = false;
347 setup();
348 }
349
350 return *this;
351 }
352
353 /// Addition.
354 template < class S >
355 SSVectorBase<R>& operator+=(const SVectorBase<S>& vec);
356
357 /// Addition.
358 template < class S >
359 SSVectorBase<R>& operator+=(const SSVectorBase<S>& vec)
360 {
361 assert(vec.isSetup());
362
363 for(int i = vec.size() - 1; i >= 0; --i)
364 VectorBase<R>::val[vec.index(i)] += vec.value(i);
365
366 if(isSetup())
367 {
368 setupStatus = false;
369 setup();
370 }
371
372 return *this;
373 }
374
375 /// Subtraction.
376 template < class S >
377 SSVectorBase<R>& operator-=(const VectorBase<S>& vec)
378 {
379 VectorBase<R>::operator-=(vec);
380
381 if(isSetup())
382 {
383 setupStatus = false;
384 setup();
385 }
386
387 return *this;
388 }
389
390 /// Subtraction.
391 template < class S >
392 SSVectorBase<R>& operator-=(const SVectorBase<S>& vec);
393
394 /// Subtraction.
395 template < class S >
396 SSVectorBase<R>& operator-=(const SSVectorBase<S>& vec)
397 {
398 if(vec.isSetup())
399 {
400 for(int i = vec.size() - 1; i >= 0; --i)
401 VectorBase<R>::val[vec.index(i)] -= vec.value(i);
402 }
403 else
404 VectorBase<R>::operator-=(VectorBase<S>(vec));
405
406 if(isSetup())
407 {
408 setupStatus = false;
409 setup();
410 }
411
412 return *this;
413 }
414
415 /// Scaling.
416 template < class S >
417 SSVectorBase<R>& operator*=(S x)
418 {
419 assert(isSetup());
420 assert(x != S(0));
421
422 for(int i = size() - 1; i >= 0; --i)
423 VectorBase<R>::val[index(i)] *= x;
424
425 assert(isConsistent());
426
427 return *this;
428 }
429
430 // Inner product.
431 template < class S >
432 R operator*(const SSVectorBase<S>& w)
433 {
434 setup();
435
436 StableSum<R> x;
437 int i = size() - 1;
438 int j = w.size() - 1;
439
440 // both *this and w non-zero vectors?
441 if(i >= 0 && j >= 0)
442 {
443 int vi = index(i);
444 int wj = w.index(j);
445
446 while(i != 0 && j != 0)
447 {
448 if(vi == wj)
449 {
450 x += VectorBase<R>::val[vi] * R(w.val[wj]);
451 vi = index(--i);
452 wj = w.index(--j);
453 }
454 else if(vi > wj)
455 vi = index(--i);
456 else
457 wj = w.index(--j);
458 }
459
460 /* check remaining indices */
461
462 while(i != 0 && vi != wj)
463 vi = index(--i);
464
465 while(j != 0 && vi != wj)
466 wj = w.index(--j);
467
468 if(vi == wj)
469 x += VectorBase<R>::val[vi] * R(w.val[wj]);
470 }
471
472 return x;
473 }
474
475 /// Addition of a scaled vector.
476 ///@todo SSVectorBase::multAdd() should be rewritten without pointer arithmetic.
477 template < class S, class T >
478 SSVectorBase<R>& multAdd(S xx, const SVectorBase<T>& vec);
479
480 /// Addition of a scaled vector.
481 template < class S, class T >
482 SSVectorBase<R>& multAdd(S x, const VectorBase<T>& vec)
483 {
484 VectorBase<R>::multAdd(x, vec);
485
486 if(isSetup())
487 {
488 setupStatus = false;
489 setup();
490 }
491
492 return *this;
493 }
494
495 /// Assigns pair wise vector product to SSVectorBase.
496 template < class S, class T >
497 SSVectorBase<R>& assignPWproduct4setup(const SSVectorBase<S>& x, const SSVectorBase<T>& y);
498
499 /// Assigns \f$x^T \cdot A\f$ to SSVectorBase.
500 template < class S, class T >
501 SSVectorBase<R>& assign2product(const SSVectorBase<S>& x, const SVSetBase<T>& A);
502
503 /// Assigns SSVectorBase to \f$A \cdot x\f$ for a setup \p x.
504 template < class S, class T >
505 SSVectorBase<R>& assign2product4setup(const SVSetBase<S>& A, const SSVectorBase<T>& x,
506 Timer* timeSparse, Timer* timeFull, int& nCallsSparse, int& nCallsFull);
507
508 public:
509
510 /// Assigns SSVectorBase to \f$A \cdot x\f$ thereby setting up \p x.
511 template < class S, class T >
512 SSVectorBase<R>& assign2productAndSetup(const SVSetBase<S>& A, SSVectorBase<T>& x);
513
514 /// Maximum absolute value, i.e., infinity norm.
515 R maxAbs() const
516 {
517 if(isSetup())
518 {
519 R maxabs = 0;
520
521 for(int i = 0; i < num; ++i)
522 {
523 R x = spxAbs(VectorBase<R>::val[idx[i]]);
524
525 if(x > maxabs)
526 maxabs = x;
527 }
528
529 return maxabs;
530 }
531 else
532 return VectorBase<R>::maxAbs();
533 }
534
535 /// Squared euclidian norm.
536 R length2() const
537 {
538 R x = 0;
539
540 if(isSetup())
541 {
542 for(int i = 0; i < num; ++i)
543 x += VectorBase<R>::val[idx[i]] * VectorBase<R>::val[idx[i]];
544 }
545 else
546 x = VectorBase<R>::length2();
547
548 return x;
549 }
550
551 /// Floating point approximation of euclidian norm (without any approximation guarantee).
552 R length() const
553 {
554 return spxSqrt(R(length2()));
555 }
556
557 ///@}
558
559 // ------------------------------------------------------------------------------------------------------------------
560 /**@name Miscellaneous */
561 ///@{
562
563 /// Dimension of VectorBase.
564 int dim() const
565 {
566 return VectorBase<R>::dim();
567 }
568
569 /// Resets dimension to \p newdim.
570 void reDim(int newdim)
571 {
572 for(int i = IdxSet::size() - 1; i >= 0; --i)
573 {
574 if(index(i) >= newdim)
575 remove(i);
576 }
577
578 VectorBase<R>::reDim(newdim);
579 setMax(VectorBase<R>::memSize() + 1);
580
581 assert(isConsistent());
582 }
583
584 /// Sets number of nonzeros (thereby unSetup SSVectorBase).
585 void setSize(int n)
586 {
587 assert(n >= 0);
588 assert(n <= IdxSet::max());
589
590 unSetup();
591 num = n;
592 }
593
594 /// Resets memory consumption to \p newsize.
595 void reMem(int newsize)
596 {
597 VectorBase<R>::reSize(newsize);
598 assert(isConsistent());
599
600 setMax(VectorBase<R>::memSize() + 1);
601 }
602
603 /// Clears vector.
604 void clear()
605 {
606 if(isSetup())
607 {
608 for(int i = 0; i < num; ++i)
609 VectorBase<R>::val[idx[i]] = 0;
610 }
611 else
612 VectorBase<R>::clear();
613
614 IdxSet::clear();
615 setupStatus = true;
616
617 assert(isConsistent());
618 }
619
620 /// consistency check.
621 bool isConsistent() const
622 {
623 #ifdef ENABLE_CONSISTENCY_CHECKS
624
625 if(VectorBase<R>::dim() > IdxSet::max())
626 return MSGinconsistent("SSVectorBase");
627
628 if(VectorBase<R>::dim() < IdxSet::dim())
629 return MSGinconsistent("SSVectorBase");
630
631 if(isSetup())
632 {
633 for(int i = 0; i < VectorBase<R>::dim(); ++i)
634 {
635 int j = pos(i);
636
637 if(j < 0 && spxAbs(VectorBase<R>::val[i]) > 0)
638 {
639 MSG_ERROR(std::cerr << "ESSVEC01 i = " << i
640 << "\tidx = " << j
641 << "\tval = " << std::setprecision(16) << VectorBase<R>::val[i]
642 << std::endl;)
643
644 return MSGinconsistent("SSVectorBase");
645 }
646 }
647 }
648
649 return VectorBase<R>::isConsistent() && IdxSet::isConsistent();
650 #else
651 return true;
652 #endif
653 }
654
655 ///@}
656
657 // ------------------------------------------------------------------------------------------------------------------
658 /**@name Constructors / Destructors */
659 ///@{
660
661 /// Default constructor.
662 explicit SSVectorBase<R>(int p_dim, R p_eps = Param::epsilon())
663 : VectorBase<R>(p_dim)
664 , IdxSet()
665 , setupStatus(true)
666 , epsilon(p_eps)
667 {
668 len = (p_dim < 1) ? 1 : p_dim;
669 spx_alloc(idx, len);
670 VectorBase<R>::clear();
671
672 assert(isConsistent());
673 }
674
675 /// Copy constructor.
676 template < class S >
677 SSVectorBase<R>(const SSVectorBase<S>& vec)
678 : VectorBase<R>(vec)
679 , IdxSet()
680 , setupStatus(vec.setupStatus)
681 , epsilon(vec.epsilon)
682 {
683 len = (vec.dim() < 1) ? 1 : vec.dim();
684 spx_alloc(idx, len);
685 IdxSet::operator=(vec);
686
687 assert(isConsistent());
688 }
689
690 /// Copy constructor.
691 /** The redundancy with the copy constructor below is necessary since otherwise the compiler doesn't realize that it
692 * could use the more general one with S = R and generates a shallow copy constructor.
693 */
694 SSVectorBase<R>(const SSVectorBase<R>& vec)
695 : VectorBase<R>(vec)
696 , IdxSet()
697 , setupStatus(vec.setupStatus)
698 , epsilon(vec.epsilon)
699 {
700 len = (vec.dim() < 1) ? 1 : vec.dim();
701 spx_alloc(idx, len);
702 IdxSet::operator=(vec);
703
704 assert(isConsistent());
705 }
706
707 /// Constructs nonsetup copy of \p vec.
708 template < class S >
709 explicit SSVectorBase<R>(const VectorBase<S>& vec, R eps = Param::epsilon())
710 : VectorBase<R>(vec)
711 , IdxSet()
712 , setupStatus(false)
713 , epsilon(eps)
714 {
715 len = (vec.dim() < 1) ? 1 : vec.dim();
716 spx_alloc(idx, len);
717
718 assert(isConsistent());
719 }
720
721 /// Sets up \p rhs vector, and assigns it.
722 template < class S >
723 void setup_and_assign(SSVectorBase<S>& rhs)
724 {
725 clear();
726 epsilon = rhs.epsilon;
727 setMax(rhs.max());
728 VectorBase<R>::reDim(rhs.dim());
729
730 if(rhs.isSetup())
731 {
732 IdxSet::operator=(rhs);
733
734 for(int i = size() - 1; i >= 0; --i)
735 {
736 int j = index(i);
737 VectorBase<R>::val[j] = rhs.val[j];
738 }
739 }
740 else
741 {
742 int d = rhs.dim();
743 num = 0;
744
745 for(int i = 0; i < d; ++i)
746 {
747 if(rhs.val[i] != 0)
748 {
749 if(spxAbs(rhs.val[i]) > epsilon)
750 {
751 rhs.idx[num] = i;
752 idx[num] = i;
753 VectorBase<R>::val[i] = rhs.val[i];
754 num++;
755 }
756 else
757 rhs.val[i] = 0;
758 }
759 }
760
761 rhs.num = num;
762 rhs.setupStatus = true;
763 }
764
765 setupStatus = true;
766
767 assert(rhs.isConsistent());
768 assert(isConsistent());
769 }
770
771 /// Assigns only the elements of \p rhs.
772 template < class S >
773 SSVectorBase<R>& assign(const SVectorBase<S>& rhs);
774
775 /// Assignment operator.
776 template < class S >
777 SSVectorBase<R>& operator=(const SSVectorBase<S>& rhs)
778 {
779 assert(rhs.isConsistent());
780
781 if(this != &rhs)
782 {
783 clear();
784 epsilon = rhs.epsilon;
785 setMax(rhs.max());
786 VectorBase<R>::reDim(rhs.dim());
787
788 if(rhs.isSetup())
789 {
790 IdxSet::operator=(rhs);
791
792 for(int i = size() - 1; i >= 0; --i)
793 {
794 int j = index(i);
795 VectorBase<R>::val[j] = rhs.val[j];
796 }
797 }
798 else
799 {
800 int d = rhs.dim();
801 num = 0;
802
803 for(int i = 0; i < d; ++i)
804 {
805 if(spxAbs(rhs.val[i]) > epsilon)
806 {
807 VectorBase<R>::val[i] = rhs.val[i];
808 idx[num] = i;
809 num++;
810 }
811 }
812 }
813
814 setupStatus = true;
815 }
816
817 assert(isConsistent());
818
819 return *this;
820 }
821
822 /// Assignment operator.
823 SSVectorBase<R>& operator=(const SSVectorBase<R>& rhs)
824 {
825 assert(rhs.isConsistent());
826
827 if(this != &rhs)
828 {
829 clear();
830 epsilon = rhs.epsilon;
831 setMax(rhs.max());
832 VectorBase<R>::reDim(rhs.dim());
833
834 if(rhs.isSetup())
835 {
836 IdxSet::operator=(rhs);
837
838 for(int i = size() - 1; i >= 0; --i)
839 {
840 int j = index(i);
841 VectorBase<R>::val[j] = rhs.val[j];
842 }
843 }
844 else
845 {
846 num = 0;
847
848 for(int i = 0; i < rhs.dim(); ++i)
849 {
850 if(spxAbs(rhs.val[i]) > epsilon)
851 {
852 VectorBase<R>::val[i] = rhs.val[i];
853 idx[num] = i;
854 num++;
855 }
856 }
857 }
858
859 setupStatus = true;
860 }
861
862 assert(isConsistent());
863
864 return *this;
865 }
866
867 /// Assignment operator.
868 template < class S >
869 SSVectorBase<R>& operator=(const SVectorBase<S>& rhs);
870
871 /// Assignment operator.
872 template < class S >
873 SSVectorBase<R>& operator=(const VectorBase<S>& rhs)
874 {
875 unSetup();
876 VectorBase<R>::operator=(rhs);
877
878 assert(isConsistent());
879
880 return *this;
881 }
882
883 /// destructor
884 ~SSVectorBase<R>()
885 {
886 if(idx)
887 spx_free(idx);
888 }
889
890 ///@}
891
892 private:
893
894 // ------------------------------------------------------------------------------------------------------------------
895 /**@name Private helpers */
896 ///@{
897
898 /// Assignment helper.
899 template < class S, class T >
900 SSVectorBase<R>& assign2product1(const SVSetBase<S>& A, const SSVectorBase<T>& x);
901
902 /// Assignment helper.
903 template < class S, class T >
904 SSVectorBase<R>& assign2productShort(const SVSetBase<S>& A, const SSVectorBase<T>& x);
905
906 /// Assignment helper.
907 template < class S, class T >
908 SSVectorBase<R>& assign2productFull(const SVSetBase<S>& A, const SSVectorBase<T>& x);
909
910 ///@}
911 };
912
913 } // namespace soplex
914 #endif // _SSVECTORBASE_H_
915