basevectors.h

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the class library                   */
/*       SoPlex --- the Sequential object-oriented simPlex.                  */
/*                                                                           */
/*  Copyright (c) 1996-2023 Zuse Institute Berlin (ZIB)                      */
/*                                                                           */
/*  Licensed under the Apache License, Version 2.0 (the "License");          */
/*  you may not use this file except in compliance with the License.         */
/*  You may obtain a copy of the License at                                  */
/*                                                                           */
/*      http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                           */
/*  Unless required by applicable law or agreed to in writing, software      */
/*  distributed under the License is distributed on an "AS IS" BASIS,        */
/*  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/*  See the License for the specific language governing permissions and      */
/*  limitations under the License.                                           */
/*                                                                           */
/*  You should have received a copy of the Apache-2.0 license                */
/*  along with SoPlex; see the file LICENSE. If not email to soplex@zib.de.  */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file  basevectors.h
 * @brief Collection of dense, sparse, and semi-sparse vectors.
 */
#ifndef _BASEVECTORS_H_
#define _BASEVECTORS_H_

/* undefine SOPLEX_DEBUG flag from including files; if SOPLEX_DEBUG should be defined in this file, do so below */
#ifdef SOPLEX_DEBUG
#define SOPLEX_DEBUG_BASEVECTORS
#undef SOPLEX_DEBUG
#endif

#include "soplex/spxdefines.h"
#include "soplex/rational.h"
#include "soplex/vectorbase.h"
#include "soplex/ssvectorbase.h"
#include "soplex/svectorbase.h"
#include "soplex/dsvectorbase.h"
#include "soplex/unitvectorbase.h"
#include "soplex/svsetbase.h"
#include "soplex/timer.h"

#define SOPLEX_VECTOR_MARKER   1e-100

namespace soplex
{

// ---------------------------------------------------------------------------------------------------------------------
//  Methods of VectorBase
// ---------------------------------------------------------------------------------------------------------------------

/// Assignment operator.
/** Assigning an SVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p vec.
 *  This is different in method assign().
 */



template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator=(const SVectorBase<S>& vec)
{
   clear();

   for(int i = 0; i < vec.size(); ++i)
   {
      assert(vec.index(i) < dim());
      val[vec.index(i)] = vec.value(i);
   }

   return *this;
}



/// Assign values of \p vec.
/** Assigns all nonzeros of \p vec to the vector.  All other values remain unchanged. */
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::assign(const SVectorBase<S>& vec)
{
   for(int i = vec.size() - 1; i >= 0; --i)
   {
      assert(vec.index(i) < dim());
      val[vec.index(i)] = vec.value(i);
   }

   return *this;
}



/// Assignment operator.
/** Assigning an SSVectorBase to a VectorBase using operator=() will set all values to 0 except the nonzeros of \p vec.
 *  This is different in method assign().
 */
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator=(const SSVectorBase<S>& vec)
{
   if(vec.isSetup())
   {
      clear();
      assign(vec);
   }
   else
      operator=(static_cast<const VectorBase<R>&>(vec));

   return *this;
}



/// Assign values of \p vec.
/** Assigns all nonzeros of \p vec to the vector.  All other values remain unchanged. */
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::assign(const SSVectorBase<S>& vec)
{
   assert(vec.dim() <= dim());

   if(vec.isSetup())
   {
      const int* idx = vec.indexMem();

      for(int i = vec.size() - 1; i >= 0; --i)
      {
         val[*idx] = vec.val[*idx];
         idx++;
      }
   }
   else
      operator=(static_cast<const VectorBase<R>&>(vec));

   return *this;
}



/// Addition.
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator+=(const SVectorBase<S>& vec)
{
   for(int i = vec.size() - 1; i >= 0; --i)
   {
      assert(vec.index(i) >= 0);
      assert(vec.index(i) < dim());
      val[vec.index(i)] += vec.value(i);
   }

   return *this;
}



/// Addition.
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator+=(const SSVectorBase<S>& vec)
{
   assert(dim() == vec.dim());

   if(vec.isSetup())
   {
      for(int i = vec.size() - 1; i >= 0 ; --i)
         val[vec.index(i)] += vec.value(i);
   }
   else
   {
      for(int i = dim() - 1; i >= 0; --i)
         val[i] += vec[i];
   }

   return *this;
}



/// Subtraction.
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator-=(const SVectorBase<S>& vec)
{
   for(int i = vec.size() - 1; i >= 0; --i)
   {
      assert(vec.index(i) >= 0);
      assert(vec.index(i) < dim());
      val[vec.index(i)] -= vec.value(i);
   }

   return *this;
}



/// Subtraction.
template < class R >
template < class S >
inline
VectorBase<R>& VectorBase<R>::operator-=(const SSVectorBase<S>& vec)
{
   assert(dim() == vec.dim());

   if(vec.isSetup())
   {
      for(int i = vec.size() - 1; i >= 0; --i)
         val[vec.index(i)] -= vec.value(i);
   }
   else
   {
      for(int i = dim() - 1; i >= 0; --i)
         val[i] -= vec[i];
   }

   return *this;
}



/// Inner product.
template < class R >
inline
R VectorBase<R>::operator*(const SVectorBase<R>& vec) const
{
   assert(dim() >= vec.dim());

   StableSum<R> x;

   for(int i = vec.size() - 1; i >= 0; --i)
      x += val[vec.index(i)] * vec.value(i);

   return x;
}



/// Inner product.
template < class R >
inline
R VectorBase<R>::operator*(const SSVectorBase<R>& vec) const
{
   assert(dim() == vec.dim());

   if(vec.isSetup())
   {
      const int* idx = vec.indexMem();

      StableSum<R> x;

      for(int i = vec.size() - 1; i >= 0; --i)
      {
         x += val[*idx] * vec.val[*idx];
         idx++;
      }

      return x;
   }
   else
      return operator*(static_cast<const VectorBase<R>&>(vec));
}



/// Addition of scaled vector.
template < class R >
template < class S, class T >
inline
VectorBase<R>& VectorBase<R>::multAdd(const S& x, const SVectorBase<T>& vec)
{
   for(int i = vec.size() - 1; i >= 0; --i)
   {
      assert(vec.index(i) < dim());
      val[vec.index(i)] += x * vec.value(i);
   }

   return *this;
}



/// Subtraction of scaled vector.
template < class R >
template < class S, class T >
inline
VectorBase<R>& VectorBase<R>::multSub(const S& x, const SVectorBase<T>& vec)
{
   for(int i = vec.size() - 1; i >= 0; --i)
   {
      assert(vec.index(i) < dim());
      val[vec.index(i)] -= x * vec.value(i);
   }

   return *this;
}



/// Addition of scaled vector.
template < class R >
template < class S, class T >
inline
VectorBase<R>& VectorBase<R>::multAdd(const S& x, const SSVectorBase<T>& vec)
{
   assert(vec.dim() <= dim());

   if(vec.isSetup())
   {
      const int* idx = vec.indexMem();

      for(int i = vec.size() - 1; i >= 0; --i)
         val[idx[i]] += x * vec[idx[i]];
   }
   else
   {
      assert(vec.dim() == dim());

      for(int i = dim() - 1; i >= 0; --i)
         val[i] += x * vec.val[i];
   }

   return *this;
}





// ---------------------------------------------------------------------------------------------------------------------
// Methods of SSVectorBase
// ---------------------------------------------------------------------------------------------------------------------



/// Addition.
template < class R >
template < class S >
inline
SSVectorBase<R>& SSVectorBase<R>::operator+=(const SVectorBase<S>& vec)
{
   VectorBase<R>::operator+=(vec);

   if(isSetup())
   {
      setupStatus = false;
      setup();
   }

   return *this;
}



/// Subtraction.
template < class R >
template < class S >
inline
SSVectorBase<R>& SSVectorBase<R>::operator-=(const SVectorBase<S>& vec)
{
   VectorBase<R>::operator-=(vec);

   if(isSetup())
   {
      setupStatus = false;
      setup();
   }

   return *this;
}



/// Addition of a scaled vector.
///@todo SSVectorBase::multAdd() should be rewritten without pointer arithmetic.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::multAdd(S xx, const SVectorBase<T>& vec)
{
   if(isSetup())
   {
      R* v = VectorBase<R>::val.data();
      R x;
      bool adjust = false;
      int j;

      for(int i = vec.size() - 1; i >= 0; --i)
      {
         j = vec.index(i);

         if(v[j] != 0)
         {
            x = v[j] + xx * vec.value(i);

            if(isNotZero(x, this->tolerances()->epsilon()))
               v[j] = x;
            else
            {
               adjust = true;
               v[j] = SOPLEX_VECTOR_MARKER;
            }
         }
         else
         {
            x = xx * vec.value(i);

            if(isNotZero(x, this->tolerances()->epsilon()))
            {
               v[j] = x;
               addIdx(j);
            }
         }
      }

      if(adjust)
      {
         int* iptr = idx;
         int* iiptr = idx;
         int* endptr = idx + num;

         for(; iptr < endptr; ++iptr)
         {
            x = v[*iptr];

            if(isNotZero(x, this->tolerances()->epsilon()))
               *iiptr++ = *iptr;
            else
               v[*iptr] = 0;
         }

         num = int(iiptr - idx);
      }
   }
   else
      VectorBase<R>::multAdd(xx, vec);

   assert(isConsistent());

   return *this;
}


/// Assigns pair wise vector product of setup x and setup y to SSVectorBase.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assignPWproduct4setup(const SSVectorBase<S>& x,
      const SSVectorBase<T>& y)
{
   assert(dim() == x.dim());
   assert(x.dim() == y.dim());
   assert(x.isSetup());
   assert(y.isSetup());

   clear();
   setupStatus = false;

   int i = 0;
   int j = 0;
   int n = x.size() - 1;
   int m = y.size() - 1;

   /* both x and y non-zero vectors? */
   if(m >= 0 && n >= 0)
   {
      int xi = x.index(i);
      int yj = y.index(j);

      while(i < n && j < m)
      {
         if(xi == yj)
         {
            VectorBase<R>::val[xi] = R(x.val[xi]) * R(y.val[xi]);
            xi = x.index(++i);
            yj = y.index(++j);
         }
         else if(xi < yj)
            xi = x.index(++i);
         else
            yj = y.index(++j);
      }

      /* check (possible) remaining indices */

      while(i < n && xi != yj)
         xi = x.index(++i);

      while(j < m && xi != yj)
         yj = y.index(++j);

      if(xi == yj)
         VectorBase<R>::val[xi] = R(x.val[xi]) * R(y.val[xi]);
   }

   setup();

   assert(isConsistent());

   return *this;
}



/// Assigns \f$x^T \cdot A\f$ to SSVectorBase.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2product(const SSVectorBase<S>& x, const SVSetBase<T>& A)
{
   assert(A.num() == dim());

   R y;

   clear();

   for(int i = dim() - 1; i >= 0; --i)
   {
      y = A[i] * x;

      if(isNotZero(y, this->tolerances()->epsilon()))
      {
         VectorBase<R>::val[i] = y;
         IdxSet::addIdx(i);
      }
   }

   assert(isConsistent());

   return *this;
}



/// Assigns SSVectorBase to \f$A \cdot x\f$ for a setup \p x.
#define SOPLEX_SHORTPRODUCT_FACTOR 0.5
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2product4setup(const SVSetBase<S>& A,
      const SSVectorBase<T>& x,
      Timer* timeSparse, Timer* timeFull,
      int& nCallsSparse, int& nCallsFull
                                                      )
{
   assert(A.num() == x.dim());
   assert(x.isSetup());
   clear();

   if(x.size() == 1)
   {
      if(timeSparse != 0)
         timeSparse->start();

      assign2product1(A, x);
      setupStatus = true;

      if(timeSparse != 0)
         timeSparse->stop();

      ++nCallsSparse;
   }
   else if(isSetup()
           && (double(x.size()) * A.memSize() <= SOPLEX_SHORTPRODUCT_FACTOR * dim() * A.num()))
   {
      if(timeSparse != 0)
         timeSparse->start();

      assign2productShort(A, x);
      setupStatus = true;

      if(timeSparse != 0)
         timeSparse->stop();

      ++nCallsSparse;
   }
   else
   {
      if(timeFull != 0)
         timeFull->start();

      assign2productFull(A, x);
      setupStatus = false;

      if(timeFull != 0)
         timeFull->stop();

      ++nCallsFull;
   }

   assert(isConsistent());

   return *this;
}



/// Assignment helper.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2product1(const SVSetBase<S>& A, const SSVectorBase<T>& x)
{
   assert(x.isSetup());
   assert(x.size() == 1);

   // get the nonzero value of x and the corresponding vector in A:
   const int nzidx = x.idx[0];
   const T nzval = x.val[nzidx];
   const SVectorBase<S>& Ai = A[nzidx];

   // compute A[nzidx] * nzval:
   if(isZero(nzval, this->tolerances()->epsilon()) || Ai.size() == 0)
      clear();    // this := zero vector
   else
   {
      num = Ai.size();

      for(int j = num - 1; j >= 0; --j)
      {
         const Nonzero<S>& Aij = Ai.element(j);
         idx[j] = Aij.idx;
         VectorBase<R>::val[Aij.idx] = nzval * Aij.val;
      }
   }

   assert(isConsistent());

   return *this;
}



/// Assignment helper.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2productShort(const SVSetBase<S>& A,
      const SSVectorBase<T>& x)
{
   assert(x.isSetup());

   if(x.size() == 0)   // x can be setup but have size 0 => this := zero vector
   {
      clear();
      return *this;
   }

   // compute x[0] * A[0]
   int curidx = x.idx[0];
   const T x0 = x.val[curidx];
   const SVectorBase<S>& A0 = A[curidx];
   int nonzero_idx = 0;
   int xsize = x.size();
   int Aisize;

   num = A0.size();

   if(isZero(x0, this->tolerances()->epsilon()) || num == 0)
   {
      // A[0] == 0 or x[0] == 0 => this := zero vector
      clear();
   }
   else
   {
      for(int j = 0; j < num; ++j)
      {
         const Nonzero<S>& elt = A0.element(j);
         const R product = x0 * elt.val;

         // store the value in any case
         idx[nonzero_idx] = elt.idx;
         VectorBase<R>::val[elt.idx] = product;

         // count only non-zero values; not 'isNotZero(product, epsilon)'
         if(product != 0)
            ++nonzero_idx;
      }
   }

   // Compute the other x[i] * A[i] and add them to the existing vector.
   for(int i = 1; i < xsize; ++i)
   {
      curidx = x.idx[i];
      const T xi     = x.val[curidx];
      const SVectorBase<S>& Ai = A[curidx];

      // If A[i] == 0 or x[i] == 0, do nothing.
      Aisize = Ai.size();

      if(isNotZero(xi, this->tolerances()->epsilon()) || Aisize == 0)
      {
         // Compute x[i] * A[i] and add it to the existing vector.
         for(int j = 0; j < Aisize; ++j)
         {
            const Nonzero<S>& elt = Ai.element(j);
            idx[nonzero_idx] = elt.idx;
            R oldval  = VectorBase<R>::val[elt.idx];

            // An old value of exactly 0 means the position is still unused.
            // It will be used now (either by a new nonzero or by a SOPLEX_VECTOR_MARKER),
            // so increase the counter. If oldval != 0, we just
            // change an existing NZ-element, so don't increase the counter.
            if(oldval == 0)
               ++nonzero_idx;

            // Add the current product x[i] * A[i][j]; if oldval was
            // SOPLEX_VECTOR_MARKER before, it does not hurt because SOPLEX_VECTOR_MARKER is really small.
            oldval += xi * elt.val;

            // If the new value is exactly 0, mark the index as used
            // by setting a value which is nearly 0; otherwise, store
            // the value. Values below epsilon will be removed later.
            if(oldval == 0)
               VectorBase<R>::val[elt.idx] = SOPLEX_VECTOR_MARKER;
            else
               VectorBase<R>::val[elt.idx] = oldval;
         }
      }
   }

   // Clean up by shifting all nonzeros (w.r.t. epsilon) to the front of idx,
   // zeroing all values which are nearly 0, and setting #num# appropriately.
   int nz_counter = 0;

   for(int i = 0; i < nonzero_idx; ++i)
   {
      curidx = idx[i];

      if(isZero(VectorBase<R>::val[curidx], this->tolerances()->epsilon()))
         VectorBase<R>::val[curidx] = 0;
      else
      {
         idx[nz_counter] = curidx;
         ++nz_counter;
      }

      num = nz_counter;
   }

   assert(isConsistent());

   return *this;
}



/// Assignment helper.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2productFull(const SVSetBase<S>& A,
      const SSVectorBase<T>& x)
{
   assert(x.isSetup());

   if(x.size() == 0)   // x can be setup but have size 0 => this := zero vector
   {
      clear();
      return *this;
   }

   bool A_is_zero = true;
   int xsize = x.size();
   int Aisize;

   for(int i = 0; i < xsize; ++i)
   {
      const int curidx = x.idx[i];
      const T xi = x.val[curidx];
      const SVectorBase<S>& Ai = A[curidx];
      Aisize = Ai.size();

      if(A_is_zero && Aisize > 0)
         A_is_zero = false;

      for(int j = 0; j < Aisize; ++j)
      {
         const Nonzero<S>& elt = Ai.element(j);
         VectorBase<R>::val[elt.idx] += xi * elt.val;
      }
   }

   if(A_is_zero)
      clear(); // case x != 0 but A == 0

   return *this;
}



/// Assigns SSVectorBase to \f$A \cdot x\f$ thereby setting up \p x.
template < class R >
template < class S, class T >
inline
SSVectorBase<R>& SSVectorBase<R>::assign2productAndSetup(const SVSetBase<S>& A, SSVectorBase<T>& x)
{
   assert(!x.isSetup());

   if(x.dim() == 0)
   {
      // x == 0 => this := zero vector
      clear();
      x.num = 0;
   }
   else
   {
      // x is not setup, so walk through its value vector
      int nzcount = 0;
      int end = x.dim();

      for(int i = 0; i < end; ++i)
      {
         // advance to the next element != 0
         T& xval = x.val[i];

         if(xval != 0)
         {
            // If x[i] is really nonzero, compute A[i] * x[i] and adapt x.idx,
            // otherwise set x[i] to 0.
            if(isNotZero(xval, this->tolerances()->epsilon()))
            {
               const SVectorBase<S>& Ai = A[i];
               x.idx[ nzcount++ ] = i;

               for(int j = Ai.size() - 1; j >= 0; --j)
               {
                  const Nonzero<S>& elt = Ai.element(j);
                  VectorBase<R>::val[elt.idx] += xval * elt.val;
               }
            }
            else
               xval = 0;
         }
      }

      x.num = nzcount;
      setupStatus = false;
   }

   x.setupStatus = true;

   assert(isConsistent());

   return *this;
}



/// Assigns only the elements of \p rhs.
template < class R >
template < class S >
inline
SSVectorBase<R>& SSVectorBase<R>::assign(const SVectorBase<S>& rhs)
{
   assert(rhs.dim() <= VectorBase<R>::dim());

   int s = rhs.size();
   num = 0;

   for(int i = 0; i < s; ++i)
   {
      int k = rhs.index(i);
      S v = rhs.value(i);

      if(isZero(v, this->tolerances()->epsilon()))
         VectorBase<R>::val[k] = 0;
      else
      {
         VectorBase<R>::val[k] = v;
         idx[num++] = k;
      }
   }

   setupStatus = true;

   assert(isConsistent());

   return *this;
}



/// Assigns only the elements of \p rhs.
template <  >
template <  >
inline
SSVectorBase<Rational>& SSVectorBase<Rational>::assign(const SVectorBase<Rational>& rhs)
{
   assert(rhs.dim() <= VectorBase<Rational>::dim());

   int s = rhs.size();
   num = 0;

   for(int i = 0; i < s; ++i)
   {
      int k = rhs.index(i);
      const Rational& v = rhs.value(i);

      if(v == 0)
         VectorBase<Rational>::val[k] = 0;
      else
      {
         VectorBase<Rational>::val[k] = v;
         idx[num++] = k;
      }
   }

   setupStatus = true;

   assert(isConsistent());

   return *this;
}



/// Assignment operator.
template < class R >
template < class S >
inline
SSVectorBase<R>& SSVectorBase<R>::operator=(const SVectorBase<S>& rhs)
{
   clear();

   return assign(rhs);
}



// ---------------------------------------------------------------------------------------------------------------------
//  Methods of SVectorBase
// ---------------------------------------------------------------------------------------------------------------------



/// Assignment operator.
template < class R >
template < class S >
inline
SVectorBase<R>& SVectorBase<R>::operator=(const VectorBase<S>& vec)
{
   int n = 0;
   Nonzero<R>* e = m_elem;

   clear();

   for(int i = vec.dim() - 1; i >= 0; --i)
   {
      if(vec[i] != 0)
      {
         assert(n < max());

         e->idx = i;
         e->val = vec[i];
         ++e;
         ++n;
      }
   }

   set_size(n);

   return *this;
}



/// Assignment operator (specialization for Real).
template <>
template < class S >
inline
SVectorBase<Real>& SVectorBase<Real>::operator=(const VectorBase<S>& vec)
{
   int n = 0;
   Nonzero<Real>* e = m_elem;

   clear();

   for(int i = vec.dim() - 1; i >= 0; --i)
   {
      if(vec[i] != 0)
      {
         assert(n < max());

         e->idx = i;
         e->val = Real(vec[i]);
         ++e;
         ++n;
      }
   }

   set_size(n);

   return *this;
}



/// Assignment operator.
template < class R >
template < class S >
inline
SVectorBase<R>& SVectorBase<R>::operator=(const SSVectorBase<S>& sv)
{
   assert(sv.isSetup());
   assert(max() >= sv.size());

   int nnz = 0;
   int idx;

   Nonzero<R>* e = m_elem;

   for(int i = 0; i < nnz; ++i)
   {
      idx = sv.index(i);

      if(sv.value(idx) != 0.0)
      {
         e->idx = idx;
         e->val = sv[idx];
         ++e;
         ++nnz;
      }
   }

   set_size(nnz);

   return *this;
}



/// Inner product.
template < class R >
inline
R SVectorBase<R>::operator*(const VectorBase<R>& w) const
{
   StableSum<R> x;
   Nonzero<R>* e = m_elem;

   for(int i = size() - 1; i >= 0; --i)
   {
      x += e->val * w[e->idx];
      e++;
   }

   return x;
}



// ---------------------------------------------------------------------------------------------------------------------
//  Methods of DSVectorBase
// ---------------------------------------------------------------------------------------------------------------------



/// Copy constructor.
template < class R >
template < class S >
inline
DSVectorBase<R>::DSVectorBase(const VectorBase<S>& vec)
   : theelem(0)
{
   allocMem((vec.dim() < 1) ? 2 : vec.dim());
   *this = vec;

   assert(isConsistent());
}



/// Copy constructor.
template < class R >
template < class S >
inline
DSVectorBase<R>::DSVectorBase(const SSVectorBase<S>& old)
   : theelem(0)
{
   allocMem(old.size() < 1 ? 2 : old.size());
   SVectorBase<R>::operator=(old);

   assert(isConsistent());
}



/// Assignment operator.
template < class R >
template < class S >
inline
DSVectorBase<R>& DSVectorBase<R>::operator=(const VectorBase<S>& vec)
{
   assert(this != (const DSVectorBase<R>*)(&vec));

   SVectorBase<R>::clear();
   setMax(vec.dim());
   SVectorBase<R>::operator=(vec);

   assert(isConsistent());

   return *this;
}



/// Assignment operator.
template < class R >
template < class S >
inline
DSVectorBase<R>& DSVectorBase<R>::operator=(const SSVectorBase<S>& vec)
{
   assert(this != &vec);

   SVectorBase<R>::clear();
   makeMem(vec.size());
   SVectorBase<R>::operator=(vec);

   return *this;
}



// ---------------------------------------------------------------------------------------------------------------------
//  Operators
// ---------------------------------------------------------------------------------------------------------------------



/// Output operator.
template < class R >
inline
std::ostream& operator<<(std::ostream& s, const VectorBase<R>& vec)
{
   int i;

   s << '(';

   for(i = 0; i < vec.dim() - 1; ++i)
      s << vec[i] << ", ";

   s << vec[i] << ')';

   return s;
}



/// Subtraction.
template < class R >
inline
VectorBase<R> operator-(const SVectorBase<R>& v, const VectorBase<R>& w)
{
   VectorBase<R> res(w.dim());

   for(int i = 0; i < res.dim(); ++i)
      res[i] = -w[i];

   res += v;

   return res;
}




/// Scaling.
template < class R >
inline
DSVectorBase<R> operator*(const SVectorBase<R>& v, R x)
{
   DSVectorBase<R> res(v.size());

   for(int i = 0; i < v.size(); ++i)
      res.add(v.index(i), v.value(i) * x);

   return res;
}



/// Scaling.
template < class R >
inline
DSVectorBase<R> operator*(R x, const SVectorBase<R>& v)
{
   return v * x;
}



template < class R >
inline
std::istream& operator>>(std::istream& s, VectorBase<R>& vec)
{
   char c;
   R val;
   int i = 0;

   while(s.get(c).good())
   {
      if(c != ' ' && c != '\t' && c != '\n')
         break;
   }

   if(c != '(')
      s.putback(c);
   else
   {
      do
      {
         s >> val;

         if(i >= vec.dim() - 1)
            vec.reDim(i + 16);

         vec[i++] = val;

         while(s.get(c).good())
         {
            if(c != ' ' && c != '\t' && c != '\n')
               break;
         }

         if(c != ',')
         {
            if(c != ')')
               s.putback(c);

            break;
         }
      }
      while(s.good());
   }

   vec.reDim(i);

   return s;
}



/// Output operator.
template < class R >
inline
std::ostream& operator<<(std::ostream& os, const SVectorBase<R>& v)
{
   for(int i = 0, j = 0; i < v.size(); ++i)
   {
      if(j)
      {
         if(v.value(i) < 0)
            os << " - " << -v.value(i);
         else
            os << " + " << v.value(i);
      }
      else
         os << v.value(i);

      os << " x" << v.index(i);
      j = 1;

      if((i + 1) % 4 == 0)
         os << "\n\t";
   }

   return os;
}



/// Output operator.
template < class R >
inline
std::ostream& operator<<(std::ostream& os, const SVSetBase<R>& s)
{
   for(int i = 0; i < s.num(); ++i)
      os << s[i] << "\n";

   return os;
}
}

/* reset the SOPLEX_DEBUG flag to its original value */
#undef SOPLEX_DEBUG
#ifdef SOPLEX_DEBUG_BASEVECTORS
#define SOPLEX_DEBUG
#undef SOPLEX_DEBUG_BASEVECTORS
#endif

#endif // _BASEVECTORS_H_