/*!
 * \file
 * \brief Sparse Vector Class definitions
 * \author Tony Ottosson and Tobias Ringstrom
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2007  (see AUTHORS file for a list of contributors)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, 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 General Public License for more details.
 *
 * You should have received a copy of the GNU 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
 *
 * -------------------------------------------------------------------------
 */

#ifndef SVEC_H
#define SVEC_H

#ifndef _MSC_VER
#  include <itpp/config.h>
#else
#  include <itpp/config_msvc.h>
#endif

#include <itpp/base/vec.h>
#include <itpp/base/math/min_max.h>
#include <cstdlib>


namespace itpp {

  // Declaration of class Vec
  template <class T> class Vec;
  // Declaration of class Sparse_Vec
  template <class T> class Sparse_Vec;

  // ----------------------- Sparse_Vec Friends -------------------------------

  //! v1+v2 where v1 and v2 are sparse vector
  template <class T>
    Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

  //! v1*v2 where v1 and v2 are sparse vectors
  template <class T>
    T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

  //! v1*v2 where v1 is a sparse vector and v2 is a dense vector
  template <class T>
    T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2);

  //! v1*v2 where v1 is a dense vector and v2 is a sparse vector
  template <class T>
    T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2);

  //! Elementwise multiplication of two sparse vectors returning a sparse vector
  template <class T>
    Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

  //! Elementwise multiplication of a sparse vector and a dense vector returning a dense vector
  template <class T>
    Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2);

  //! Elementwise multiplication of a sparse vector and a dense vector returning a sparse vector
  template <class T>
    Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2);

  //! Elementwise multiplication of a dense vector and a sparse vector returning a dense vector
  template <class T>
    Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2);

  //! Elementwise multiplication of a dense vector and a sparse vector returning a sparse vector
  template <class T>
    Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2);

  /*!
    \brief Templated sparse vector class
    \author Tony Ottosson and Tobias Ringstrom

    A sparse vector is a vector where most elements are zero. The
    maximum number of none-zero elements is a parameter to the
    constructor. The elements are stored in random order, i.e. they
    are not sorted.

  */
  template <class T>
    class Sparse_Vec {
    public:

    //! Default constructor
    Sparse_Vec();

    /*!
      \brief Initiate an empty sparse vector

      \param sz Size of the sparse vector (i.e. maximum index is (\c sz - 1))
      \param data_init Maximum number of non-zero elements in the sparse vector (default value 200)
    */
    Sparse_Vec(int sz, int data_init=200);

    /*!
      \brief Initiate a new sparse vector.

      \param v The elements of \c v are copied into the new sparse vector
    */
    Sparse_Vec(const Sparse_Vec<T> &v);

    /*!
      \brief Initiate a new sparse vector from a dense vector.

      \param v The elements of \c v are copied into the new sparse vector
    */
    Sparse_Vec(const Vec<T> &v);

    /*!
      \brief Initiate a new sparse vector from a dense vector. Elements of \c v larger than \c epsilon are copied into the new sparse vector.

      \note If the type T is complex<double>, then the elements of \c v larger than \c abs(epsilon) are copied into the new sparse vector.
    */
    Sparse_Vec(const Vec<T> &v, T epsilon);

    //! Destructor
    ~Sparse_Vec();

    /*!
      \brief Set the size \c sz of the sparse vector. Default value \c data_init=-1 \c => allocated size for the data is not changed.

      \param sz Size of the sparse vector (i.e. maximum index is (\c sz - 1))
      \param data_init Maximum number of non-zero elements in the sparse vector (default value -1 \c => allocated size for the data is not changed)
    */
    void set_size(int sz, int data_init=-1);

    //! Returns the size of the sparse vector
    int size() const { return v_size; }

    //! Number of non-zero elements in the sparse vector
    inline int nnz()
      {
	if (check_small_elems_flag) {
	  remove_small_elements();
	}
	return used_size;
      }

    //! Returns the density of the sparse vector: (number of non-zero elements)/(size of the vector)
    double density();

    //! Set that all elements smaller than \a epsilon should be set to zero.
    void set_small_element(const T& epsilon);

    /*!
      Removes all elements that are smaller than \a epsilon from the non-zero elements.

      \note The small element \a epsilon can be set by the member function set_small_element. If no small value is set, the default value is always \c epsilon=0.
    */
    void remove_small_elements();

    /*!
      \brief Set the maximum number of non-zero elements to \c new_size

      \param new_size The new maximum number of non-zero elements.
    */
    void resize_data(int new_size);

    //! Set the maximum number of non-zero elements equal to the actual number of non-zero elements
    void compact();

    //! Returns a full, dense vector in \c v
    void full(Vec<T> &v) const;

    //! Returns a full, dense vector
    Vec<T> full() const;

    //! Returns the element with index \c i
    T operator()(int i) const;

    //! Set element \c i equal to \c v
    void set(int i, T v);

    //! Set the elements of the sparse vector with indices \c index_vec to the values in \c v
    void set(const ivec &index_vec, const Vec<T> &v);

    //! Set a new element with index \c i equal to \c v
    void set_new(int i, T v);

    //! Set new elements with indices \c index_vec equal to the values in \c v (no check whether the same index is used several times)
    void set_new(const ivec &index_vec, const Vec<T> &v);

    //! Add element \c i with \c v
    void add_elem(const int i, const T v);

    //! Add \c v to the elements specified by \c index_vec with \c v
    void add(const ivec &index_vec, const Vec<T> &v);

    //! Set the sparse vector to the all zero vector (removes all non-zero elements)
    void zeros();

    //! Set the i-th element to zero (i.e. clear that element if it contains a non-zero value)
    void zero_elem(const int i);

    //! Clear all non-zero elements of the sparse vector
    void clear();

    //! Clear the i-th element (if it contains a non-zero value)
    void clear_elem(const int i);

    /*!
      \brief Extract the reference to the p-th non-zero data element
    */
    inline void get_nz_data(int p, T& data_out)
      {
	if (check_small_elems_flag) {
	  remove_small_elements();
	}
	data_out = data[p];
      }

    //! Returns the p-th non-zero data element
    inline T get_nz_data(int p)
      {
	if (check_small_elems_flag) {
	  remove_small_elements();
	}
	return data[p];
      }

    //! Returns the vector index of the p-th non-zero element
    inline int get_nz_index(int p)
      {
	if (check_small_elems_flag) {
	  remove_small_elements();
	}
	return index[p];
      }

    //! Returns the p-th non-zero value in \c dat and the corresponding vector index in \c idx
    inline void get_nz(int p, int &idx, T &dat)
      {
	if (check_small_elems_flag) {
	  remove_small_elements();
	}
	idx=index[p];
	dat=data[p];
      }

    //! Retrun the indices of non-zero values
    ivec get_nz_indices();

    //! Return sparse subvector from index \c i1 to index \c i2
    Sparse_Vec<T> get_subvector(int i1, int i2) const;

    //! Returns the sum of all values squared
    T sqr() const;

    //! Assign sparse vector the value and length of the sparse vector \c v
    void operator=(const Sparse_Vec<T> &v);
    //! Assign sparse vector the value and length of the dense vector \c v
    void operator=(const Vec<T> &v);

   //! Returns the sign inverse of all elements in the sparse vector
    Sparse_Vec<T> operator-() const;

    //! Compare two sparse vectors. False if wrong sizes or different values
    bool operator==(const Sparse_Vec<T> &v);

    //! Add sparse vector \c v to all non-zero elements of the sparse vector
    void operator+=(const Sparse_Vec<T> &v);

    //! Add vector \c v to all non-zero elements of the sparse vector
    void operator+=(const Vec<T> &v);

    //! Subtract sparse vector \c v from all non-zero elements of the sparse vector
    void operator-=(const Sparse_Vec<T> &v);

    //! Subtract vector \c v from all non-zero elements of the sparse vector
    void operator-=(const Vec<T> &v);

    //! Multiply the scalar \c v to all non-zero elements of the sparse vector
    void operator*=(const T &v);

    //! Divide all non-zero elements of the sparse vector with the scalar \c v
    void operator/=(const T &v);

   //! Addition v1+v2 where v1 and v2 are sparse vector
    friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
    //! Scalar product v1*v2 where v1 and v2 are sparse vectors
    friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
    //! Scalar product v1*v2 where v1 is a sparse vector and v2 is a dense vector
    friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
    //! Scalar product v1*v2 where v1 is a dense vector and v2 is a sparse vector
    friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2);

    //! Element wise multiplication of two sparse vectors
    friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

    //! Element wise multiplication of a sparse vector and a dense vector
    friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);

    //! Element wise multiplication of a sparse vector and a dense vector returning a sparse vector
    friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);

    //! Element wise multiplication of a a dense vector and a sparse vector
    friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);

    //! Element wise multiplication of a a dense vector and a sparse vector returning a sparse vector
    friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);

    private:
    void init();
    void alloc();
    void free();

    int v_size, used_size, data_size;
    T *data;
    int *index;
    T eps;
    bool check_small_elems_flag;
  };


  /*!
    \relates Sparse_Vec
    \brief Type definition of an integer sparse vector
  */
  typedef Sparse_Vec<int> sparse_ivec;

  /*!
    \relates Sparse_Vec
    \brief Type definition of a double sparse vector
  */
  typedef Sparse_Vec<double> sparse_vec;

  /*!
    \relates Sparse_Vec
    \brief Type definition of a complex<double> sparse vector
  */
  typedef Sparse_Vec<std::complex<double> > sparse_cvec;

  // ----------------------- Implementation starts here --------------------------------

  template <class T>
    void Sparse_Vec<T>::init()
    {
      v_size = 0;
      used_size = 0;
      data_size = 0;
      data = 0;
      index = 0;
      eps = 0;
      check_small_elems_flag = true;
    }

  template <class T>
    void Sparse_Vec<T>::alloc()
    {
      if (data_size != 0) {
	data = new T[data_size];
	index = new int[data_size];
      }
    }

  template <class T>
    void Sparse_Vec<T>::free()
    {
      delete [] data;
      data = 0;
      delete [] index;
      index = 0;
    }

  template <class T>
    Sparse_Vec<T>::Sparse_Vec()
    {
      init();
    }

  template <class T>
    Sparse_Vec<T>::Sparse_Vec(int sz, int data_init)
    {
      init();
      v_size = sz;
      used_size = 0;
      data_size = data_init;
      alloc();
    }

  template <class T>
    Sparse_Vec<T>::Sparse_Vec(const Sparse_Vec<T> &v)
    {
      init();
      v_size = v.v_size;
      used_size = v.used_size;
      data_size = v.data_size;
      eps = v.eps;
      check_small_elems_flag = v.check_small_elems_flag;
      alloc();

      for (int i=0; i<used_size; i++) {
	data[i] = v.data[i];
	index[i] = v.index[i];
      }
    }

  template <class T>
    Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v)
    {
      init();
      v_size = v.size();
      used_size = 0;
      data_size = std::min(v.size(), 10000);
      alloc();

      for (int i=0; i<v_size; i++) {
	if (v(i) != T(0)) {
	  if (used_size == data_size)
	    resize_data(data_size*2);
	  data[used_size] = v(i);
	  index[used_size] = i;
	  used_size++;
	}
      }
      compact();
    }

  template <class T>
    Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon)
    {
      init();
      v_size = v.size();
      used_size = 0;
      data_size = std::min(v.size(), 10000);
      eps = epsilon;
      alloc();

    double e = std::abs(epsilon);
    for (int i=0; i<v_size; i++) {
	if (std::abs(v(i)) > e) {
	  if (used_size == data_size)
	    resize_data(data_size*2);
	  data[used_size] = v(i);
	  index[used_size] = i;
	  used_size++;
	}
      }
      compact();
    }

  template <class T>
    Sparse_Vec<T>::~Sparse_Vec()
    {
      free();
    }

  template <class T>
    void Sparse_Vec<T>::set_size(int new_size, int data_init)
    {
      v_size = new_size;
      used_size = 0;
      if (data_init != -1){
	free();
	data_size = data_init;
	alloc();
      }
    }

  template <class T>
    double Sparse_Vec<T>::density()
    {
      if (check_small_elems_flag) {
	remove_small_elements();
      }
      //return static_cast<double>(used_size) / v_size;
      return double(used_size) / v_size;
    }

  template <class T>
    void Sparse_Vec<T>::set_small_element(const T& epsilon)
    {
      eps=epsilon;
      remove_small_elements();
    }

  template <class T>
    void Sparse_Vec<T>::remove_small_elements()
    {
      int i;
      int nrof_removed_elements = 0;
      double e;

      //Remove small elements
      e = std::abs(eps);
      for (i=0;i<used_size;i++) {
	if (std::abs(data[i]) <= e) {
	  nrof_removed_elements++;
	}
	else if (nrof_removed_elements > 0) {
	  data[i-nrof_removed_elements] = data[i];
	  index[i-nrof_removed_elements] = index[i];
	}
      }

      //Set new size after small elements have been removed
      used_size -= nrof_removed_elements;

      //Set the flag to indicate that all small elements have been removed
      check_small_elems_flag = false;
    }


  template <class T>
    void Sparse_Vec<T>::resize_data(int new_size)
    {
      it_assert(new_size>=used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small");

      if (new_size != data_size) {
	if (new_size == 0)
	  free();
	else {
	  T *tmp_data = data;
	  int *tmp_pos = index;
	  data_size = new_size;
	  alloc();
	  for (int p=0; p<used_size; p++) {
	    data[p] = tmp_data[p];
	    index[p] = tmp_pos[p];
	  }
	  delete [] tmp_data;
	  delete [] tmp_pos;
	}
      }
    }

  template <class T>
    void Sparse_Vec<T>::compact()
    {
      if (check_small_elems_flag) {
	remove_small_elements();
      }
      resize_data(used_size);
    }

  template <class T>
    void Sparse_Vec<T>::full(Vec<T> &v) const
    {
      v.set_size(v_size);

      v = T(0);
      for (int p=0; p<used_size; p++)
	v(index[p]) = data[p];
    }

  template <class T>
    Vec<T> Sparse_Vec<T>::full() const
    {
      Vec<T> r(v_size);
      full(r);
      return r;
    }

  // This is slow. Implement a better search
  template <class T>
    T Sparse_Vec<T>::operator()(int i) const
    {
      it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");

      bool found = false;
      int p;
      for (p=0; p<used_size; p++) {
	if (index[p] == i) {
	  found = true;
	  break;
	}
      }
      return found ? data[p] : T(0);
    }

  template <class T>
    void Sparse_Vec<T>::set(int i, T v)
    {
      it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");

      bool found = false;
      bool larger_than_eps;
      int p;

      for (p=0; p<used_size; p++) {
	if (index[p] == i) {
	  found = true;
	  break;
	}
      }

      larger_than_eps = (std::abs(v) > std::abs(eps));

      if (found && larger_than_eps)
	data[p] = v;
      else if (larger_than_eps) {
	if (used_size == data_size)
	  resize_data(data_size*2+100);
	data[used_size] = v;
	index[used_size] = i;
	used_size++;
      }

      //Check if the stored element is smaller than eps. In that case it should be removed.
      if (std::abs(v) <= std::abs(eps)) {
	remove_small_elements();
      }

    }

  template <class T>
    void Sparse_Vec<T>::set_new(int i, T v)
    {
      it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");

      //Check that the new element is larger than eps!
      if (std::abs(v) > std::abs(eps)) {
	if (used_size == data_size)
	  resize_data(data_size*2+100);
	data[used_size] = v;
	index[used_size] = i;
	used_size++;
      }
    }

  template <class T>
    void Sparse_Vec<T>::add_elem(const int i, const T v)
    {
      bool found = false;
      int p;

      it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");

      for (p=0; p<used_size; p++) {
	if (index[p] == i) {
	  found = true;
	  break;
	}
      }
      if (found)
	data[p] += v;
      else {
	if (used_size == data_size)
	  resize_data(data_size*2+100);
	data[used_size] = v;
	index[used_size] = i;
	used_size++;
      }

      check_small_elems_flag = true;

    }

  template <class T>
    void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v)
    {
      bool found = false;
      int i,p,q;
      int nrof_nz=v.size();

      it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");

      //Elements are added if they have identical indices
      for (q=0; q<nrof_nz; q++){
	i=index_vec(q);
	for (p=0; p<used_size; p++) {
	  if (index[p] == i) {
	    found = true;
	    break;
	  }
	}
	if (found)
	  data[p] += v(q);
	else {
	  if (used_size == data_size)
	    resize_data(data_size*2+100);
	  data[used_size] = v(q);
	  index[used_size] = i;
	  used_size++;
	}
	found = false;
      }

      check_small_elems_flag = true;

    }

  template <class T>
    void Sparse_Vec<T>::zeros()
    {
      used_size = 0;
      check_small_elems_flag = false;
    }

  template <class T>
    void Sparse_Vec<T>::zero_elem(const int i)
    {
      bool found = false;
      int p;

      it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");

      for (p=0; p<used_size; p++) {
	if (index[p] == i) {
	  found = true;
	  break;
	}
      }
      if (found) {
	data[p] = data[used_size-1];
	index[p] = index[used_size-1];
	used_size--;
      }
    }

  template <class T>
    void Sparse_Vec<T>::clear()
    {
      used_size = 0;
      check_small_elems_flag = false;
    }

  template <class T>
    void Sparse_Vec<T>::clear_elem(const int i)
    {
      bool found = false;
      int p;

      it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");

      for (p=0; p<used_size; p++) {
	if (index[p] == i) {
	  found = true;
	  break;
	}
      }
      if (found) {
	data[p] = data[used_size-1];
	index[p] = index[used_size-1];
	used_size--;
      }
    }

  template <class T>
    void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v)
    {
      it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");

      //Clear all old non-zero elements
      clear();

      //Add the new non-zero elements
      add(index_vec,v);
    }

  template <class T>
    void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v)
    {
      int q;
      int nrof_nz=v.size();

      it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");

      //Clear all old non-zero elements
      clear();

      for (q=0; q<nrof_nz; q++){
	if (std::abs(v[q]) > std::abs(eps)) {
	  if (used_size == data_size)
	    resize_data(data_size*2+100);
	  data[used_size] = v(q);
	  index[used_size] = index_vec(q);
	  used_size++;
	}
      }
    }

  template <class T>
  ivec Sparse_Vec<T>::get_nz_indices()
  {
    int n = nnz();
    ivec r(n);
    for (int i = 0; i < n; i++) {
      r(i) = get_nz_index(i);
    }
    return r;
  }

  template <class T>
    Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const
    {
      it_assert_debug(v_size > i1 && v_size > i2 && i1<=i2 && i1>=0, "The index of the element exceeds the size of the sparse vector");

      Sparse_Vec<T> r(i2-i1+1);

      for (int p=0; p<used_size; p++) {
	if (index[p] >= i1 && index[p] <= i2) {
	  if (r.used_size == r.data_size)
	    r.resize_data(r.data_size*2+100);
	  r.data[r.used_size] = data[p];
	  r.index[r.used_size] = index[p]-i1;
	  r.used_size++;
	}
      }
      r.eps = eps;
      r.check_small_elems_flag = check_small_elems_flag;
      r.compact();

      return r;
    }

  template <class T>
    T Sparse_Vec<T>::sqr() const
    {
      T sum(0);
      for (int p=0; p<used_size; p++)
	sum += data[p] * data[p];

      return sum;
    }

  template <class T>
    void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v)
    {
      free();
      v_size = v.v_size;
      used_size = v.used_size;
      data_size = v.data_size;
      eps = v.eps;
      check_small_elems_flag = v.check_small_elems_flag;
      alloc();

      for (int i=0; i<used_size; i++) {
	data[i] = v.data[i];
	index[i] = v.index[i];
      }
    }

  template <class T>
    void Sparse_Vec<T>::operator=(const Vec<T> &v)
    {
      free();
      v_size = v.size();
      used_size = 0;
      data_size = std::min(v.size(), 10000);
      eps = T(0);
      check_small_elems_flag = false;
      alloc();

      for (int i=0; i<v_size; i++) {
	if (v(i) != T(0)) {
	  if (used_size == data_size)
	    resize_data(data_size*2);
	  data[used_size] = v(i);
	  index[used_size] = i;
	  used_size++;
	}
      }
      compact();
    }

  template <class T>
    Sparse_Vec<T> Sparse_Vec<T>::operator-() const
    {
      Sparse_Vec r(v_size, used_size);

      for (int p=0; p<used_size; p++) {
	r.data[p] = -data[p];
	r.index[p] = index[p];
      }
      r.used_size = used_size;

      return r;
    }

  template <class T>
    bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v)
    {
      int p,q;
      bool found=false;

      //Remove small elements before comparing the two sparse_vectors
      if (check_small_elems_flag)
	remove_small_elements();

      if (v_size!=v.v_size) {
	//Return false if vector sizes are unequal
	return false;
      } else {
	for(p=0;p<used_size;p++) {
	  for(q=0;q<v.used_size;q++) {
	    if (index[p] == v.index[q]) {
	      found = true;
	      break;
	    }
	  }
	  if (found==false)
	    //Return false if non-zero element not found, or if elements are unequal
	    return false;
	  else if (data[p]!=v.data[q])
	    //Return false if non-zero element not found, or if elements are unequal
	    return false;
	  else
	    //Check next non-zero element
	    found = false;
	}
      }

      /*Special handling if sizes do not match.
	Required since v may need to do remove_small_elements() for true comparison*/
      if (used_size!=v.used_size) {
	if (used_size > v.used_size) {
	  //Return false if number of non-zero elements is less in v
	  return false;
	}
	else {
	  //Ensure that the remaining non-zero elements in v are smaller than v.eps
	  int nrof_small_elems = 0;
	  for(q=0;q<v.used_size;q++) {
	    if (std::abs(v.data[q]) <= std::abs(v.eps))
	      nrof_small_elems++;
	  }
	  if (v.used_size-nrof_small_elems != used_size)
	    //Return false if the number of "true" non-zero elements are unequal
	    return false;
	}
      }

      //All elements checks => return true
      return true;
    }

  template <class T>
    void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v)
    {
      int i,p;
      T tmp_data;
      int nrof_nz_v=v.used_size;

      it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");

      for (p=0; p<nrof_nz_v; p++) {
	i = v.index[p];
	tmp_data = v.data[p];
	//get_nz(p,i,tmp_data);
	add_elem(i,tmp_data);
      }

      check_small_elems_flag = true;
    }

  template <class T>
    void Sparse_Vec<T>::operator+=(const Vec<T> &v)
    {
      int i;

      it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");

      for (i=0; i<v.size(); i++)
	if (v(i)!=T(0))
	  add_elem(i,v(i));

      check_small_elems_flag = true;
    }


  template <class T>
    void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v)
    {
      int i,p;
      T tmp_data;
      int nrof_nz_v=v.used_size;

      it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");

      for (p=0; p<nrof_nz_v; p++) {
	i = v.index[p];
	tmp_data = v.data[p];
	//v.get_nz(p,i,tmp_data);
	add_elem(i,-tmp_data);
      }

      check_small_elems_flag = true;
    }

  template <class T>
    void Sparse_Vec<T>::operator-=(const Vec<T> &v)
    {
      int i;

      it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");

      for (i=0; i<v.size(); i++)
	if (v(i)!=T(0))
	  add_elem(i,-v(i));

      check_small_elems_flag = true;
    }

  template <class T>
    void Sparse_Vec<T>::operator*=(const T &v)
    {
      int p;

      for (p=0; p<used_size; p++) {
	data[p]*=v;}

      check_small_elems_flag = true;
    }

  template <class T>
    void Sparse_Vec<T>::operator/=(const T &v)
    {
      int p;
      for (p=0; p<used_size; p++) {
	data[p]/=v;}

      if (std::abs(eps) > 0) {
	check_small_elems_flag = true;
      }
    }

  template <class T>
    T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.v_size == v2.v_size, "Sparse_Vec<T> * Sparse_Vec<T>");

      T sum(0);
      Vec<T> v1f(v1.v_size);
      v1.full(v1f);
      for (int p=0; p<v2.used_size; p++) {
	if (v1f[v2.index[p]] != T(0))
	  sum += v1f[v2.index[p]] * v2.data[p];
      }

      return sum;
    }

  template <class T>
    T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2)
    {
      it_assert_debug(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");

      T sum(0);
      for (int p1=0; p1<v1.used_size; p1++)
	sum += v1.data[p1] * v2[v1.index[p1]];

      return sum;
    }

  template <class T>
    T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted");

      T sum(0);
      for (int p2=0; p2<v2.used_size; p2++)
	sum += v1[v2.index[p2]] * v2.data[p2];

      return sum;
    }

  template <class T>
    Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.v_size == v2.v_size, "elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)");

      Sparse_Vec<T> r(v1.v_size);
      ivec pos(v1.v_size);
      pos = -1;
      for (int p1=0; p1<v1.used_size; p1++)
	pos[v1.index[p1]] = p1;
      for (int p2=0; p2<v2.used_size; p2++) {
	if (pos[v2.index[p2]] != -1) {
	  if (r.used_size == r.data_size)
	    r.resize_data(r.used_size*2+100);
	  r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2];
	  r.index[r.used_size] = v2.index[p2];
	  r.used_size++;
	}
      }
      r.compact();

      return r;
    }

  template <class T>
    Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2)
    {
      it_assert_debug(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");

      Vec<T> r(v1.v_size);
      r = T(0);
      for (int p1=0; p1<v1.used_size; p1++)
	r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]];

      return r;
    }

  template <class T>
    Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2)
    {
      it_assert_debug(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)");

      Sparse_Vec<T> r(v1.v_size);
      for (int p1=0; p1<v1.used_size; p1++) {
	if (v2[v1.index[p1]] != T(0)) {
	  if (r.used_size == r.data_size)
	    r.resize_data(r.used_size*2+100);
	  r.data[r.used_size] = v1.data[p1] * v2[v1.index[p1]];
	  r.index[r.used_size] = v1.index[p1];
	  r.used_size++;
	}
      }
      r.compact();

      return r;
    }

  template <class T>
    Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");

      Vec<T> r(v2.v_size);
      r = T(0);
      for (int p2=0; p2<v2.used_size; p2++)
	r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2];

      return r;
    }

  template <class T>
    Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)");

      Sparse_Vec<T> r(v2.v_size);
      for (int p2=0; p2<v2.used_size; p2++) {
	if (v1[v2.index[p2]] != T(0)) {
	  if (r.used_size == r.data_size)
	    r.resize_data(r.used_size*2+100);
	  r.data[r.used_size] = v1[v2.index[p2]] * v2.data[p2];
	  r.index[r.used_size] = v2.index[p2];
	  r.used_size++;
	}
      }
      r.compact();

      return r;
    }

  template <class T>
    Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2)
    {
      it_assert_debug(v1.v_size == v2.v_size, "Sparse_Vec<T> + Sparse_Vec<T>");

      Sparse_Vec<T> r(v1);
      ivec pos(v1.v_size);
      pos = -1;
      for (int p1=0; p1<v1.used_size; p1++)
	pos[v1.index[p1]] = p1;
      for (int p2=0; p2<v2.used_size; p2++) {
	if (pos[v2.index[p2]] == -1) {// A new entry
	  if (r.used_size == r.data_size)
	    r.resize_data(r.used_size*2+100);
	  r.data[r.used_size] = v2.data[p2];
	  r.index[r.used_size] = v2.index[p2];
	  r.used_size++;
	}
	else
	  r.data[pos[v2.index[p2]]] += v2.data[p2];
      }
      r.check_small_elems_flag = true;  // added dec 7, 2006
      r.compact();

      return r;
    }

  //! Convert a dense vector \c v into its sparse representation
  template <class T>
    inline Sparse_Vec<T> sparse(const Vec<T> &v)
    {
      Sparse_Vec<T> s(v);
      return s;
    }

  //! Convert a dense vector \c v into its sparse representation
  template <class T>
    inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon)
    {
      Sparse_Vec<T> s(v, epsilon);
      return s;
    }

  //! Convert a sparse vector \c s into its dense representation
  template <class T>
    inline Vec<T> full(const Sparse_Vec<T> &s)
    {
      Vec<T> v;
      s.full(v);
      return v;
    }

  //! \cond

  // ---------------------------------------------------------------------
  // Instantiations
  // ---------------------------------------------------------------------

#ifdef HAVE_EXTERN_TEMPLATE

  extern template class Sparse_Vec<int>;
  extern template class Sparse_Vec<double>;
  extern template class Sparse_Vec<std::complex<double> >;

  extern template sparse_ivec operator+(const sparse_ivec &,
                                        const sparse_ivec &);
  extern template sparse_vec operator+(const sparse_vec &,
                                       const sparse_vec &);
  extern template sparse_cvec operator+(const sparse_cvec &,
                                        const sparse_cvec &);

  extern template int operator*(const sparse_ivec &, const sparse_ivec &);
  extern template double operator*(const sparse_vec &, const sparse_vec &);
  extern template std::complex<double> operator*(const sparse_cvec &,
                                                 const sparse_cvec &);

  extern template int operator*(const sparse_ivec &, const ivec &);
  extern template double operator*(const sparse_vec &, const vec &);
  extern template std::complex<double> operator*(const sparse_cvec &,
                                                 const cvec &);

  extern template int operator*(const ivec &, const sparse_ivec &);
  extern template double operator*(const vec &, const sparse_vec &);
  extern template std::complex<double> operator*(const cvec &,
                                                 const sparse_cvec &);

  extern template sparse_ivec elem_mult(const sparse_ivec &,
                                        const sparse_ivec &);
  extern template sparse_vec elem_mult(const sparse_vec &, const sparse_vec &);
  extern template sparse_cvec elem_mult(const sparse_cvec &,
                                        const sparse_cvec &);

  extern template ivec elem_mult(const sparse_ivec &, const ivec &);
  extern template vec elem_mult(const sparse_vec &, const vec &);
  extern template cvec elem_mult(const sparse_cvec &, const cvec &);

  extern template sparse_ivec elem_mult_s(const sparse_ivec &, const ivec &);
  extern template sparse_vec elem_mult_s(const sparse_vec &, const vec &);
  extern template sparse_cvec elem_mult_s(const sparse_cvec &, const cvec &);

  extern template ivec elem_mult(const ivec &, const sparse_ivec &);
  extern template vec elem_mult(const vec &, const sparse_vec &);
  extern template cvec elem_mult(const cvec &, const sparse_cvec &);

  extern template sparse_ivec elem_mult_s(const ivec &, const sparse_ivec &);
  extern template sparse_vec elem_mult_s(const vec &, const sparse_vec &);
  extern template sparse_cvec elem_mult_s(const cvec &, const sparse_cvec &);

#endif // HAVE_EXTERN_TEMPLATE

  //! \endcond

} // namespace itpp

#endif // #ifndef SVEC_H