svec.h

Go to the documentation of this file.

#ifndef SVEC_H
#define SVEC_H

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

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 -------------------------------

template <class T>
Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

template <class T>
T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

template <class T>
T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2);

template <class T>
T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2);

template <class T>
Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

template <class T>
Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2);

template <class T>
Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2);

template <class T>
Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2);

template <class T>
Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2);

namespace details
{
        //this template selects appropriate type for Eps value used to remove small elements from
        //sparse containers
        template <typename NumT> struct Sparse_Eps_Type_Selector;
        template <> struct Sparse_Eps_Type_Selector<double> {typedef double eps_type;};
        template <> struct Sparse_Eps_Type_Selector<std::complex<double> > {typedef double eps_type;};
        template <> struct Sparse_Eps_Type_Selector<int> {typedef int eps_type;};
        template <> struct Sparse_Eps_Type_Selector<short> {typedef short eps_type;};
        template <> struct Sparse_Eps_Type_Selector<itpp::bin> {typedef int eps_type;};
}


template <class T>
class Sparse_Vec
{
public:

  Sparse_Vec();

  Sparse_Vec(int sz, int data_init = 200);

  Sparse_Vec(const Sparse_Vec<T> &v);

  Sparse_Vec(const Vec<T> &v);

  Sparse_Vec(const Vec<T> &v, T epsilon);

  ~Sparse_Vec();

  void set_size(int sz, int data_init = -1);

  int size() const { return v_size; }

  inline int nnz() {
    if (check_small_elems_flag) {
      remove_small_elements();
    }
    return used_size;
  }

  double density();

  void set_small_element(const T& epsilon);

  void remove_small_elements();

  void resize_data(int new_size);

  void compact();

  void full(Vec<T> &v) const;

  Vec<T> full() const;

  T operator()(int i) const;

  void set(int i, T v);

  void set(const ivec &index_vec, const Vec<T> &v);

  void set_new(int i, T v);

  void set_new(const ivec &index_vec, const Vec<T> &v);

  void add_elem(const int i, const T v);

  void add(const ivec &index_vec, const Vec<T> &v);

  void zeros();

  void zero_elem(const int i);

  void clear();

  void clear_elem(const int i);

  inline void get_nz_data(int p, T& data_out) {
    if (check_small_elems_flag) {
      remove_small_elements();
    }
    data_out = data[p];
  }

  inline T get_nz_data(int p) {
    if (check_small_elems_flag) {
      remove_small_elements();
    }
    return data[p];
  }

  inline int get_nz_index(int p) {
    if (check_small_elems_flag) {
      remove_small_elements();
    }
    return index[p];
  }

  inline void get_nz(int p, int &idx, T &dat) {
    if (check_small_elems_flag) {
      remove_small_elements();
    }
    idx = index[p];
    dat = data[p];
  }

  ivec get_nz_indices();

  Sparse_Vec<T> get_subvector(int i1, int i2) const;

  T sqr() const;

  void operator=(const Sparse_Vec<T> &v);
  void operator=(const Vec<T> &v);

  Sparse_Vec<T> operator-() const;

  bool operator==(const Sparse_Vec<T> &v);

  void operator+=(const Sparse_Vec<T> &v);

  void operator+=(const Vec<T> &v);

  void operator-=(const Sparse_Vec<T> &v);

  void operator-=(const Vec<T> &v);

  void operator*=(const T &v);

  void operator/=(const T &v);

  friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
  friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
  friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
  friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2);

  friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);

  friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);

  friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);

  friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);

  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;
  typename details::Sparse_Eps_Type_Selector<T>::eps_type eps;
  bool check_small_elems_flag;
};


typedef Sparse_Vec<int> sparse_ivec;

typedef Sparse_Vec<double> sparse_vec;

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 = std::abs(epsilon);
  alloc();

  for (int i = 0; i < v_size; i++) {
    if (std::abs(v(i)) > eps) {
      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 = std::abs(epsilon);
  remove_small_elements();
}

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

  //Remove small elements
  for (i = 0;i < used_size;i++) {
    if (std::abs(data[i]) <= eps) {
      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) > 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) <= 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) > 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]) > 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 = std::abs(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]) <= 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 (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;
}

template <class T>
inline Sparse_Vec<T> sparse(const Vec<T> &v)
{
  Sparse_Vec<T> s(v);
  return s;
}

template <class T>
inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon)
{
  Sparse_Vec<T> s(v, epsilon);
  return s;
}

template <class T>
inline Vec<T> full(const Sparse_Vec<T> &s)
{
  Vec<T> v;
  s.full(v);
  return v;
}


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

ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Vec<int>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Vec<double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Sparse_Vec<std::complex<double> >;

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

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

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

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

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_ivec elem_mult(const sparse_ivec &,
                                        const sparse_ivec &);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_vec elem_mult(const sparse_vec &, const sparse_vec &);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT sparse_cvec elem_mult(const sparse_cvec &,
                                        const sparse_cvec &);

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

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

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

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


} // namespace itpp

#endif // #ifndef SVEC_H