/Users/mourrain/Devel/mmx/basix/include/basix/sparse_vector.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : sparse_vector.hpp
* DESCRIPTION: Hash sparse_vectors
* COPYRIGHT  : (C) 2004  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_SPARSE_VECTOR_HPP
#define __MMX_SPARSE_VECTOR_HPP
#include <basix/pair.hpp>
#include <basix/symbol.hpp>


namespace mmx {
struct eq_sparse_vector;
#define TMPL_DEF \
  template<typename C, typename T, typename V= exact_eq_sparse_vector>
#define TMPL template<typename C, typename T, typename V>
#define Sparse_vector sparse_vector<C,T,V>
#define Sparse_vector_rep sparse_vector_rep<C,T,V>
#define Pair pair<T,C>
TMPL class sparse_vector_rep;
TMPL class sparse_vector;
TMPL inline nat N (const Sparse_vector& v);

/******************************************************************************
* Variants of sparse vectors
******************************************************************************/

struct hard_eq_sparse_vector {
  typedef hard_eq_op key_op;
  typedef equal_op val_op;
  typedef hard_less_op l_op;
};

struct exact_eq_sparse_vector {
  typedef exact_eq_op key_op;
  typedef equal_op val_op;
  typedef less_op l_op;
};

/******************************************************************************
* Sparse_vector class and its representation class
******************************************************************************/

TMPL_DEF
class sparse_vector_rep REP_STRUCT {
  Pair* a;   // the entries
  nat   n;   // number of entries
  nat   l;   // allocated number of entries

public:
  inline sparse_vector_rep (Pair* a2, nat n2):
    a (a2), n (n2), l (n2) {}
  inline sparse_vector_rep (Pair* a2, nat n2, nat l2):
    a (a2), n (n2), l (l2) {}
  inline ~sparse_vector_rep () {
    mmx_delete<Pair > (a, l); }

  friend class Sparse_vector;
  friend nat N LESSGTR (const Sparse_vector& v);
};

TMPL_DEF
class sparse_vector {
INDIRECT_PROTO_3 (sparse_vector, sparse_vector_rep, C, T, V)
  inline sparse_vector () {
    rep= new Sparse_vector_rep (mmx_new<Pair > (0), 0); }
  inline sparse_vector (const T& k, const C& v) {
    typedef typename V::val_op Eq;
    if (Eq::op (v, C (0))) rep= new Sparse_vector_rep (mmx_new<Pair > (0), 0);
    else rep= new Sparse_vector_rep (mmx_new<Pair > (1, Pair (k, v)), 1); }
  inline sparse_vector (const Pair& p) {
    rep= new Sparse_vector_rep (mmx_new<Pair > (1, p), 1); }
  inline sparse_vector (Pair* a, nat n) {
    rep= new Sparse_vector_rep (a, n); }
  inline sparse_vector (Pair* a, nat n, nat l) {
    rep= new Sparse_vector_rep (a, n, l); }
  inline Pair operator [] (nat i) const {
    VERIFY (i<rep->n, "index out of range");
    return rep->a[i]; }
  inline Pair& operator [] (nat i) {
    VERIFY (i<rep->n, "index out of range");
    secure ();
    return rep->a[i]; }
  friend nat N LESSGTR (const Sparse_vector& v);
};
INDIRECT_IMPL_3 (sparse_vector, sparse_vector_rep,
                 typename C, C, typename T, T, typename V, V)

TMPL inline nat N (const Sparse_vector& v) { return v.rep->n; }

/******************************************************************************
* Printing sparse_vectors
******************************************************************************/

TMPL syntactic
flatten (const Sparse_vector& v) {
  vector<syntactic> r;
  for (nat i=0; i<N(v); i++)
    r << flatten (v[i]);
  return apply ("sparse_vector", r);
}

TMPL
struct binary_helper<Sparse_vector >:
  public void_binary_helper<Sparse_vector >
{
  static inline string short_type_name () {
    return "Sv" * Short_type_name (C); }
  static inline generic full_type_name () {
    return gen ("Sparse_vector", Full_type_name (C)); }
  static inline nat size (const Sparse_vector& v) {
    return N(v); }
  static inline generic access (const Sparse_vector& v, nat i) {
    return as<generic> (v[i]); }
  static inline generic disassemble (const Sparse_vector& v) {
    vector<generic> r;
    for (nat i=0; i<N(v); i++) r << as<generic> (v[i]);
    return as<generic> (r); }
  static inline Sparse_vector assemble (const generic& x) {
    vector<generic> v= as<vector<generic> > (x);
    Pair* r= mmx_new<Pair > (N(v));
    for (nat i=0; i<N(v); i++) r[i]= as<Pair > (v);
    return Sparse_vector (r, N(v)); }
  static inline void write (const port& out, const Sparse_vector& v) {
    binary_write<nat> (out, N(v));
    for (nat i=0; i<N(v); i++)
      binary_write<Pair > (out, v[i]); }
  static Sparse_vector read (const port& in, nat n) {
    if (n == 0) return Sparse_vector ();
    else if (n == 1) return Sparse_vector (binary_read<Pair > (in));
    else return read (in, (n+1) >> 1) + read (in, n >> 1); }
  static inline Sparse_vector read (const port& in) {
    nat n= binary_read<nat> (in);
    return read (in, n); }
};

/******************************************************************************
* Hashing and equality testing
******************************************************************************/

template<typename Op, typename C, typename T, typename V> nat
unary_hash (const Sparse_vector& v) {
  nat h=1236321;
  for (nat i=0; i<N(v); i++)
    h = (h << 5) ^ (h >> 27) ^ Op::op (v[i]);
  return h;
}

template<typename Op, typename C, typename T, typename V> bool
binary_test (const Sparse_vector& v1, const Sparse_vector& v2) {
  nat n1= N(v1), n2= N(v2);
  if (n1 != n2) return false;
  for (nat i=0; i<n1; i++)
    if (!Op::op (v1[i], v2[i]))
      return false;
  return true;
}

TRUE_IDENTITY_OP_SUGAR(TMPL,Sparse_vector)
EXACT_IDENTITY_OP_SUGAR(TMPL,Sparse_vector)

/******************************************************************************
* Unary operations
******************************************************************************/

template<typename Op, typename C, typename T, typename V> Sparse_vector
unary_map (const Sparse_vector& v) {
  typedef typename V::val_op Eq;
  nat i, j, n= N(v);
  Pair* r= mmx_new<Pair > (n);
  for (i= j= 0; i < n; i++) {
    Pair e (v[i].x1, Op::op (v[i].x2));
    if (Eq::not_op (e.x2, C (0))) r[j++]= e;
  }
  return Sparse_vector (r, j, n);
}

TMPL inline Sparse_vector
operator - (const Sparse_vector& t) {
  return unary_map<neg_op,C,T,V> (t); }

/******************************************************************************
* Binary operations
******************************************************************************/

template<typename Op, typename C, typename T, typename V> Sparse_vector
binary_map (const Sparse_vector& v1, const Sparse_vector& v2) {
  typedef typename V::val_op Eq;
  typedef typename V::l_op Less;
  nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
  Pair* r= mmx_new<Pair > (n);
  for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
    if (Less::op (v1[i1].x1, v2[i2].x1)) {
      Pair e (v1[i1].x1, Op::op (v1[i1].x2, C(0)));
      i1++; if (Eq::not_op (e.x2, C(0))) r[j++]= e;
    }
    else if (Less::op (v2[i2].x1, v1[i1].x1)) {
      Pair e (v2[i2].x1, Op::op (C(0), v2[i2].x2));
      i2++; if (Eq::not_op (e.x2, C(0))) r[j++]= e;
    }
    else {
      Pair e (v2[i2].x1, Op::op (v1[i1].x2, v2[i2].x2));
      i1++; i2++; if (Eq::not_op (e.x2, C (0))) r[j++]= e;
    }
  }
  while (i1 < n1) {
    Pair e (v1[i1].x1, Op::op (v1[i1].x2, C(0)));
    i1++; if (Eq::not_op (e.x2, C(0))) r[j++]= e;
  }
  while (i2 < n2) {
    Pair e (v2[i2].x1, Op::op (C(0), v2[i2].x2));
    i2++; if (Eq::not_op (e.x2, C(0))) r[j++]= e;
  }
  return Sparse_vector (r, j, n);
}

template<typename Op, typename C, typename T, typename V> Sparse_vector
binary_map_optimized (const Sparse_vector& v1, const Sparse_vector& v2) {
  // With optimizations for addition and subtraction
  typedef typename V::val_op Eq;
  typedef typename V::l_op Less;
  nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
  Pair* r= mmx_new<Pair > (n);
  for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
    if (Less::op (v1[i1].x1, v2[i2].x1)) {
      typedef typename Op::lop Lop;
      r[j++]= Pair (v1[i1].x1, Lop::op (v1[i1].x2)); i1++;
    }
    else if (Less::op (v2[i2].x1, v1[i1].x1)) {
      typedef typename Op::rop Rop;
      r[j++]= Pair (v2[i2].x1, Rop::op (v2[i2].x2)); i2++;
    }
    else {
      Pair e (v2[i2].x1, Op::op (v1[i1].x2, v2[i2].x2));
      i1++; i2++; if (Eq::not_op (e.x2, C (0))) r[j++]= e;
    }
  }
  while (i1 < n1) {
    typedef typename Op::lop Lop;
    r[j++]= Pair (v1[i1].x1, Lop::op (v1[i1].x2)); i1++;
  }
  while (i2 < n2) {
    typedef typename Op::rop Rop;
    r[j++]= Pair (v2[i2].x1, Rop::op (v2[i2].x2)); i2++;
  }
  return Sparse_vector (r, j, n);
}

TMPL inline Sparse_vector
operator + (const Sparse_vector& t, const Sparse_vector& u) {
  return binary_map_optimized<add_op,C,T,V> (t, u); }
TMPL inline Sparse_vector
operator - (const Sparse_vector& t, const Sparse_vector& u) {
  return binary_map_optimized<sub_op,C,T,V> (t, u); }
TMPL inline Sparse_vector
operator * (const Sparse_vector& t, const Sparse_vector& u) {
  return binary_map<mul_op,C,T,V> (t, u); }
TMPL inline Sparse_vector
sup (const Sparse_vector& t, const Sparse_vector& u) {
  return binary_map<sup_op,C,T,V> (t, u); }
TMPL inline Sparse_vector
inf (const Sparse_vector& t, const Sparse_vector& u) {
  return binary_map<inf_op,C,T,V> (t, u); }

/******************************************************************************
* Scalar operations
******************************************************************************/

template<typename Op, typename C, typename T, typename V>
Sparse_vector
binary_map_scalar (const Sparse_vector& v, const C& x) {
  typedef typename V::val_op Eq;
  nat i, j, n= N(v);
  Pair* r= mmx_new<Pair > (n);
  for (i=0, j=0; i<n; i++) {
    Pair e (v[i].x1, Op::op (v[i].x2, x));
    if (Eq::not_op (e.x2, C (0))) r[j++]= e;
  }
  return Sparse_vector (r, j, n);
}

TMPL inline Sparse_vector
operator * (const Sparse_vector& t, const C& sc) {
  return binary_map_scalar<rmul_op> (t, sc); }
TMPL inline Sparse_vector
operator * (const C& sc, const Sparse_vector& t) {
  return binary_map_scalar<lmul_op> (t, sc); }
TMPL inline Sparse_vector
operator / (const Sparse_vector& t, const C& sc) {
  return binary_map_scalar<rdiv_op> (t, sc); }
TMPL inline Sparse_vector
operator / (const C& sc, const Sparse_vector& t) {
  return binary_map_scalar<ldiv_op> (t, sc); }

/******************************************************************************
* Composed operations
******************************************************************************/

template<typename Mul, typename C, typename T, typename V> Sparse_vector
_mul_add (const Sparse_vector& v1, const Sparse_vector& v2, const C& x) {
  // v1 + x * v2
  typedef typename V::val_op Eq;
  typedef typename V::l_op Less;
  if (is_exact_zero (x)) return v1;
  nat i1, i2, n1= N(v1), n2= N(v2), j, n= n1 + n2;
  Pair* r= mmx_new<Pair > (n);
  for (i1= i2= j= 0; i1 < n1 && i2 < n2; ) {
    if (Less::op (v1[i1].x1, v2[i2].x1)) {
      r[j++]= v1[i1]; i1++;
    }
    else if (Less::op (v2[i2].x1, v1[i1].x1)) {
      r[j++]= Pair (v2[i2].x1, Mul::op (x, v2[i2].x2)); i2++;
    }
    else {
      Pair e (v2[i2].x1, v1[i1].x2 + Mul::op (x, v2[i2].x2));
      i1++; i2++; if (Eq::not_op (e.x2, C (0))) r[j++]= e;
    }
  }
  while (i1 < n1) {
    r[j++]= Pair (v1[i1].x1, v1[i1].x2); i1++;
  }
  while (i2 < n2) {
    r[j++]= Pair (v2[i2].x1,  Mul::op (x, v2[i2].x2)); i2++;
  }
  return Sparse_vector (r, j, n);
}

template<typename C, typename T, typename V> inline Sparse_vector
sparse_mul_add (const Sparse_vector& v1, const Sparse_vector& v2, const C& x) {
  return _mul_add<mul_op> (v1, v2, x); }

#undef TMPL_DEF
#undef TMPL
#undef Sparse_vector
#undef Sparse_vector_rep
#undef Pair
} // namespace mmx
#endif // __MMX_SPARSE_VECTOR_HPP