Scippy

SoPlex

Sequential object-oriented simPlex

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