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