/* -*- c++ -*- (enables emacs c++ mode) */
/*===========================================================================

 Copyright (C) 2002-2017 Yves Renard

 This file is a part of GetFEM++

 GetFEM++  is  free software;  you  can  redistribute  it  and/or modify it
 under  the  terms  of the  GNU  Lesser General Public License as published
 by  the  Free Software Foundation;  either version 3 of the License,  or
 (at your option) any later version along with the GCC Runtime Library
 Exception either version 3.1 or (at your option) any later version.
 This program  is  distributed  in  the  hope  that it will be useful,  but
 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 or  FITNESS  FOR  A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 License and GCC Runtime Library Exception for more details.
 You  should  have received a copy of the GNU Lesser General Public License
 along  with  this program;  if not, write to the Free Software Foundation,
 Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301, USA.

 As a special exception, you  may use  this file  as it is a part of a free
 software  library  without  restriction.  Specifically,  if   other  files
 instantiate  templates  or  use macros or inline functions from this file,
 or  you compile this  file  and  link  it  with other files  to produce an
 executable, this file  does  not  by itself cause the resulting executable
 to be covered  by the GNU Lesser General Public License.  This   exception
 does not  however  invalidate  any  other  reasons why the executable file
 might be covered by the GNU Lesser General Public License.

===========================================================================*/
/**@file gmm_vector.h
   @author  Yves Renard <Yves.Renard@insa-lyon.fr>
   @date October 13, 2002.
   @brief Declaration of the vector types (gmm::rsvector, gmm::wsvector,
     gmm::slvector ,..)
*/
#ifndef GMM_VECTOR_H__
#define GMM_VECTOR_H__

#include <map>
#include "gmm_interface.h"

namespace gmm {

  /*************************************************************************/
  /*                                                                       */
  /* Class ref_elt_vector: reference on a vector component.                */
  /*                                                                       */
  /*************************************************************************/


  template<typename T, typename V> class ref_elt_vector {

    V *pm;
    size_type l;
    
    public :

    operator T() const { return pm->r(l); }
    ref_elt_vector(V *p, size_type ll) : pm(p), l(ll) {}
    inline bool operator ==(T v) const { return ((*pm).r(l) == v); }
    inline bool operator !=(T v) const { return ((*pm).r(l) != v); }
    inline bool operator ==(std::complex<T> v) const
    { return ((*pm).r(l) == v); }
    inline bool operator !=(std::complex<T> v) const
    { return ((*pm).r(l) != v); }
    inline ref_elt_vector &operator +=(T v)
    { (*pm).wa(l, v); return *this; }
    inline ref_elt_vector &operator -=(T v)
    { (*pm).wa(l, -v); return *this; }
    inline ref_elt_vector &operator /=(T v)
    { (*pm).w(l,(*pm).r(l) / v); return *this; }
    inline ref_elt_vector &operator *=(T v)
    { (*pm).w(l,(*pm).r(l) * v); return *this; }
    inline ref_elt_vector &operator =(const ref_elt_vector &re)
    { *this = T(re); return *this; }
    inline ref_elt_vector &operator =(T v)
    { (*pm).w(l,v); return *this; }
    T operator +()    { return  T(*this);   }
    T operator -()    { return -T(*this);   }
    T operator +(T v) { return T(*this)+ v; }
    T operator -(T v) { return T(*this)- v; }
    T operator *(T v) { return T(*this)* v; }
    T operator /(T v) { return T(*this)/ v; }
    std::complex<T> operator +(std::complex<T> v) { return T(*this)+ v; }
    std::complex<T> operator -(std::complex<T> v) { return T(*this)- v; }
    std::complex<T> operator *(std::complex<T> v) { return T(*this)* v; }
    std::complex<T> operator /(std::complex<T> v) { return T(*this)/ v; }
  };

  template<typename T, typename V> class ref_elt_vector<std::complex<T>,V> {

    V *pm;
    size_type l;
    
    public :

    operator std::complex<T>() const { return pm->r(l); }
    ref_elt_vector(V *p, size_type ll) : pm(p), l(ll) {}
    inline bool operator ==(std::complex<T> v) const
    { return ((*pm).r(l) == v); }
    inline bool operator !=(std::complex<T> v) const
    { return ((*pm).r(l) != v); }
    inline bool operator ==(T v) const { return ((*pm).r(l) == v); }
    inline bool operator !=(T v) const { return ((*pm).r(l) != v); }
    inline ref_elt_vector &operator +=(std::complex<T> v)
    { (*pm).w(l,(*pm).r(l) + v); return *this; }
    inline ref_elt_vector &operator -=(std::complex<T> v)
    { (*pm).w(l,(*pm).r(l) - v); return *this; }
    inline ref_elt_vector &operator /=(std::complex<T> v)
    { (*pm).w(l,(*pm).r(l) / v); return *this; }
    inline ref_elt_vector &operator *=(std::complex<T> v)
    { (*pm).w(l,(*pm).r(l) * v); return *this; }
    inline ref_elt_vector &operator =(const ref_elt_vector &re)
    { *this = T(re); return *this; }
    inline ref_elt_vector &operator =(std::complex<T> v)
    { (*pm).w(l,v); return *this; }
    inline ref_elt_vector &operator =(T v)
    { (*pm).w(l,std::complex<T>(v)); return *this; }
    inline ref_elt_vector &operator +=(T v)
    { (*pm).w(l,(*pm).r(l) + v); return *this; }
    inline ref_elt_vector &operator -=(T v)
    { (*pm).w(l,(*pm).r(l) - v); return *this; }
    inline ref_elt_vector &operator /=(T v)
    { (*pm).w(l,(*pm).r(l) / v); return *this; }
    inline ref_elt_vector &operator *=(T v)
    { (*pm).w(l,(*pm).r(l) * v); return *this; }
    std::complex<T> operator +()    { return  std::complex<T>(*this);   }
    std::complex<T> operator -()    { return -std::complex<T>(*this);   }
    std::complex<T> operator +(T v) { return std::complex<T>(*this)+ v; }
    std::complex<T> operator -(T v) { return std::complex<T>(*this)- v; }
    std::complex<T> operator *(T v) { return std::complex<T>(*this)* v; }
    std::complex<T> operator /(T v) { return std::complex<T>(*this)/ v; }
    std::complex<T> operator +(std::complex<T> v)
    { return std::complex<T>(*this)+ v; }
    std::complex<T> operator -(std::complex<T> v)
    { return std::complex<T>(*this)- v; }
    std::complex<T> operator *(std::complex<T> v)
    { return std::complex<T>(*this)* v; }
    std::complex<T> operator /(std::complex<T> v)
    { return std::complex<T>(*this)/ v; }
  };  

  
  template<typename T, typename V> inline
  bool operator ==(T v, const ref_elt_vector<T, V> &re) { return (v==T(re)); }
  template<typename T, typename V> inline
  bool operator !=(T v, const ref_elt_vector<T, V> &re) { return (v!=T(re)); }
  template<typename T, typename V> inline
  T &operator +=(T &v, const ref_elt_vector<T, V> &re) 
  { v += T(re); return v; }
  template<typename T, typename V> inline
  T &operator -=(T &v, const ref_elt_vector<T, V> &re)
  { v -= T(re); return v; }
  template<typename T, typename V> inline
  T &operator *=(T &v, const ref_elt_vector<T, V> &re) 
  { v *= T(re); return v; }
  template<typename T, typename V> inline
  T &operator /=(T &v, const ref_elt_vector<T, V> &re)
  { v /= T(re); return v; }
  template<typename T, typename V> inline
  T operator +(T v, const ref_elt_vector<T, V> &re) { return v+ T(re); }
  template<typename T, typename V> inline
  T operator -(T v, const ref_elt_vector<T, V> &re) { return v- T(re); }
  template<typename T, typename V> inline
  T operator *(T v, const ref_elt_vector<T, V> &re) { return v* T(re); }
  template<typename T, typename V> inline
  T operator /(T v, const ref_elt_vector<T, V> &re) { return v/ T(re); }
  template<typename T, typename V> inline
  std::complex<T> operator +(std::complex<T> v, const ref_elt_vector<T, V> &re)
  { return v+ T(re); }
  template<typename T, typename V> inline
  std::complex<T> operator -(std::complex<T> v, const ref_elt_vector<T, V> &re)
  { return v- T(re); }
  template<typename T, typename V> inline
  std::complex<T> operator *(std::complex<T> v, const ref_elt_vector<T, V> &re)
  { return v* T(re); }
  template<typename T, typename V> inline
  std::complex<T> operator /(std::complex<T> v, const ref_elt_vector<T, V> &re)
  { return v/ T(re); }
  template<typename T, typename V> inline
  std::complex<T> operator +(T v, const ref_elt_vector<std::complex<T>, V> &re)
  { return v+ std::complex<T>(re); }
  template<typename T, typename V> inline
  std::complex<T> operator -(T v, const ref_elt_vector<std::complex<T>, V> &re)
  { return v- std::complex<T>(re); }
  template<typename T, typename V> inline
  std::complex<T> operator *(T v, const ref_elt_vector<std::complex<T>, V> &re)
  { return v* std::complex<T>(re); }
  template<typename T, typename V> inline
  std::complex<T> operator /(T v, const ref_elt_vector<std::complex<T>, V> &re)
  { return v/ std::complex<T>(re); }
  template<typename T, typename V> inline
  typename number_traits<T>::magnitude_type
  abs(const ref_elt_vector<T, V> &re) { return gmm::abs(T(re)); }
  template<typename T, typename V> inline
  T sqr(const ref_elt_vector<T, V> &re) { return gmm::sqr(T(re)); }
  template<typename T, typename V> inline
  typename number_traits<T>::magnitude_type
  abs_sqr(const ref_elt_vector<T, V> &re) { return gmm::abs_sqr(T(re)); }
  template<typename T, typename V> inline
  T conj(const ref_elt_vector<T, V> &re) { return gmm::conj(T(re)); }
  template<typename T, typename V> std::ostream &operator <<
  (std::ostream &o, const ref_elt_vector<T, V> &re) { o << T(re); return o; }
  template<typename T, typename V> inline
  typename number_traits<T>::magnitude_type
  real(const ref_elt_vector<T, V> &re) { return gmm::real(T(re)); }
  template<typename T, typename V> inline
  typename number_traits<T>::magnitude_type
  imag(const ref_elt_vector<T, V> &re) { return gmm::imag(T(re)); }

  /*************************************************************************/
  /*                                                                       */
  /* Class dsvector: sparse vector optimized for random write operations   */
  /* with constant complexity for read and write operations.               */
  /* Based on distribution sort principle.                                 */
  /* Cheap for densely populated vectors.                                  */
  /*                                                                       */
  /*************************************************************************/

  template<typename T> class dsvector;

  template<typename T> struct dsvector_iterator {
    size_type i;    // Current index.
    T* p;           // Pointer to the current position.
    dsvector<T> *v; // Pointer to the vector.
    
    typedef T                   value_type;
    typedef value_type*         pointer;
    typedef const value_type*   const_pointer;
    typedef value_type&         reference;
    // typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    typedef dsvector_iterator<T> iterator;
    
    reference operator *() const { return *p; }
    pointer operator->() const { return &(operator*()); }

    iterator &operator ++() {
      for (size_type k = (i & 15); k < 15; ++k)
	{ ++p; ++i; if (*p != T(0)) return *this; }
      v->next_pos(*(const_cast<const_pointer *>(&(p))), i);
      return *this;
    }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() {
      for (size_type k = (i & 15); k > 0; --k)
	{ --p; --i; if (*p != T(0)) return *this; }
      v->previous_pos(p, i);
      return *this;
    }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &it) const
    { return (i == it.i && p == it.p && v == it.v); }
    bool operator !=(const iterator &it) const
    { return !(it == *this); }
    
    size_type index(void) const { return i; }

    dsvector_iterator(void) : i(size_type(-1)), p(0), v(0) {}
    dsvector_iterator(dsvector<T> &w) : i(size_type(-1)), p(0), v(&w) {};
  };


  template<typename T> struct dsvector_const_iterator {
    size_type i;          // Current index.
    const T* p;           // Pointer to the current position.
    const dsvector<T> *v; // Pointer to the vector.
    
    typedef T                   value_type;
    typedef const value_type*   pointer;
    typedef const value_type&   reference;
    // typedef size_t           size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    typedef dsvector_const_iterator<T> iterator;
   
    reference operator *() const { return *p; }
    pointer operator->() const { return &(operator*()); }
    iterator &operator ++() {
      for (size_type k = (i & 15); k < 15; ++k)
	{ ++p; ++i; if (*p != T(0)) return *this; }
      v->next_pos(p, i);
      return *this;
    }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() {
      for (size_type k = (i & 15); k > 0; --k)
	{ --p; --i; if (*p != T(0)) return *this; }
      v->previous_pos(p, i);
      return *this;
    }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &it) const
    { return (i == it.i && p == it.p && v == it.v); }
    bool operator !=(const iterator &it) const
    { return !(it == *this); }
    
    size_type index(void) const { return i; }

    dsvector_const_iterator(void) : i(size_type(-1)), p(0) {}
    dsvector_const_iterator(const dsvector_iterator<T> &it)
      : i(it.i), p(it.p), v(it.v) {}
    dsvector_const_iterator(const dsvector<T> &w)
      : i(size_type(-1)), p(0), v(&w) {};
  };

  
  /**
     Sparse vector built on distribution sort principle.
     Read and write access have a constant complexity depending only on the
     vector size.
  */
  template<typename T> class dsvector {

    typedef dsvector_iterator<T>       iterator;
    typedef dsvector_const_iterator<T> const_iterator;
    typedef dsvector<T>                this_type;
    typedef T *                        pointer;
    typedef const T *                  const_pointer;
    typedef void *                     void_pointer;
    typedef const void *               const_void_pointer;
 
  protected:
    size_type    n;         // Potential vector size
    size_type    depth;     // Number of row of pointer arrays
    size_type    mask;      // Mask for the first pointer array
    size_type    shift;     // Shift for the first pointer array
    void_pointer root_ptr;  // Root pointer

    const T *read_access(size_type i) const {
      GMM_ASSERT1(i < n, "index out of range");
      size_type my_mask = mask, my_shift = shift;
      void_pointer p = root_ptr;
      if (!p) return 0;
      for (size_type k = 0; k < depth; ++k) {
	p = ((void **)(p))[(i & my_mask) >> my_shift];
	if (!p) return 0;
	my_mask = (my_mask >> 4);
	my_shift -= 4;
      }
      GMM_ASSERT1(my_shift == 0, "internal error");
      GMM_ASSERT1(my_mask == 15, "internal error");
      return &(((const T *)(p))[i & 15]);
    }

    T *write_access(size_type i) {
      GMM_ASSERT1(i < n, "index " << i << " out of range (size " << n << ")");
      size_type my_mask = mask, my_shift = shift;
      if (!root_ptr) {
	if (depth) {
	  root_ptr = new void_pointer[16];
	  std::memset(root_ptr, 0, 16*sizeof(void_pointer));
	} else {
	  root_ptr = new T[16];
	  for (size_type l = 0; l < 16; ++l) ((T *)(root_ptr))[l] = T(0);
	}
      }

      void_pointer p = root_ptr;
      for (size_type k = 0; k < depth; ++k) {
	size_type j = (i & my_mask) >> my_shift;
	void_pointer q = ((void_pointer *)(p))[j];
	if (!q) {
	  if (k+1 != depth) {
	    q = new void_pointer[16];
	    std::memset(q, 0, 16*sizeof(void_pointer));
	  } else {
	    q = new T[16];
	    for (size_type l = 0; l < 16; ++l) ((T *)(q))[l] = T(0);
	  }
	  ((void_pointer *)(p))[j] = q;
	}
	p = q;
	my_mask = (my_mask >> 4);
	my_shift -= 4;
      }
      GMM_ASSERT1(my_shift == 0, "internal error");
      GMM_ASSERT1(my_mask == 15, "internal error " << my_mask);
      return &(((T *)(p))[i & 15]);
    }

    void init(size_type n_) {
      n = n_; depth = 0; shift = 0; mask = 1; if (n_) --n_;
      while (n_) { n_ /= 16; ++depth; shift += 4; mask *= 16; }
      mask--; if (shift) shift -= 4; if (depth) --depth;
      root_ptr = 0;
    }

    void rec_del(void_pointer p, size_type my_depth) {
      if (my_depth) {
	for (size_type k = 0; k < 16; ++k)
	  if (((void_pointer *)(p))[k])
	    rec_del(((void_pointer *)(p))[k], my_depth-1);
	delete[] ((void_pointer *)(p));
      } else {
	delete[] ((T *)(p));
      }
    }

    void rec_clean(void_pointer p, size_type my_depth, double eps) {
      if (my_depth) {
	for (size_type k = 0; k < 16; ++k)
	  if (((void_pointer *)(p))[k])
	    rec_clean(((void_pointer *)(p))[k], my_depth-1, eps);
      } else {
	for (size_type k = 0; k < 16; ++k)
	  if (gmm::abs(((T *)(p))[k]) <= eps) ((T *)(p))[k] = T(0);
      }
    }

    void rec_clean_i(void_pointer p, size_type my_depth, size_type my_mask,
		     size_type i, size_type base) {
      if (my_depth) {
	my_mask = (my_mask >> 4);
	for (size_type k = 0; k < 16; ++k)
	  if (((void_pointer *)(p))[k] && (base + (k+1)*(mask+1)) >= i)
	    rec_clean_i(((void_pointer *)(p))[k], my_depth-1, my_mask,
			i, base + k*(my_mask+1));
      } else {
	for (size_type k = 0; k < 16; ++k)
	  if (base+k > i) ((T *)(p))[k] = T(0);
      }
    }
 
      
    size_type rec_nnz(void_pointer p, size_type my_depth) const {
      size_type nn = 0;
      if (my_depth) {
	for (size_type k = 0; k < 16; ++k)
	  if (((void_pointer *)(p))[k])
	    nn += rec_nnz(((void_pointer *)(p))[k], my_depth-1);
      } else {
	for (size_type k = 0; k < 16; ++k)
	  if (((const T *)(p))[k] != T(0)) nn++;
      }
      return nn;
    }

    void copy_rec(void_pointer &p, const_void_pointer q, size_type my_depth) {
      if (my_depth) {
	p = new void_pointer[16];
	std::memset(p, 0, 16*sizeof(void_pointer));
	for (size_type l = 0; l < 16; ++l)
	  if (((const const_void_pointer *)(q))[l])
	    copy_rec(((void_pointer *)(p))[l],
		     ((const const_void_pointer *)(q))[l], my_depth-1);
      } else {
	p = new T[16];
	for (size_type l = 0; l < 16; ++l) ((T *)(p))[l] = ((const T *)(q))[l];
      }
    }

    void copy(const dsvector<T> &v) {
      if (root_ptr) rec_del(root_ptr, depth);
      root_ptr = 0;
      mask = v.mask; depth = v.depth; n = v.n; shift = v.shift;
      if (v.root_ptr) copy_rec(root_ptr, v.root_ptr, depth);
    }

    void next_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
		      const_pointer &pp, size_type &i, size_type base) const {
      size_type ii = i;
      if (my_depth) {
	my_mask = (my_mask >> 4);
	for (size_type k = 0; k < 16; ++k)
	  if (((void_pointer *)(p))[k] && (base + (k+1)*(my_mask+1)) >= i) {
	    next_pos_rec(((void_pointer *)(p))[k], my_depth-1, my_mask,
			 pp, i, base + k*(my_mask+1));
	    if (i != size_type(-1)) return; else i = ii;
	}
	i = size_type(-1); pp = 0;
      } else {
	for (size_type k = 0; k < 16; ++k)
	  if (base+k > i && ((const_pointer)(p))[k] != T(0))
	    { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
	i = size_type(-1); pp = 0;
      }
    }

    void previous_pos_rec(void_pointer p, size_type my_depth, size_type my_mask,
			  const_pointer &pp, size_type &i,
			  size_type base) const {
      size_type ii = i;
      if (my_depth) {
	my_mask = (my_mask >> 4);
	for (size_type k = 15; k != size_type(-1); --k)
	  if (((void_pointer *)(p))[k] && ((base + k*(my_mask+1)) < i)) {
	    previous_pos_rec(((void_pointer *)(p))[k], my_depth-1,
			     my_mask, pp, i, base + k*(my_mask+1));
	    if (i != size_type(-1)) return; else i = ii;
	}
	i = size_type(-1); pp = 0;
      } else {
	for (size_type k = 15; k != size_type(-1); --k)
	  if (base+k < i && ((const_pointer)(p))[k] != T(0))
	    { i = base+k; pp = &(((const_pointer)(p))[k]); return; }
	i = size_type(-1); pp = 0;
      }
    }
    
    
  public:
    void clean(double eps) { if (root_ptr) rec_clean(root_ptr, depth); }
    void resize(size_type n_) {
      if (n_ != n) {
	n = n_;
	if (n_ < n) { // Depth unchanged (a choice)
	  if (root_ptr) rec_clean_i(root_ptr, depth, mask, n_, 0);
	} else {
	  // may change the depth (add some levels)
	  size_type my_depth = 0, my_shift = 0, my_mask = 1; if (n_) --n_;
	  while (n_) { n_ /= 16; ++my_depth; my_shift += 4; my_mask *= 16; }
	  my_mask--; if (my_shift) my_shift -= 4; if (my_depth) --my_depth;
	  if (my_depth > depth || depth == 0) {
	    if (root_ptr) {
	      for (size_type k = depth; k < my_depth; ++k) {
		void_pointer *q = new void_pointer [16];
		std::memset(q, 0, 16*sizeof(void_pointer));
		q[0] = root_ptr; root_ptr = q;
	      }
	    }
	    mask = my_mask; depth = my_depth; shift = my_shift;
	  }
	}
      }
    }
    
    void clear(void) { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
    
    void next_pos(const_pointer &pp, size_type &i) const {
      if (!root_ptr || i >= n) { pp = 0, i = size_type(-1); return; }
      next_pos_rec(root_ptr, depth, mask, pp, i, 0);
    }

    void previous_pos(const_pointer &pp, size_type &i) const {
      if (!root_ptr) { pp = 0, i = size_type(-1); return; }
      if (i == size_type(-1)) { i = n; }
      previous_pos_rec(root_ptr, depth, mask, pp, i, 0);
    }

    iterator begin(void) {
      iterator it(*this); 
      if (n && root_ptr) {
	it.i = 0; it.p = const_cast<T *>(read_access(0));
	if (!(it.p) || *(it.p) == T(0))
	  next_pos(*(const_cast<const_pointer *>(&(it.p))), it.i);
      }
      return it;
    }

    iterator end(void) { return iterator(*this); }

    const_iterator begin(void) const {
      const_iterator it(*this);
      if (n && root_ptr) {
	it.i = 0; it.p = read_access(0);
	if (!(it.p) || *(it.p) == T(0)) next_pos(it.p, it.i);
      }
      return it;
    }

    const_iterator end(void) const { return const_iterator(*this); }
    
    inline ref_elt_vector<T, dsvector<T> > operator [](size_type c)
    { return ref_elt_vector<T, dsvector<T> >(this, c); }

    inline void w(size_type c, const T &e) {
      if (e == T(0)) { if (read_access(c)) *(write_access(c)) = e; }
      else *(write_access(c)) = e;
    }

    inline void wa(size_type c, const T &e)
    { if (e != T(0)) { *(write_access(c)) += e; } }

    inline T r(size_type c) const
    { const T *p = read_access(c); if (p) return *p; else return T(0); }

    inline T operator [](size_type c) const { return r(c); }
    
    size_type nnz(void) const
    { if (root_ptr) return rec_nnz(root_ptr, depth); else return 0; }
    size_type size(void) const { return n; }

    void swap(dsvector<T> &v) {
      std::swap(n, v.n); std::swap(root_ptr, v.root_ptr);
      std::swap(depth, v.depth); std::swap(shift, v.shift);
      std::swap(mask, v.mask);
    }
    
    /* Constructors */
    dsvector(const dsvector<T> &v) { init(0); copy(v); }
    dsvector<T> &operator =(const dsvector<T> &v) { copy(v); return *this; }
    explicit dsvector(size_type l){ init(l); }
    dsvector(void) { init(0); }
    ~dsvector() { if (root_ptr) rec_del(root_ptr, depth); root_ptr = 0; }
  };

  template <typename T> struct linalg_traits<dsvector<T>> {
    typedef dsvector<T> this_type;
    typedef this_type origin_type;
    typedef linalg_false is_reference;
    typedef abstract_vector linalg_type;
    typedef T value_type;
    typedef ref_elt_vector<T, dsvector<T> > reference;
    typedef dsvector_iterator<T>  iterator;
    typedef dsvector_const_iterator<T> const_iterator;
    typedef abstract_sparse storage_type;
    typedef linalg_true index_sorted;
    static size_type size(const this_type &v) { return v.size(); }
    static iterator begin(this_type &v) { return v.begin(); }
    static const_iterator begin(const this_type &v) { return v.begin(); }
    static iterator end(this_type &v) { return v.end(); }
    static const_iterator end(const this_type &v) { return v.end(); }
    static origin_type* origin(this_type &v) { return &v; }
    static const origin_type* origin(const this_type &v) { return &v; }
    static void clear(origin_type* o, const iterator &, const iterator &)
    { o->clear(); }
    static void do_clear(this_type &v) { v.clear(); }
    static value_type access(const origin_type *o, const const_iterator &,
			     const const_iterator &, size_type i)
    { return (*o)[i]; }
    static reference access(origin_type *o, const iterator &, const iterator &,
			    size_type i)
    { return (*o)[i]; }
    static void resize(this_type &v, size_type n) { v.resize(n); }
  };

  template<typename T> std::ostream &operator <<
  (std::ostream &o, const dsvector<T>& v) { gmm::write(o,v); return o; }

  /******* Optimized operations for dsvector<T> ****************************/

  template <typename T> inline void copy(const dsvector<T> &v1,
 					 dsvector<T> &v2) {
    GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
    v2 = v1;
  }
  template <typename T> inline void copy(const dsvector<T> &v1,
					 const dsvector<T> &v2) {
    GMM_ASSERT2(v1.size() == v2.size(), "dimensions mismatch");
    v2 = const_cast<dsvector<T> &>(v1);
  }
 template <typename T> inline
  void copy(const dsvector<T> &v1, const simple_vector_ref<dsvector<T> *> &v2){
    simple_vector_ref<dsvector<T> *>
      *svr = const_cast<simple_vector_ref<dsvector<T> *> *>(&v2);
    dsvector<T>
      *pv = const_cast<dsvector<T> *>((v2.origin));
    GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
    *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
  }
  template <typename T> inline
  void copy(const simple_vector_ref<const dsvector<T> *> &v1,
	    dsvector<T> &v2)
  { copy(*(v1.origin), v2); }
  template <typename T> inline
  void copy(const simple_vector_ref<dsvector<T> *> &v1, dsvector<T> &v2)
  { copy(*(v1.origin), v2); }
  template <typename T> inline
  void copy(const simple_vector_ref<dsvector<T> *> &v1,
	    const simple_vector_ref<dsvector<T> *> &v2)
  { copy(*(v1.origin), v2); }
  template <typename T> inline
  void copy(const simple_vector_ref<const dsvector<T> *> &v1,
	    const simple_vector_ref<dsvector<T> *> &v2)
  { copy(*(v1.origin), v2); }
  
  template <typename T>
  inline size_type nnz(const dsvector<T>& l) { return l.nnz(); }
  
  /*************************************************************************/
  /*                                                                       */
  /* Class wsvector: sparse vector optimized for random write operations,  */
  /* with log(n) complexity for read and write operations.                 */
  /* Based on std::map                                                     */
  /*                                                                       */
  /*************************************************************************/
  
  template<typename T> struct wsvector_iterator
    : public std::map<size_type, T>::iterator {
    typedef typename std::map<size_type, T>::iterator base_it_type;
    typedef T                   value_type;
    typedef value_type*         pointer;
    typedef value_type&         reference;
    // typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    
    reference operator *() const { return (base_it_type::operator*()).second; }
    pointer operator->() const { return &(operator*()); }
    size_type index(void) const { return (base_it_type::operator*()).first; }

    wsvector_iterator(void) {}
    wsvector_iterator(const base_it_type &it) : base_it_type(it) {}
  };

  template<typename T> struct wsvector_const_iterator
    : public std::map<size_type, T>::const_iterator {
    typedef typename std::map<size_type, T>::const_iterator base_it_type;
    typedef T                   value_type;
    typedef const value_type*   pointer;
    typedef const value_type&   reference;
    // typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    
    reference operator *() const { return (base_it_type::operator*()).second; }
    pointer operator->() const { return &(operator*()); }
    size_type index(void) const { return (base_it_type::operator*()).first; }

    wsvector_const_iterator(void) {}
    wsvector_const_iterator(const wsvector_iterator<T> &it)
      : base_it_type(it) {}
    wsvector_const_iterator(const base_it_type &it) : base_it_type(it) {}
  };


  /**
     sparse vector built upon std::map.
     Read and write access are quite fast (log n)
  */
  template<typename T> class wsvector : public std::map<size_type, T> {
  public:
    
    typedef typename std::map<int, T>::size_type size_type;
    typedef std::map<size_type, T> base_type;
    typedef typename base_type::iterator iterator;
    typedef typename base_type::const_iterator const_iterator;

  protected:
    size_type nbl;
    
  public:
    void clean(double eps);
    void resize(size_type);
    
    inline ref_elt_vector<T, wsvector<T> > operator [](size_type c)
    { return ref_elt_vector<T, wsvector<T> >(this, c); }

    inline void w(size_type c, const T &e) {
      GMM_ASSERT2(c < nbl, "out of range");
      if (e == T(0)) { this->erase(c); }
      else base_type::operator [](c) = e;
    }

    inline void wa(size_type c, const T &e) {
      GMM_ASSERT2(c < nbl, "out of range");
      if (e != T(0)) {
	iterator it = this->lower_bound(c);
	if (it != this->end() && it->first == c) it->second += e;
	else base_type::operator [](c) = e;
      }
    }

    inline T r(size_type c) const {
      GMM_ASSERT2(c < nbl, "out of range");
      const_iterator it = this->lower_bound(c);
      if (it != this->end() && c == it->first) return it->second;
      else return T(0);
    }

    inline T operator [](size_type c) const { return r(c); }
    
    size_type nb_stored(void) const { return base_type::size(); }
    size_type size(void) const { return nbl; }

    void swap(wsvector<T> &v)
    { std::swap(nbl, v.nbl); std::map<size_type, T>::swap(v); }
				       

    /* Constructors */
    void init(size_type l) { nbl = l; this->clear(); }
    explicit wsvector(size_type l){ init(l); }
    wsvector(void) { init(0); }
  };

  template<typename T>  void wsvector<T>::clean(double eps) {
    iterator it = this->begin(), itf = it, ite = this->end();
    while (it != ite) {
      ++itf; if (gmm::abs(it->second) <= eps) this->erase(it); it = itf;
    }
  }

  template<typename T>  void wsvector<T>::resize(size_type n) {
    if (n < nbl) {
      iterator it = this->begin(), itf = it, ite = this->end();
      while (it != ite) { ++itf; if (it->first >= n) this->erase(it); it=itf; }
    }
    nbl = n;
  }

  template <typename T> struct linalg_traits<wsvector<T> > {
    typedef wsvector<T> this_type;
    typedef this_type origin_type;
    typedef linalg_false is_reference;
    typedef abstract_vector linalg_type;
    typedef T value_type;
    typedef ref_elt_vector<T, wsvector<T> > reference;
    typedef wsvector_iterator<T>  iterator;
    typedef wsvector_const_iterator<T> const_iterator;
    typedef abstract_sparse storage_type;
    typedef linalg_true index_sorted;
    static size_type size(const this_type &v) { return v.size(); }
    static iterator begin(this_type &v) { return v.begin(); }
    static const_iterator begin(const this_type &v) { return v.begin(); }
    static iterator end(this_type &v) { return v.end(); }
    static const_iterator end(const this_type &v) { return v.end(); }
    static origin_type* origin(this_type &v) { return &v; }
    static const origin_type* origin(const this_type &v) { return &v; }
    static void clear(origin_type* o, const iterator &, const iterator &)
    { o->clear(); }
    static void do_clear(this_type &v) { v.clear(); }
    static value_type access(const origin_type *o, const const_iterator &,
			     const const_iterator &, size_type i)
    { return (*o)[i]; }
    static reference access(origin_type *o, const iterator &, const iterator &,
			    size_type i)
    { return (*o)[i]; }
    static void resize(this_type &v, size_type n) { v.resize(n); }
  };

  template<typename T> std::ostream &operator <<
  (std::ostream &o, const wsvector<T>& v) { gmm::write(o,v); return o; }

  /******* Optimized BLAS for wsvector<T> **********************************/

  template <typename T> inline void copy(const wsvector<T> &v1,
					 wsvector<T> &v2) {
    GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
    v2 = v1;
  }
  template <typename T> inline
  void copy(const wsvector<T> &v1, const simple_vector_ref<wsvector<T> *> &v2){
    simple_vector_ref<wsvector<T> *>
      *svr = const_cast<simple_vector_ref<wsvector<T> *> *>(&v2);
    wsvector<T>
      *pv = const_cast<wsvector<T> *>(v2.origin);
    GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
    *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
  }
  template <typename T> inline
  void copy(const simple_vector_ref<const wsvector<T> *> &v1,
	    wsvector<T> &v2)
  { copy(*(v1.origin), v2); }
  template <typename T> inline
  void copy(const simple_vector_ref<wsvector<T> *> &v1, wsvector<T> &v2)
  { copy(*(v1.origin), v2); }

  template <typename T> inline void clean(wsvector<T> &v, double eps) {
    typedef typename number_traits<T>::magnitude_type R;
    typename wsvector<T>::iterator it = v.begin(), ite = v.end(), itc;
    while (it != ite) 
      if (gmm::abs((*it).second) <= R(eps))
	{ itc=it; ++it; v.erase(itc); } else ++it; 
  }

  template <typename T>
  inline void clean(const simple_vector_ref<wsvector<T> *> &l, double eps) {
    simple_vector_ref<wsvector<T> *>
      *svr = const_cast<simple_vector_ref<wsvector<T> *> *>(&l);
    wsvector<T>
      *pv = const_cast<wsvector<T> *>((l.origin));
    clean(*pv, eps);
    svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
  }

  template <typename T>
  inline size_type nnz(const wsvector<T>& l) { return l.nb_stored(); }

  /*************************************************************************/
  /*                                                                       */
  /*    rsvector: sparse vector optimized for linear algebra operations.   */
  /*                                                                       */
  /*************************************************************************/

  template<typename T> struct elt_rsvector_ {
    size_type c; T e;
    /* e is initialized by default to avoid some false warnings of valgrind.
       (from http://valgrind.org/docs/manual/mc-manual.html:
      
       When memory is read into the CPU's floating point registers, the
       relevant V bits are read from memory and they are immediately
       checked. If any are invalid, an uninitialised value error is
       emitted. This precludes using the floating-point registers to copy
       possibly-uninitialised memory, but simplifies Valgrind in that it
       does not have to track the validity status of the floating-point
       registers.
    */
    elt_rsvector_(void) : e(0) {}
    elt_rsvector_(size_type cc) : c(cc), e(0) {}
    elt_rsvector_(size_type cc, const T &ee) : c(cc), e(ee) {}
    bool operator < (const elt_rsvector_ &a) const { return c < a.c; }
    bool operator == (const elt_rsvector_ &a) const { return c == a.c; }
    bool operator != (const elt_rsvector_ &a) const { return c != a.c; }
  };

  template<typename T> struct rsvector_iterator {
    typedef typename std::vector<elt_rsvector_<T> >::iterator IT;
    typedef T                   value_type;
    typedef value_type*         pointer;
    typedef value_type&         reference;
    typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::bidirectional_iterator_tag iterator_category;
    typedef rsvector_iterator<T> iterator;

    IT it;

    reference operator *() const { return it->e; }
    pointer operator->() const { return &(operator*()); }

    iterator &operator ++() { ++it; return *this; }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() { --it; return *this; }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &i) const { return it == i.it; }
    bool operator !=(const iterator &i) const { return !(i == *this); }

    size_type index(void) const { return it->c; }
    rsvector_iterator(void) {}
    rsvector_iterator(const IT &i) : it(i) {}
  };

  template<typename T> struct rsvector_const_iterator {
    typedef typename std::vector<elt_rsvector_<T> >::const_iterator IT;
    typedef T                   value_type;
    typedef const value_type*   pointer;
    typedef const value_type&   reference;
    typedef size_t              size_type;
    typedef ptrdiff_t           difference_type;
    typedef std::forward_iterator_tag iterator_category;
    typedef rsvector_const_iterator<T> iterator;

    IT it;

    reference operator *() const { return it->e; }
    pointer operator->() const { return &(operator*()); }
    size_type index(void) const { return it->c; }

    iterator &operator ++() { ++it; return *this; }
    iterator operator ++(int) { iterator tmp = *this; ++(*this); return tmp; }
    iterator &operator --() { --it; return *this; }
    iterator operator --(int) { iterator tmp = *this; --(*this); return tmp; }

    bool operator ==(const iterator &i) const { return it == i.it; }
    bool operator !=(const iterator &i) const { return !(i == *this); }

    rsvector_const_iterator(void) {}
    rsvector_const_iterator(const rsvector_iterator<T> &i) : it(i.it) {}
    rsvector_const_iterator(const IT &i) : it(i) {}
  };

  /**
     sparse vector built upon std::vector. Read access is fast,
     but insertion is O(n) 
  */
  template<typename T> class rsvector : public std::vector<elt_rsvector_<T> > {
  public:
    
    typedef std::vector<elt_rsvector_<T> > base_type_;
    typedef typename base_type_::iterator iterator;
    typedef typename base_type_::const_iterator const_iterator;
    typedef typename base_type_::size_type size_type;
    typedef T value_type;

  protected:
    size_type nbl;    	/* size of the vector.	          	  */
    
  public:

    void sup(size_type j);
    void base_resize(size_type n) { base_type_::resize(n); }
    void resize(size_type);
    
    ref_elt_vector<T, rsvector<T> > operator [](size_type c)
    { return ref_elt_vector<T, rsvector<T> >(this, c); }

    void w(size_type c, const T &e);
    void wa(size_type c, const T &e);
    T r(size_type c) const;
    void swap_indices(size_type i, size_type j);

    inline T operator [](size_type c) const { return r(c); }
    
    size_type nb_stored(void) const { return base_type_::size(); }
    size_type size(void) const { return nbl; }
    void clear(void) { base_type_::resize(0); }
    void swap(rsvector<T> &v)
    { std::swap(nbl, v.nbl); std::vector<elt_rsvector_<T> >::swap(v); }

    /* Constructeurs */
    explicit rsvector(size_type l) : nbl(l) { }
    rsvector(void) : nbl(0) { }
  };

  template <typename T>
  void rsvector<T>::swap_indices(size_type i, size_type j) {
    if (i > j) std::swap(i, j);
    if (i != j) {
      int situation = 0;
      elt_rsvector_<T> ei(i), ej(j), a;
      iterator it, ite, iti, itj;
      iti = std::lower_bound(this->begin(), this->end(), ei);
      if (iti != this->end() && iti->c == i) situation += 1;
      itj = std::lower_bound(this->begin(), this->end(), ej);
      if (itj != this->end() && itj->c == j) situation += 2;

      switch (situation) {
      case 1 : a = *iti; a.c = j; it = iti; ++it; ite = this->end();
	       for (; it != ite && it->c <= j; ++it, ++iti) *iti = *it;
	       *iti = a;
	       break;
      case 2 : a = *itj; a.c = i; it = itj; ite = this->begin();
	if (it != ite) {
	  --it;
	  while (it->c >= i) { *itj = *it;  --itj; if (it==ite) break; --it; }
	}
	*itj = a;
	break;
      case 3 : std::swap(iti->e, itj->e);
	       break;
      }
    }
  }

  template <typename T> void rsvector<T>::sup(size_type j) {
    if (nb_stored() != 0) {
      elt_rsvector_<T> ev(j);
      iterator it = std::lower_bound(this->begin(), this->end(), ev);
      if (it != this->end() && it->c == j) {
	for (iterator ite = this->end() - 1; it != ite; ++it) *it = *(it+1);
	base_resize(nb_stored()-1);
      }
    }
  }

  template<typename T>  void rsvector<T>::resize(size_type n) {
    if (n < nbl) {
      for (size_type i = 0; i < nb_stored(); ++i)
	if (base_type_::operator[](i).c >= n) { base_resize(i); break; }
    }
    nbl = n;
  }

  template <typename T> void rsvector<T>::w(size_type c, const T &e) {
    GMM_ASSERT2(c < nbl, "out of range");
    if (e == T(0)) sup(c);
    else {
      elt_rsvector_<T> ev(c, e);
      if (nb_stored() == 0) {
	base_type_::push_back(ev);
      }
      else {
	iterator it = std::lower_bound(this->begin(), this->end(), ev);
	if (it != this->end() && it->c == c) it->e = e;
	else {
	  size_type ind = it - this->begin(), nb = this->nb_stored();
          if (nb - ind > 1100)
            GMM_WARNING2("Inefficient addition of element in rsvector with "
                         << this->nb_stored() - ind << " non-zero entries");
	  base_type_::push_back(ev);
	  if (ind != nb) {
	    it = this->begin() + ind;
	    iterator ite = this->end(); --ite; iterator itee = ite; 
	    for (; ite != it; --ite) { --itee; *ite = *itee; }
	    *it = ev;
	  }
	}
      }
    }
  }

  template <typename T> void rsvector<T>::wa(size_type c, const T &e) {
    GMM_ASSERT2(c < nbl, "out of range");
    if (e != T(0)) {
      elt_rsvector_<T> ev(c, e);
      if (nb_stored() == 0) {
	base_type_::push_back(ev);
      }
      else {
	iterator it = std::lower_bound(this->begin(), this->end(), ev);
	if (it != this->end() && it->c == c) it->e += e;
	else {
	  size_type ind = it - this->begin(), nb = this->nb_stored();
          if (nb - ind > 1100)
            GMM_WARNING2("Inefficient addition of element in rsvector with "
                         << this->nb_stored() - ind << " non-zero entries");
	  base_type_::push_back(ev);
	  if (ind != nb) {
	    it = this->begin() + ind;
	    iterator ite = this->end(); --ite; iterator itee = ite; 
	    for (; ite != it; --ite) { --itee; *ite = *itee; }
	    *it = ev;
	  }
	}
      }
    }
  }
  
  template <typename T> T rsvector<T>::r(size_type c) const {
    GMM_ASSERT2(c < nbl, "out of range. Index " << c 
		<< " for a length of " << nbl);
    if (nb_stored() != 0) {
      elt_rsvector_<T> ev(c);
      const_iterator it = std::lower_bound(this->begin(), this->end(), ev);
      if (it != this->end() && it->c == c) return it->e;
    }
    return T(0);
  }

  template <typename T> struct linalg_traits<rsvector<T> > {
    typedef rsvector<T> this_type;
    typedef this_type origin_type;
    typedef linalg_false is_reference;
    typedef abstract_vector linalg_type;
    typedef T value_type;
    typedef ref_elt_vector<T, rsvector<T> > reference;
    typedef rsvector_iterator<T>  iterator;
    typedef rsvector_const_iterator<T> const_iterator;
    typedef abstract_sparse storage_type;
    typedef linalg_true index_sorted;
    static size_type size(const this_type &v) { return v.size(); }
    static iterator begin(this_type &v) { return iterator(v.begin()); }
    static const_iterator begin(const this_type &v)
    { return const_iterator(v.begin()); }
    static iterator end(this_type &v) { return iterator(v.end()); }
    static const_iterator end(const this_type &v)
      { return const_iterator(v.end()); }
    static origin_type* origin(this_type &v) { return &v; }
    static const origin_type* origin(const this_type &v) { return &v; }
    static void clear(origin_type* o, const iterator &, const iterator &)
    { o->clear(); }
    static void do_clear(this_type &v) { v.clear(); }
    static value_type access(const origin_type *o, const const_iterator &,
			     const const_iterator &, size_type i)
    { return (*o)[i]; }
    static reference access(origin_type *o, const iterator &, const iterator &,
			    size_type i)
    { return (*o)[i]; }
    static void resize(this_type &v, size_type n) { v.resize(n); }
  };

  template<typename T> std::ostream &operator <<
  (std::ostream &o, const rsvector<T>& v) { gmm::write(o,v); return o; }

  /******* Optimized operations for rsvector<T> ****************************/

  template <typename T> inline void copy(const rsvector<T> &v1,
 					 rsvector<T> &v2) {
    GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
    v2 = v1;
  }
  template <typename T> inline
  void copy(const rsvector<T> &v1, const simple_vector_ref<rsvector<T> *> &v2){
    simple_vector_ref<rsvector<T> *>
      *svr = const_cast<simple_vector_ref<rsvector<T> *> *>(&v2);
    rsvector<T>
      *pv = const_cast<rsvector<T> *>((v2.origin));
    GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
    *pv = v1; svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
  }
  template <typename T> inline
  void copy(const simple_vector_ref<const rsvector<T> *> &v1,
	    rsvector<T> &v2)
  { copy(*(v1.origin), v2); }
  template <typename T> inline
  void copy(const simple_vector_ref<rsvector<T> *> &v1, rsvector<T> &v2)
  { copy(*(v1.origin), v2); }

  template <typename V, typename T> inline void add(const V &v1,
						    rsvector<T> &v2) {
    if ((const void *)(&v1) != (const void *)(&v2)) {
      GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
	add_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
    }
  }

  template <typename V, typename T> 
  inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
  { add(v1, v2, abstract_dense(), abstract_sparse()); }

  template <typename V, typename T> 
  inline void add_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
  { add(v1, v2, abstract_skyline(), abstract_sparse()); }

  template <typename V, typename T> 
  void add_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
    add_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
  }

  template <typename V, typename T> 
  void add_rsvector(const V &v1, rsvector<T> &v2, linalg_false) {
    add(v1, v2, abstract_sparse(), abstract_sparse());
  }

  template <typename V, typename T> 
  void add_rsvector(const V &v1, rsvector<T> &v2, linalg_true) {
    typename linalg_traits<V>::const_iterator it1 = vect_const_begin(v1),
      ite1 = vect_const_end(v1);
    typename rsvector<T>::iterator it2 = v2.begin(), ite2 = v2.end(), it3;
    size_type nbc = 0, old_nbc = v2.nb_stored();
    for (; it1 != ite1 && it2 != ite2 ; ++nbc)
      if (it1.index() == it2->c) { ++it1; ++it2; }
      else if (it1.index() < it2->c) ++it1; else ++it2;
    for (; it1 != ite1; ++it1) ++nbc;
    for (; it2 != ite2; ++it2) ++nbc;

    v2.base_resize(nbc);
    it3 = v2.begin() + old_nbc;
    it2 = v2.end(); ite2 = v2.begin();
    it1 = vect_end(v1); ite1 = vect_const_begin(v1);
    while (it1 != ite1 && it3 != ite2) {
      --it3; --it1; --it2;
      if (it3->c > it1.index()) { *it2 = *it3; ++it1; }
      else if (it3->c == it1.index()) { *it2=*it3; it2->e+=*it1; }
      else { it2->c = it1.index(); it2->e = *it1; ++it3; }
    }
    while (it1 != ite1) { --it1; --it2; it2->c = it1.index(); it2->e = *it1; }
  }

  template <typename V, typename T> void copy(const V &v1, rsvector<T> &v2) {
    if ((const void *)(&v1) != (const void *)(&v2)) {
      GMM_ASSERT2(vect_size(v1) == vect_size(v2), "dimensions mismatch");
      if (same_origin(v1, v2))
	GMM_WARNING2("a conflict is possible in vector copy\n");
      copy_rsvector(v1, v2, typename linalg_traits<V>::storage_type());
    }
  }

  template <typename V, typename T> 
  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_dense)
  { copy_vect(v1, v2, abstract_dense(), abstract_sparse()); }

  template <typename V, typename T> 
  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_skyline)
  { copy_vect(v1, v2, abstract_skyline(), abstract_sparse()); }

  template <typename V, typename T>
  void copy_rsvector(const V &v1, rsvector<T> &v2, abstract_sparse) {
    copy_rsvector(v1, v2, typename linalg_traits<V>::index_sorted());
  }
  
  template <typename V, typename T2>
  void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_true) {
    typedef typename linalg_traits<V>::value_type T1;
    typename linalg_traits<V>::const_iterator it = vect_const_begin(v1),
      ite = vect_const_end(v1);
    v2.base_resize(nnz(v1));
    typename rsvector<T2>::iterator it2 = v2.begin();
    size_type nn = 0;
    for (; it != ite; ++it)
      if ((*it) != T1(0)) { it2->c = it.index(); it2->e = *it; ++it2; ++nn; }
    v2.base_resize(nn);
  }

  template <typename V, typename T2>
  void copy_rsvector(const V &v1, rsvector<T2> &v2, linalg_false) {
    typedef typename linalg_traits<V>::value_type T1;
    typename linalg_traits<V>::const_iterator it = vect_const_begin(v1),
      ite = vect_const_end(v1);
    v2.base_resize(nnz(v1));
    typename rsvector<T2>::iterator it2 = v2.begin();
    size_type nn = 0;
    for (; it != ite; ++it)
      if ((*it) != T1(0)) { it2->c = it.index(); it2->e = *it; ++it2; ++nn; }
    v2.base_resize(nn);
    std::sort(v2.begin(), v2.end());
  }
  
  template <typename T> inline void clean(rsvector<T> &v, double eps) {
    typedef typename number_traits<T>::magnitude_type R;
    typename rsvector<T>::iterator it = v.begin(), ite = v.end();
    for (; it != ite; ++it) if (gmm::abs((*it).e) <= eps) break;
    if (it != ite) {
      typename rsvector<T>::iterator itc = it;
      size_type erased = 1;
      for (++it; it != ite; ++it)
	{ *itc = *it; if (gmm::abs((*it).e) <= R(eps)) ++erased; else ++itc; }
      v.base_resize(v.nb_stored() - erased);
    }
  }

  template <typename T>
  inline void clean(const simple_vector_ref<rsvector<T> *> &l, double eps) {
    simple_vector_ref<rsvector<T> *>
      *svr = const_cast<simple_vector_ref<rsvector<T> *> *>(&l);
    rsvector<T>
      *pv = const_cast<rsvector<T> *>((l.origin));
    clean(*pv, eps);
    svr->begin_ = vect_begin(*pv); svr->end_ = vect_end(*pv);
  }
  
  template <typename T>
  inline size_type nnz(const rsvector<T>& l) { return l.nb_stored(); }

  /*************************************************************************/
  /*                                                                       */
  /* Class slvector: 'sky-line' vector.                                    */
  /*                                                                       */
  /*************************************************************************/

  template<typename T> struct slvector_iterator {
    typedef T value_type;
    typedef T *pointer;
    typedef T &reference;
    typedef ptrdiff_t difference_type;
    typedef std::random_access_iterator_tag iterator_category;
    typedef size_t size_type;
    typedef slvector_iterator<T> iterator;
    typedef typename std::vector<T>::iterator base_iterator;

    base_iterator it;
    size_type shift;
    
   
    iterator &operator ++()
    { ++it; ++shift; return *this; }
    iterator &operator --()
    { --it; --shift; return *this; }
    iterator operator ++(int)
    { iterator tmp = *this; ++(*(this)); return tmp; }
    iterator operator --(int)
    { iterator tmp = *this; --(*(this)); return tmp; }
    iterator &operator +=(difference_type i)
    { it += i; shift += i; return *this; }
    iterator &operator -=(difference_type i)
    { it -= i; shift -= i; return *this; }
    iterator operator +(difference_type i) const
    { iterator tmp = *this; return (tmp += i); }
    iterator operator -(difference_type i) const
    { iterator tmp = *this; return (tmp -= i); }
    difference_type operator -(const iterator &i) const
    { return it - i.it; }
	
    reference operator *() const
    { return *it; }
    reference operator [](int ii)
    { return *(it + ii); }
    
    bool operator ==(const iterator &i) const
    { return it == i.it; }
    bool operator !=(const iterator &i) const
    { return !(i == *this); }
    bool operator < (const iterator &i) const
    { return it < i.it; }
    size_type index(void) const { return shift; }

    slvector_iterator(void) {}
    slvector_iterator(const base_iterator &iter, size_type s)
      : it(iter), shift(s) {}
  };

  template<typename T> struct slvector_const_iterator {
    typedef T value_type;
    typedef const T *pointer;
    typedef value_type reference;
    typedef ptrdiff_t difference_type;
    typedef std::random_access_iterator_tag iterator_category;
    typedef size_t size_type;
    typedef slvector_const_iterator<T> iterator;
    typedef typename std::vector<T>::const_iterator base_iterator;

    base_iterator it;
    size_type shift;
    
   
    iterator &operator ++()
    { ++it; ++shift; return *this; }
    iterator &operator --()
    { --it; --shift; return *this; }
    iterator operator ++(int)
    { iterator tmp = *this; ++(*(this)); return tmp; }
    iterator operator --(int)
    { iterator tmp = *this; --(*(this)); return tmp; }
    iterator &operator +=(difference_type i)
    { it += i; shift += i; return *this; }
    iterator &operator -=(difference_type i)
    { it -= i; shift -= i; return *this; }
    iterator operator +(difference_type i) const
    { iterator tmp = *this; return (tmp += i); }
    iterator operator -(difference_type i) const
    { iterator tmp = *this; return (tmp -= i); }
    difference_type operator -(const iterator &i) const
    { return it - i.it; }
	
    value_type operator *() const
    { return *it; }
    value_type operator [](int ii)
    { return *(it + ii); }
    
    bool operator ==(const iterator &i) const
    { return it == i.it; }
    bool operator !=(const iterator &i) const
    { return !(i == *this); }
    bool operator < (const iterator &i) const
    { return it < i.it; }
    size_type index(void) const { return shift; }

    slvector_const_iterator(void) {}
    slvector_const_iterator(const slvector_iterator<T>& iter)
      : it(iter.it), shift(iter.shift) {}
    slvector_const_iterator(const base_iterator &iter, size_type s)
      : it(iter), shift(s) {}
  };


  /** skyline vector.
   */
  template <typename T> class slvector {
    
  public :
    typedef slvector_iterator<T> iterators;
    typedef slvector_const_iterator<T> const_iterators;
    typedef typename std::vector<T>::size_type size_type;
    typedef T value_type;

  protected :
    std::vector<T> data;
    size_type shift;
    size_type size_;


  public :

    size_type size(void) const { return size_; }
    size_type first(void) const { return shift; }
    size_type last(void) const { return shift + data.size(); }
    ref_elt_vector<T, slvector<T> > operator [](size_type c)
    { return ref_elt_vector<T, slvector<T> >(this, c); }

    typename std::vector<T>::iterator data_begin(void) { return data.begin(); }
    typename std::vector<T>::iterator data_end(void) { return data.end(); }
    typename std::vector<T>::const_iterator data_begin(void) const
      { return data.begin(); }
    typename std::vector<T>::const_iterator data_end(void) const
      { return data.end(); }

    void w(size_type c, const T &e);
    void wa(size_type c, const T &e);
    T r(size_type c) const {
      GMM_ASSERT2(c < size_, "out of range");
      if (c < shift || c >= shift + data.size()) return T(0);
      return data[c - shift];
    }

    inline T operator [](size_type c) const { return r(c); }
    void resize(size_type);
    void clear(void) { data.resize(0); shift = 0; }
    void swap(slvector<T> &v) {
      std::swap(data, v.data);
      std::swap(shift, v.shift);
      std::swap(size_, v.size_);
    }


    slvector(void) : data(0), shift(0), size_(0) {}
    explicit slvector(size_type l) : data(0), shift(0), size_(l) {}
    slvector(size_type l, size_type d, size_type s)
      : data(d), shift(s), size_(l) {}

  };

  template<typename T>  void slvector<T>::resize(size_type n) {
    if (n < last()) {
      if (shift >= n) clear(); else { data.resize(n-shift); }
    }
    size_ = n;
  }

  template<typename T>  void slvector<T>::w(size_type c, const T &e) {
    GMM_ASSERT2(c < size_, "out of range");
    size_type s = data.size();
    if (!s) { data.resize(1); shift = c; }
    else if (c < shift) {
      data.resize(s + shift - c); 
      typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
      typename std::vector<T>::iterator it3 = it2 - shift + c;
      for (; it3 >= it; --it3, --it2) *it2 = *it3;
      std::fill(it, it + shift - c, T(0));
      shift = c;
    }
    else if (c >= shift + s) {
      data.resize(c - shift + 1, T(0));
      // std::fill(data.begin() + s, data.end(), T(0));
    }
    data[c - shift] = e;
  }

  template<typename T>  void slvector<T>::wa(size_type c, const T &e) {
    GMM_ASSERT2(c < size_, "out of range");
    size_type s = data.size();
    if (!s) { data.resize(1, e); shift = c; return; }
    else if (c < shift) {
      data.resize(s + shift - c); 
      typename std::vector<T>::iterator it = data.begin(),it2=data.end()-1;
      typename std::vector<T>::iterator it3 = it2 - shift + c;
      for (; it3 >= it; --it3, --it2) *it2 = *it3;
      std::fill(it, it + shift - c, T(0));
      shift = c;
      data[c - shift] = e;
      return;
    }
    else if (c >= shift + s) {
      data.resize(c - shift + 1, T(0));
      data[c - shift] = e;
      return;
      // std::fill(data.begin() + s, data.end(), T(0));
    }
    data[c - shift] += e;
  }
  
  
  template <typename T> struct linalg_traits<slvector<T> > {
    typedef slvector<T> this_type;
    typedef this_type origin_type;
    typedef linalg_false is_reference;
    typedef abstract_vector linalg_type;
    typedef T value_type;
    typedef ref_elt_vector<T, slvector<T> > reference;
    typedef slvector_iterator<T>  iterator;
    typedef slvector_const_iterator<T> const_iterator;
    typedef abstract_skyline storage_type;
    typedef linalg_true index_sorted;
    static size_type size(const this_type &v) { return v.size(); }
    static iterator begin(this_type &v)
      { return iterator(v.data_begin(), v.first()); }
    static const_iterator begin(const this_type &v)
      { return const_iterator(v.data_begin(), v.first()); }
    static iterator end(this_type &v)
      { return iterator(v.data_end(), v.last()); }
    static const_iterator end(const this_type &v)
      { return const_iterator(v.data_end(), v.last()); }
    static origin_type* origin(this_type &v) { return &v; }
    static const origin_type* origin(const this_type &v) { return &v; }
    static void clear(origin_type* o, const iterator &, const iterator &)
    { o->clear(); }
    static void do_clear(this_type &v) { v.clear(); }
    static value_type access(const origin_type *o, const const_iterator &,
			     const const_iterator &, size_type i)
    { return (*o)[i]; }
    static reference access(origin_type *o, const iterator &, const iterator &,
			    size_type i)
    { return (*o)[i]; }
    static void resize(this_type &v, size_type n) { v.resize(n); }
  };

  template<typename T> std::ostream &operator <<
  (std::ostream &o, const slvector<T>& v) { gmm::write(o,v); return o; }

  template <typename T>
  inline size_type nnz(const slvector<T>& l) { return l.last() - l.first(); }

}

namespace std {
  template <typename T> void swap(gmm::wsvector<T> &v, gmm::wsvector<T> &w)
  { v.swap(w);}
  template <typename T> void swap(gmm::rsvector<T> &v, gmm::rsvector<T> &w)
  { v.swap(w);}
  template <typename T> void swap(gmm::slvector<T> &v, gmm::slvector<T> &w)
  { v.swap(w);}
}



#endif /* GMM_VECTOR_H__ */