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