uservector.h

Go to the documentation of this file.

// Copyright (C) 2006 Ivo Nowak and Stefan Vigerske
// All Rights Reserved.
// This code is published under the Common Public License.
//
// Author: Stefan Vigerske

#ifndef USERVECTOR_H
#define USERVECTOR_H

// TNT-stuff
#ifndef TNT_NO_BOUNDS_CHECK
#define TNT_BOUNDS_CHECK
#endif
#include <tnt_vec.h>

template <class Type> class DenseVector;

template <class Type> class UserVector {
        friend ostream& operator << (ostream& out, const UserVector<Type>& v) {
                v.print(out);
                return out;
        }
  
        friend void operator>>(const UserVector<Type>& v, Type* out) {
                for (int i=0; i<v.dim(); ++i) out[i]=v(i);
        }

        friend istream& operator >> (istream& in, UserVector<Type>& v) {
                for (int i=0; i<v.dim(); i++) in >> v[i];
                return in;
        }

  public:
    UserVector() { }

    virtual ~UserVector() { }

    virtual Pointer<UserVector<Type> > getemptycopy(int n) const=0;

    virtual Pointer<UserVector<Type> > getemptycopy() const {
      return getemptycopy(dim());
    }

    virtual Pointer<UserVector<Type> > getcopy() const {
      Pointer<UserVector<Type> > ret(getemptycopy());
      *ret=*this;
      return ret;
    }

    virtual Pointer<UserVector<Type> > getcopy(const UserVector<int>& block) const {
      Pointer<UserVector<Type> > ret(getemptycopy(block.size()));
                        for (int i=0; i<block.size(); i++) ret->SetElement(i, (*this)(block(i)));
      return ret;
    }

    virtual int dim() const=0;

    int size() const { return dim(); }

    virtual void resize(const int n) {
      assert(n>=0);
    }

    virtual Type& operator[](const int i)=0;

    virtual Type operator[](const int i) const {
      return (*this)(i);
    }

    virtual Type operator()(const int i) const=0;

    virtual DenseVector<Type> operator()(const UserVector<int>& block) const {
      return DenseVector<Type>(*this, block);
    }

    virtual DenseVector<Type> operator()(const int low, const int up) const {
      return DenseVector<Type>(*this, low, up);
    }

    virtual void SetElement(const int i, const Type& v) {
      (*this)[i]=v;
    }

    virtual void set_block(const UserVector<Type>& v, const UserVector<int>& block) {
      for (int i=0; i<block.size(); i++) (*this)[block[i]]=v(i);
    }

    virtual operator const Pointer<Type>() const=0;

    virtual UserVector<Type>& operator=(const UserVector<Type>& v)=0;

    virtual UserVector<Type>& operator=(const Type& v)=0;

    virtual void set(const Type* v, int n) {
      for (int i=0; i<n; i++) SetElement(i, v[i]);
    }

    void set(const Type* v) { set(v, dim()); }

    virtual UserVector<Type>& AddMult(Type a, const UserVector<Type>& v) {
      return (*this)+=a*v;
    }

    virtual UserVector<Type>& operator+=(const UserVector<Type>& v)=0;

    virtual UserVector<Type>& operator-=(const UserVector<Type>& v)=0;

    virtual UserVector<Type>& operator*=(const Type& v)=0;

    virtual UserVector<Type>& operator/=(const Type& v)=0;

    virtual bool operator>=(const Type& v) const {
      for (int i=0; i<dim(); i++) if ((*this)(i)<v) return false;
      return true;
    }

    virtual bool operator==(const Type& v) const {
      for (int i=0; i<dim(); i++) if ((*this)(i)!=v) return false;
      return true;
    }

    virtual bool operator==(const UserVector<Type>& v) const {
      assert(v.dim()==dim());
      for (int i=0; i<dim(); i++) if ((*this)(i)!=v(i)) return false;
      return true;
    }

                virtual Type operator*(const UserVector<Type>& v) const {
      assert(v.dim()==dim());
      Type val=Type();
      for (int i=0; i<dim(); i++) val+=(*this)(i)*v(i);
      return val;
    }

    virtual DenseVector<Type> operator*(const Type& v) const {
      return DenseVector<Type>(*this)*v;
    }

    friend DenseVector<Type> operator*(const Type& v1, const UserVector<Type>& v2) {
      return v2 * v1;
    }

    virtual DenseVector<Type> operator+(const UserVector<Type>& v) const {
      return DenseVector<Type>(*this) + v;
    }

    virtual const UserVector<Type>& operator+() const {
      return *(const UserVector<Type>*)this;
    }

    virtual UserVector<Type>& operator+() {
      return *this;
    }

    virtual DenseVector<Type> operator-(const UserVector<Type>& v) const {
      return DenseVector<Type>(*this)-v;
    }

    virtual DenseVector<Type> operator-() const {
      return -DenseVector<Type>(*this);
    }

    virtual void diagmult(UserVector<Type>& y, const UserVector<Type>& v) const {
      y=diagmult(v);
    }

    virtual DenseVector<Type> diagmult(const UserVector<Type>& v) const {
      return DenseVector<Type>(*this).diagmult(v);
    }

    virtual Type sq_norm2() const {
      Type ret=0;
      for (int i=0; i<dim(); i++) ret+=(*this)(i)*(*this)(i);
      return ret;
    }

    virtual Type dist(const UserVector<Type>& v) const {
                        return (Type)sqrt((*this-v).sq_norm2());
    }

    virtual Type mean_value() const {
      Type val=0;
      for (int i=0; i<dim(); i++) val+=(*this)(i);
                        val*=(Type)(1./dim());
      return val;
    }

    virtual Type standard_deviation() const {
      Type mval(mean_value());
      Type val=0;
      for (int i=0; i<dim(); i++) val+=(mval-(*this)(i))*(mval-(*this)(i));
                        val*=(Type)(1./dim());
      return (Type)sqrt(val);
    }

    void set_random(const Type lb, const Type ub) {
      assert(lb<=ub);
      for (int i=0; i<dim(); i++) (*this)[i]=random(lb, ub);
    }

    void set_random(const UserVector<Type>& low, const UserVector<Type>& up) {
        for (int i=0; i<dim(); i++) (*this)[i]=random(low(i), up(i));
    }

    void set_random(const Type lb, const Type ub, const double sparsity) {
      assert(0<=sparsity && sparsity<=1);
      assert(lb<=ub);
      (*this)=0;
      int ne=(dim()*sparsity);
      for (int i=0; i<ne; i++)
        (*this)[random((int)(dim()/ne)*i,(int)(dim()/ne)*(i+1))]=random(lb, ub);
    }

    virtual void print(ostream& out) const {
      for (int i=0; i<dim(); i++) out << (*this)(i) << " ";
      out << endl;
    }

};

template <class Type>
class DenseVector : public UserVector<Type> {
  public:
    static Pointer<DenseVector<Type> > random(const int n, const Type lb, const Type ub) {
      Pointer<DenseVector<Type> > ret=new DenseVector<Type>(n);
      ret->set_random(lb, ub);
      return ret;
    }

    static Pointer<DenseVector<Type> > random(const int n, const Type lb, const Type ub, const double sparsity) {
      Pointer<DenseVector<Type> > ret=new DenseVector<Type>(n);
      ret->set_random(lb, ub, sparsity);
      return ret;
    }

    static Pointer<DenseVector<Type> > random(const DenseVector<Type> low, const DenseVector<Type> up) {
        Pointer<DenseVector<Type> > ret(new DenseVector<Type>(low.dim()));
        ret->set_random(low, up);
        return ret;
    }

  protected:
    TNT::Vector<Type> x;

  public:
    DenseVector(int n=0, const Type& v=Type())
    : x(n, v)
    { }

    DenseVector(const DenseVector<Type>& v)
    : x(v.x)
    { }

    DenseVector(const UserVector<Type>& v)
    : x(v.dim())
    { for (int i=0; i<dim(); ++i) x[i]=v(i);
    }

    DenseVector(const UserVector<Type>& v, const int low, const int up)
    : x(up-low+1)
    { for (int i=low; i<=up; i++) x[i-low]=v(i);
    }

    DenseVector(const UserVector<Type>& v, const UserVector<int>& block)
    : x(block.dim())
    { for (int i=0; i<block.size(); i++) x[i]=v(block(i));
    }

    DenseVector(const vector<DenseVector<Type>* >& v, const vector<DenseVector<int> >& block)
    : x(0)
    { int s=0; for (int i=0; i<block.size(); i++) s+=block[i].size();
      resize(s);
      for (int k=0; k<block.size(); k++) if (v[k]) set_block(*v[k], block[k]);
    }

    DenseVector(const vector<Pointer<UserVector<Type> > >& v, const vector<DenseVector<int> >& block)
    : x(0)
    { int s=0; for (int i=0; i<block.size(); i++) s+=block[i].size();
      resize(s);
      for (int k=0; k<block.size(); k++) if (v[k]) set_block(*v[k], block[k]);
    }

    DenseVector(const vector<Pointer<DenseVector<Type> > >& v, const vector<DenseVector<int> >& block)
    : x(0)
    { int s=0; for (int i=0; i<block.size(); i++) s+=block[i].size();
      resize(s);
      for (int k=0; k<block.size(); k++) if (v[k]) set_block(*v[k], block[k]);
    }

    DenseVector(const vector<DenseVector<Type> >& v, const vector<DenseVector<int> >& block)
    : x(0)
    { int s=0; for (int i=0; i<block.size(); i++) s+=block[i].size();
      resize(s);
      for (int k=0; k<block.size(); k++) set_block(v[k], block[k]);
    }

    DenseVector(const Type* v, const int n)
    : x(n, v)
    { }

    DenseVector(const TNT::Vector<Type>& x_)
    : x(x_)
    { }

    Pointer<UserVector<Type> > getemptycopy(int n) const {
      return new DenseVector<Type>(n);
    }

    Pointer<UserVector<Type> > getemptycopy() const {
      return new DenseVector<Type>(dim());
    }

    Pointer<UserVector<Type> > getcopy() const {
      return new DenseVector<Type>(*this);
    }

    Pointer<UserVector<Type> > getcopy(const UserVector<int>& block) const {
      return new DenseVector<Type>(*this, block);
    }

    int dim() const {
      return x.dim();
    }

    void resize(const int n) {
      UserVector<Type>::resize(n);
      TNT::Vector<Type> tmp(x);
      x.newsize(n);
      int i=0;
      for (; i<MIN(n, tmp.size()); i++) x[i]=tmp[i];
      for (; i<n; i++) x[i]=Type();
    }

    Type& operator[](const int i) {
      return x[i];
    }

    Type operator[](const int i) const {
      return x[i];
    }

    Type operator()(const int i) const {
      return x[i];
    }

    void SetElement(const int i, const Type& v) {
      x[i]=v;
    }

    DenseVector<Type> operator()(const UserVector<int>& block) const {
      return DenseVector<Type>(*this, block);
    }

    DenseVector<Type> operator()(const int low, const int up) const {
      return DenseVector<Type>(*this, low, up);
    }

    DenseVector<Type>& operator=(const DenseVector<Type>& v) {
      x=v.x;
      return *this;
    }

    DenseVector<Type>& operator=(const UserVector<Type>& v) {
      assert(dim()==v.dim());
      for (int i=0; i<dim(); i++) x[i]=v(i);
      return *this;
    }

    DenseVector<Type>& operator=(const Type& v) {
      x=v;
      return *this;
    }

    virtual UserVector<Type>& AddMult(Type a, const UserVector<Type>& v) {
      assert(dim()==v.dim());
      for (int i=0; i<dim(); i++) x[i]+=a*v(i);
      return *this;
    }

    DenseVector<Type>& operator+=(const DenseVector<Type>& v) {
      x=x+v.x;
      return *this;
    }

    DenseVector<Type>& operator+=(const UserVector<Type>& v) {
      assert(dim()==v.dim());
      for (int i=0; i<dim(); i++) x[i]+=v(i);
      return *this;
    }

    DenseVector<Type>& operator-=(const DenseVector<Type>& v) {
      x=x-v.x;
      return *this;
    }

    DenseVector<Type>& operator-=(const UserVector<Type>& v) {
      assert(dim()==v.dim());
      for (int i=0; i<dim(); i++) x[i]-=v(i);
      return *this;
    }

    DenseVector<Type>& operator*=(const Type& v) {
      for (int i=0; i<dim(); i++) x[i]*=v;
      return *this;
    }

    DenseVector<Type>& operator/=(const Type& v) {
      for (int i=0; i<dim(); i++) x[i]/=v;
      return *this;
    }

    DenseVector<Type> operator*(const Type& v) const {
      DenseVector<Type> ret(*this);
      ret*=v;
      return ret;
    }

    friend DenseVector<Type> operator*(const Type& v1, const DenseVector<Type>& v2) {
      return v2 * v1;
    }

    Type operator*(const DenseVector<Type>& v) const {
      return TNT::dot_prod(x, v.x);
    }
                using UserVector<Type>::operator*;

    DenseVector<Type> operator+(const DenseVector<Type>& v) const {
      return DenseVector<Type>(x+v.x);
    }

    DenseVector<Type> operator+(const UserVector<Type>& v) const {
      DenseVector<Type> ret(*this);
      return (ret+=v);
    }

    const UserVector<Type>& operator+() const {
      return *(const DenseVector<Type>*)this;
    }
                using UserVector<Type>::operator+;

    DenseVector<Type> operator-(const DenseVector<Type>& v) const {
      return DenseVector<Type>(x-v.x);
    }

    DenseVector<Type> operator-(const UserVector<Type>& v) const {
      DenseVector<Type> ret(*this);
      return (ret-=v);
    }

    DenseVector<Type> operator-() const {
      DenseVector<Type> ret(*this);
      ret*=-1;
      return ret;
    }

    DenseVector<Type> diagmult(const DenseVector<Type>& v) const {
      return DenseVector<Type>(x*v.x);
    }

    DenseVector<Type> diagmult(const UserVector<Type>& v) const {
      DenseVector<Type> ret(*this);
      for (int i=0; i<dim(); i++) ret[i]*=v(i);
      return ret;
    }
                using UserVector<Type>::diagmult;

    operator const Pointer<Type>() const {
      return (const Pointer<Type>)Pointer<Type>(x.begin(), false);
    }

/*    operator Pointer<Type>() {
      return Pointer<Type>(x.begin(), false);
    }
*/
    void print(ostream& out) const {
      for (int i=0; i<dim(); i++) out << x[i] << " ";
      out << endl;
    }

};

typedef DenseVector<double> dvector;
typedef DenseVector<int> ivector;

class SparseMatrix2;
class IntervalVector;

template <class Type>
class SparseVector : public UserVector<Type> {
  friend class SparseMatrix;
        friend class IntervalVector;
  public:
    static Pointer<SparseVector<Type> > random(const int n, const Type lb, const Type ub) {
      Pointer<SparseVector<Type> > ret=new SparseVector<Type>(n);
      ret->set_random(lb, ub);
      return ret;
    }

    static Pointer<SparseVector<Type> > random(const int n, const Type lb, const Type ub, const double sparsity) {
      Pointer<SparseVector<Type> > ret=new SparseVector<Type>(n);
      ret->set_random(lb, ub, sparsity);
      return ret;
    }

    struct VectorElement {
      int index;
      Type value;
      VectorElement* next;
    };

  protected:
    VectorElement* head;
    VectorElement* last;

    int dim_;

    VectorElement* GetElement(const int i, bool create_if_not_exists=true) {
      assert(i<dim());
      if (last && (last->index<i)) {
        if (! create_if_not_exists) return NULL;
        last=last->next=new VectorElement;
        last->index=i;
        last->value=Type();
        last->next=0;
        return last;
      }
      VectorElement* p=head;
      for (VectorElement* q=p->next; q; p=q, q=q->next) {
        if (q->index==i) return q; // element exists
        if (q->index>i) { // create new one, when at right position
          if (! create_if_not_exists) return NULL;
          p->next=new VectorElement;
          p->next->index=i;
          p->next->value=Type();
          p->next->next=q;
          return p->next;
        }
      } // this is only reached, when last is 0.
      if (! create_if_not_exists) return NULL;
      last=p->next=new VectorElement;
      last->index=i;
      last->value=Type();
      last->next=0;
      return last;
    }

    void clear(VectorElement* v) {
                        VectorElement* next=v->next;
      for (VectorElement* p=next; next; p=next) { next=next->next; delete p; }
      last=v;
      last->next=NULL;
    }

  public:
                const VectorElement* gethead() { return head; }

    SparseVector(int n=0)
    : dim_(n), head(new VectorElement), last(0)
    { assert(n>=0);
      head->next=NULL;
    }

    SparseVector(int n, int i, Type v)
    : dim_(n), head(new VectorElement), last(0)
    { assert(n>=0);
      head->next=NULL;
      SetElement(i, v);
    }

    SparseVector(int n, int i, int j, Type vi, Type vj)
    : dim_(n), head(new VectorElement), last(0)
    { assert(n>=0);
      head->next=NULL;
      SetElement(i, vi);
      SetElement(j, vj);
    }

    SparseVector(const SparseVector<Type>& v)
    : dim_(v.dim()), head(new VectorElement), last(0)
    { head->next=NULL;
      for (VectorElement* p=v.head->next; p; p=p->next) {
        if (! p->value) continue;
        if (last) last=last->next=new VectorElement;
        else { last=new VectorElement; last->next=NULL; }
        if (!head->next) head->next=last;  // only in first iteration
        last->index=p->index;
        last->value=p->value;
      }
      if (last) last->next=0;
    }

    SparseVector(const UserVector<Type>& v)
    : dim_(v.dim()), head(new VectorElement), last(0)
    { head->next=NULL;
      double value;
      for (int i=0; i<v.dim(); i++) {
        if (fabs(value=v(i))<rtol) continue;
        if (last) last=last->next=new VectorElement;
        else { last=new VectorElement; last->next=NULL; }
        if (!head->next) head->next=last;
        last->index=i;
        last->value=value;
      }
      if (last) last->next=0;
    }

    SparseVector(const UserVector<Type>& v, const int low, const int up)
    : dim_(up-low+1), head(new VectorElement), last(NULL)
    { head->next=NULL;
      double value;
      for (int i=low; i<=up; i++) {
        if (fabs(value=v(i))<rtol) continue;
        if (last) last=last->next=new VectorElement;
        else { last=new VectorElement; last->next=NULL; }
        if (!head->next) head->next=last;
        last->index=i-low;
        last->value=value;
      }
      if (last) last->next=0;
    }

    SparseVector(const UserVector<Type>& v, const UserVector<int>& block)
    : dim_(block.size()), head(new VectorElement), last(NULL)
    { head->next=NULL;
        for (int i=0; i<block.size(); i++) SetElement(i, v(block(i)));
    }

    SparseVector(const vector<Pointer<UserVector<Type> > >& v, const vector<DenseVector<int> >& block)
    : dim_(0), head(new VectorElement), last(NULL)
    { head->next=NULL;
        for (int k=0; k<block.size(); k++) dim_+=block[k].size();
                        for (int k=0; k<block.size(); k++)
        if (v[k])
                                        for (int i=0; i<block[k].size(); i++) SetElement(block[k](i), (*v[k])(i));
    }

    SparseVector(const Type* v, const int& n)
    : dim_(n), head(new VectorElement), last(head)
    { head->next=NULL;
      for (int i=0; i<n; i++) {
        if (fabs(v[i])<rtol) continue;
        last=last->next=new VectorElement;
        if (!head->next) head->next=last;
        last->index=i;
        last->value=v[i];
      }
      if (last) last->next=NULL;
    }
    
    SparseVector(const int& n, const int& nz, const int* indices, const double* elements)
    : dim_(n), head(new VectorElement), last(head)
    { head->next=NULL;
      for (int i=0; i<nz; ++i) {
        if (i) assert(indices[i-1]<indices[i]);
        last=last->next=new VectorElement;
        if (!head->next) head->next=last;
        last->index=indices[i];
        last->value=elements[i];
      }
      if (last) last->next=NULL;
    }

    Pointer<UserVector<Type> > getemptycopy(int n) const {
      return new SparseVector<Type>(n);
    }

        Pointer<UserVector<Type> > getemptycopy() const {
                return new SparseVector<Type>(dim());
        }

    Pointer<UserVector<Type> > getcopy() const {
      return new SparseVector<Type>(*this);
    }

    Pointer<UserVector<Type> > getcopy(const UserVector<int>& block) const {
      return new SparseVector<Type>(*this, block);
    }

    ~SparseVector() {
       clear(head);
       delete head;
    }

    int dim() const { return dim_; }

    void resize(const int n) {
      UserVector<Type>::resize(n);
      dim_=n;
      if ((last && last->index<n) || (!head->next)) return;
      for (VectorElement* p=head; p; p=p->next)
        if (p->next && p->next->index>=n) {
          clear(p);
          last=p;
          return;
        }
    }

    void SetElement(int i, const Type& v, bool test_for_zero) {
      if (test_for_zero && fabs(v)<rtol) {
        DelElement(i);
        return;
      }
      GetElement(i)->value=v;
    }

    void SetElement(const int i, const Type& v) { SetElement(i, v, true); }

    void DelElement(int i) {
      for (VectorElement* p=head->next, *q=head; p && p->index<=i; q=p, p=p->next)
        if (p->index==i) {
          q->next=p->next;
          delete p;
          if (last==p) last=q->next;
          return;
        }
    }

    Type& operator[](const int i) {
      return GetElement(i)->value;
    }

    Type operator[](const int i) const {
      return (*this)(i);
    }

    Type operator()(const int i) const {
      assert(i<dim());
      if (last && (last->index<i)) return Type();
      for (VectorElement* p=head->next; p && p->index<=i; p=p->next)
        if (p->index==i) return p->value;
      return Type();
    }
                using UserVector<Type>::operator();

    // could be improved
    UserVector<Type>& operator=(const UserVector<Type>& v) {
      assert(v.dim()==dim());
      for (int i=0; i<v.dim(); i++) SetElement(i, v(i));
      return *this;
    }

    SparseVector<Type>& operator=(const SparseVector<Type>& v) {
      assert(v.dim()==dim());
                        clear(head);
      for (VectorElement* p=v.head->next; p; p=p->next) {
        if (! p->value) continue;
        last=last->next=new VectorElement; // last = head after clear
        last->index=p->index;
        last->value=p->value;
      }
                        if (last==head) last=NULL;
                        else last->next=NULL;
/*      VectorElement* p=head;
      VectorElement* q=v.head->next;
      for (; p->next && q; )
        if (p->next->index == q->index) {
          p->next->value=q->value;
          p=p->next; q=q->next;
        }
        else if (p->next->index < q->index) {
          VectorElement* r=p->next;
          p=p->next=r->next;
          delete r;
        }
        else {
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=q->index;
          p->value=q->value;
          p->next=r;
          q=q->next;
        }
      for (; q; q=q->next) {
        p=p->next=new VectorElement;
        p->index=q->index;
        p->value=q->value;
      }
      clear(p);
*/      return *this;
    }

    UserVector<Type>& operator=(const Type& v) {
      if (fabs(v)<rtol) {
        clear(head); last=NULL;
        return *this;
      }
      for (int i=0; i<dim(); i++) SetElement(i, v, false);
      return *this;
    }

    void set(const Type* v, int n) {
      clear(head);
      VectorElement* p=head;
      for (int i=0; i<n; i++, v++)
        if (fabs(*v)>rtol) {
          p=p->next=new VectorElement;
          p->index=i;
          p->value=*v;
        }
      p->next=NULL;
      last=p;
    }

    void set(const Type* v) { set(v, dim()); }

    UserVector<Type>& AddMult(Type a, const UserVector<Type>& v) {
      assert(dim()==v.dim());
      int i=0;
      VectorElement* p=head;
      if (p->next)
        for (; p->next->index < i && i < dim(); i++)
          if (fabs(v(i))>0) {
            VectorElement* r=p->next;
            p=p->next=new VectorElement;
            p->index=i;
            p->value=a*v(i);
            p->next=r;
          }
      for (; p->next && i<v.dim(); i++) {
        if (p->next->index == i) {
          p->next->value+=a*v(i);
          p=p->next;
        }
        else if (fabs(v(i))>0) {  // here is p->next->index > i
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=i;
          p->value=a*v(i);
          p->next=r;
        }
      }
      for (; i<v.dim(); i++)
        if (fabs(v(i))>0) {
          p=p->next=new VectorElement;
          p->index=i;
          p->value=a*v(i);
        }
      p->next=NULL;
                        if (p!=head) last=p;
      return *this;
    }

    UserVector<Type>& operator+=(const UserVector<Type>& v) {
      assert(dim()==v.dim());
      int i=0;
      VectorElement* p=head;
      if (p->next)
        for (; p->next->index < i && i < dim(); i++)
          if (fabs(v(i))>0) {
            VectorElement* r=p->next;
            p=p->next=new VectorElement;
            p->index=i;
            p->value=v(i);
            p->next=r;
          }
      for (; p->next && i<v.dim(); i++) {
        if (p->next->index == i) {
          p->next->value+=v(i);
          p=p->next;
        }
        else if (fabs(v(i))>0) {  // here is p->next->index > i
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=i;
          p->value=v(i);
          p->next=r;
        }
      }
      for (; i<v.dim(); i++)
        if (fabs(v(i))>0) {
          p=p->next=new VectorElement;
          p->index=i;
          p->value=v(i);
        }
      p->next=NULL;
      last=p;
      return *this;
    }

    SparseVector<Type>& AddMult(Type a, const SparseVector<Type>& v) {
      assert(dim()==v.dim());
      VectorElement* p=head;
      VectorElement* q=v.head->next;
      while(p->next && q)
        if (p->next->index == q->index) {
          p->next->value+=a*q->value;
          p=p->next; q=q->next;
        }
        else
          if (p->next->index < q->index) p=p->next;
        else {
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=q->index;
          p->value=a*q->value;
          p->next=r;
          q=q->next;
        }
      for(; q; q=q->next) {
        p=p->next=new VectorElement;
        p->index=q->index;
        p->value=a*q->value;
        p->next=NULL;
      }
      if (!p->next) last=p;
      return *this;
    }

    SparseVector<Type>& operator+=(const SparseVector<Type>& v) {
      assert(dim()==v.dim());
      VectorElement* p=head;
      VectorElement* q=v.head->next;
      while(p->next && q)
        if (p->next->index == q->index) {
          p->next->value+=q->value;
          p=p->next; q=q->next;
        }
        else
          if (p->next->index < q->index) p=p->next;
        else {
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=q->index;
          p->value=q->value;
          p->next=r;
          q=q->next;
        }
      for(; q; q=q->next) {
        p=p->next=new VectorElement;
        p->index=q->index;
        p->value=q->value;
        p->next=NULL;
      }
      if (!p->next) last=p;
      return *this;
    }

    UserVector<Type>& operator-=(const UserVector<Type>& v) {
      assert(dim()==v.dim());
      int i=0;
      VectorElement* p=head;
      if (p->next)
        for (; p->next->index < i && i < dim(); i++)
          if (fabs(v(i))>0) {
            VectorElement* r=p->next;
            p=p->next=new VectorElement;
            p->index=i;
            p->value=-v(i);
            p->next=r;
          }
      for (; p->next && i<v.dim(); i++) {
        if (p->next->index == i) {
          p->next->value-=v(i);
          p=p->next;
        }
        else if (fabs(v(i))>0) {  // here is p->next->index > i
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=i;
          p->value=-v(i);
          p->next=r;
        }
      }
      for (; i<v.dim(); i++)
        if (fabs(v(i))>0) {
          p=p->next=new VectorElement;
          p->index=i;
          p->value=-v(i);
        }
      p->next=NULL;
      last=p;
      return *this;
    }

    SparseVector<Type>& operator-=(const SparseVector<Type>& v) {
      assert(dim()==v.dim());
      VectorElement* p=head;
      VectorElement* q=v.head->next;
      while(p->next && q)
        if (p->next->index == q->index) {
          p->next->value-=q->value;
          p=p->next; q=q->next;
        }
        else
          if (p->next->index < q->index) p=p->next;
        else {
          VectorElement* r=p->next;
          p=p->next=new VectorElement;
          p->index=q->index;
          p->value=-q->value;
          p->next=r;
          q=q->next;
        }
      for(; q; q=q->next) {
        p=p->next=new VectorElement;
        p->index=q->index;
        p->value=-q->value;
        p->next=NULL;
      }
      if (!p->next) last=p;
      return *this;
    }

    UserVector<Type>& operator*=(const Type& v) {
      if (v==0.) {
        clear(head);
        last=NULL;
        return *this;
      }
      for (VectorElement* p=head->next; p; p=p->next) p->value*=v;
      return *this;
    }

    UserVector<Type>& operator/=(const Type& v) {
      for (VectorElement* p=head->next; p; p=p->next) p->value/=v;
      return *this;
    }

    DenseVector<Type> operator*(const Type& v) const {
      DenseVector<Type> ret(dim());
      for (VectorElement* p=head->next; p; p=p->next) ret[p->index]=v*p->value;
      return ret;
    }

    bool operator>=(const Type& v) const {
      if (v<-rtol && last && last->index<dim()-1) return false;
      for (VectorElement* p=head->next; p; p=p->next) if (p->value+rtol<v) return false;
      return true;
    }

    bool operator==(const Type& v) const {
      if (fabs(v)>rtol && last && last->index<dim()-1) return false;
      for (VectorElement* p=head->next; p; p=p->next)
        if (fabs(p->value-v)>rtol) return false;
      return true;
    }

    bool operator==(const UserVector<Type>& v) const {
      assert(v.dim()==dim());
      VectorElement* p=head->next;
      int i=0;
      if (p)
        while (p->index>i && i<dim())
          if (fabs(v(i++))>rtol) return false;
      for (; p && i<dim(); i++)
        if (p->index==i)
          if (fabs(p->value-v(i))>=rtol) return false;
          else p=p->next;
        else
          if (fabs(v(i))>=rtol) return false;
      while (i<dim()) if (fabs(v(i++))>=rtol) return false;
      return true;
    }

    Type operator*(const UserVector<Type>& v) const {
      assert(v.dim()==dim());
      Type val=0;
      for (VectorElement* p=head->next; p; p=p->next)
        val+=p->value*v(p->index);
      return val;
    }

    Type operator*(const SparseVector<Type>& v) const {
      assert(v.dim()==dim());
      Type val=0;
      for (VectorElement* p=head->next, *q=v.head; p && q; )
        if (p->index==q->index) {
          val+=p->value * q->value;
          p=p->next; q=q->next;
        }
        else if (p->index < q->index) p=p->next;
        else q=q->next;
      return val;
    }

    friend SparseVector<Type> operator*(const Type& v1, const SparseVector<Type>& v2) {
      SparseVector<Type> ret(v2);
      for (VectorElement* p=ret.head->next; p; p=p->next)
                                p->value*=v1;
                        return ret;
    }

    SparseVector<Type> operator+(const SparseVector<Type>& v) const {
      assign(v.dim()==dim());
      VectorElement* p1=head->next;
      VectorElement* p2=v.head->next;
      SparseVector<Type> ret(dim());
      ret.last=ret.head;
      while(p1 && p2) {
        ret.last=ret.last->next=new VectorElement;
        if (p1->index==p2->index) {
          ret.last->index=p1->index;
          ret.last->value=p1->value+p2->value;
          p1=p1->next; p2=p2->next;
        }
        else if (p1->index<p2->index) {
          ret.last->index=p1->index;
          ret.last->value=p1->value;
          p1=p1->next;
        }
        else {
          ret.last->index=p2->index;
          ret.last->value=p2->value;
          p2=p2->value;
        }
      }
      for (p1+=p2; p1; p1=p1->next) {  // at least one of them is NULL, so continue with the one, which is not NULL
        ret.last=ret.last->next=new VectorElement;
        ret.last->index=p1->index;
        ret.last->value=p1->value;
      }
      return ret;
    }
                using UserVector<Type>::operator+;

    SparseVector<Type> operator-(const SparseVector<Type>& v) const {
      assign(v.dim()==dim());
      VectorElement* p1=head->next;
      VectorElement* p2=v.head->next;
      SparseVector<Type> ret(dim());
      if (p1 || p2) ret.last=ret.head;
      else return ret;
      while(p1 && p2) {
        ret.last=ret.last->next=new VectorElement;
        if (p1->index==p2->index) {
          ret.last->index=p1->index;
          ret.last->value=p1->value-p2->value;
          p1=p1->next; p2=p2->next;
        }
        else if (p1->index<p2->index) {
          ret.last->index=p1->index;
          ret.last->value=p1->value;
          p1=p1->next;
        }
        else {
          ret.last->index=p2->index;
          ret.last->value=-p2->value;
          p2=p2->value;
        }
      }
      for (; p1; p1=p1->next) {
        ret.last=ret.last->next=new VectorElement;
        ret.last->index=p1->index;
        ret.last->value=p1->value;
      }
      for (; p2; p2=p2->next) {
        ret.last=ret.last->next=new VectorElement;
        ret.last->index=p2->index;
        ret.last->value=-p2->value;
      }
      return ret;
    }
                using UserVector<Type>::operator-;

    void diagmult(SparseVector<Type> y, const SparseVector<Type>& v) const {
      assert(dim()==v.dim());
      assert(dim()==y.dim());
      VectorElement* p1=head->next;
      VectorElement* p2=v.head->next;
      y=0;
      if (p1 && p2) y.last=y.head;
      else return;
      while (p1 && p2)
        if (p1->index==p2->index) {
          y.last=y.last->next=new VectorElement;
          y.last->index=p1->index;
          y.last->value=p1->value * p2->value;
          p1=p1->next; p2=p2->next;
        }
        else if (p1->index < p2->index)
          p1=p1->next;
        else p2=p2->next;
    }

    void diagmult(UserVector<Type>& y, const UserVector<Type>& v) const {
      assert(dim()==v.dim());
      assert(dim()==y.dim());
      y=0;
      for (VectorElement* p=head->next; p; p=p->next)
        y[p->index]=p->value * v(p->index);
    }
                using UserVector<Type>::diagmult;

    void set_block(const UserVector<Type>& v, const UserVector<int>& block) {
      for (int i=0; i<block.size(); i++) SetElement(block(i), v(i));
    }

    operator const Pointer<Type>() const {
      Pointer<Type> ret(new Type[dim()]);
      VectorElement* p=head->next;
      for (int i=0; i<dim(); i++)
        if (p && p->index==i) {
          ret[i]=p->value;
          p=p->next;
        } else ret[i]=Type();
      return ret;
    }

    Type sq_norm2() const {
      Type val=0;
      for (VectorElement* p=head->next; p; p=p->next) val+=p->value*p->value;
      return val;
    }

    void print(ostream& out) const {
      for (VectorElement* p=head->next; p; p=p->next) out << p->index << ": " << p->value << " ";
      out << endl;
    }

};


#ifdef FILIB_AVAILABLE
class IntervalVector : public DenseVector<interval<double> > {
        public:
                IntervalVector(int n=0, const double& v=0.)
                : DenseVector<interval<double> >(n, v)
                { }

                IntervalVector(const IntervalVector& v)
                : DenseVector<interval<double> >(v)
                { }

    IntervalVector(const UserVector<double>& v)
    : DenseVector<interval<double> >(v.dim())
    { for (int i=0; i<dim(); i++) x[i]=v(i);
    }

                IntervalVector(const UserVector<double>& v1, const UserVector<double>& v2)
    : DenseVector<interval<double> >(v1.dim())
    { for (int i=0; i<dim(); i++) x[i]=interval<double>(v1(i), v2(i));
    }

    IntervalVector(const IntervalVector& v, const int low, const int up)
    : DenseVector<interval<double> >(v, low, up)
    { }

    IntervalVector(const IntervalVector& v, const UserVector<int>& block)
    : DenseVector<interval<double> >(v, block)
    { }

    IntervalVector(const vector<IntervalVector>& v, const vector<DenseVector<int> >& block)
    {   int s=0; for (int i=0; i<block.size(); i++) s+=block[i].size();
      resize(s);
      for (int k=0; k<block.size(); k++) set_block(v[k], block[k]);
                }

                IntervalVector& operator=(const IntervalVector& v) {
                        DenseVector<interval<double> >::operator=(v);
                        return *this;
                }
                using DenseVector<interval<double> >::operator=;

                interval<double> operator*(const UserVector<double>& v) const {
      assert(v.dim()==dim());
      interval<double> val;
      for (int i=0; i<dim(); i++) val+=(*this)(i)*v(i);
      return val;
    }

                interval<double> operator*(const SparseVector<double>& v) const {
      assert(v.dim()==dim());
      interval<double> val;
      for (SparseVector<double>::VectorElement* p=v.head->next; p; p=p->next)
        val+=(*this)(p->index)*p->value;
      return val;
    }
                using DenseVector<interval<double> >::operator*;

                void AddMult(const double& a, const UserVector<double>& v) {
                        for (int i=0; i<v.dim(); i++) (*this)[i]+=a*v(i);
                }

                void AddMult(const double& a, const SparseVector<double>& v) {
                        for (SparseVector<double>::VectorElement* p=v.head->next; p; p=p->next)
                (*this)[p->index]+=a*p->value;
                }
                using DenseVector<interval<double> >::AddMult;

    IntervalVector diagmult(const UserVector<double>& v) const {
      IntervalVector ret(*this);
      for (int i=0; i<dim(); i++) ret[i]*=v(i);
      return ret;
    }
                using DenseVector<interval<double> >::diagmult;
};
#endif

class Project {
        public:
                static dvector project(const dvector& x, const dvector& lower, const dvector& upper) {
                        dvector xp(x.dim());
                        project(xp, x, lower, upper);
                        return xp;
                }
                        
                static void project(dvector& xp, const dvector& x, const dvector& lower, const dvector& upper) {
                for (int j=0; j<x.dim(); j++)
                          if (x(j) < lower(j)) xp[j]=lower(j);
                          else if (x(j) > upper(j)) xp[j]=upper(j);
                          else xp[j]=x(j);
                }
                
};

class Round {
        public:
                static dvector round(const dvector& x, vector<int>& indices, const dvector& lower, const dvector& upper) {
                        dvector xr(x.dim());
                        round(xr, x, indices, lower, upper);
                        return xr;
                }
                
                static void round(dvector& xr, const dvector& x, vector<int>& indices, const dvector& lower, const dvector& upper) {
                        xr=x;
                for (vector<int>::iterator i(indices.begin()); i!=indices.end(); i++)
                          xr[*i] = project(closestint(x(*i)), lower(*i), upper(*i)); 
                }
};


#endif // USERVECTOR_H

---