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