bvector.h

Go to the documentation of this file.

//-------------------------------------------------------------
// Attila Michael Bilgic, 1/97
// last update  10/04/97 AMB
//              12/14/97 KG
// $Id: bvector.h,v 1.36.2.72 2022/11/03 17:28:10 garloff Exp $
//-------------------------------------------------------------
 
 
#ifndef TBCI_BVECTOR_H
#define TBCI_BVECTOR_H  
 
#include "tbci/basics.h"
#include "tbci/vector_sig.h"
#include "tbci/tbci_traits.h"
 
// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/bvector_gd.h"
#endif
 
NAMESPACE_TBCI
 
#ifndef TBCI_DISABLE_EXCEPT
// except.h is included by basics.h
 
class VecErr : public NumErr
{
    public:
        VecErr() 
        : NumErr("Error in Vector library")  {}
        VecErr(const char* t, const long i = 0) 
        : NumErr(t, i) {}
        VecErr(const VecErr& ve)
        : NumErr(ve) {}
        virtual ~VecErr() throw() {}
};
#endif
 
#ifdef PRAGMA_I
# pragma interface "bvector.h"
#endif
 
template <typename T> class   TMatrix;
template <typename T> class    Matrix;
template <typename T> class F_TMatrix;
template <typename T> class  F_Matrix;
template <typename T> class   TVector;
template <typename T> class  TSVector;
template <typename T> class    Vector;
 
//template <typename T> class CTensor; 
//template <typename T> class Tensor; 
 

 
template <typename T>
class BVector : public BVector_Sig<T>
{
 protected:
  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
  T* vec;
  unsigned long dim;
  bool keep;  // whether to destroy or not
 
 public:
  friend class   TMatrix<T>;
  friend class F_TMatrix<T>;
  friend class  F_Matrix<T>;
#ifdef HAVE_BCXX_INHER_BREAKS_FRIEND_BUG
  friend class   TVector<T>;
  friend class    Vector<T>;
  //friend class CTensor<T>;
#endif
 
  // constructor, destructor
  explicit BVector (const unsigned long = 0);
  BVector (const T&, const unsigned long);
  BVector (const BVector<T>&) HOT;              
#ifndef NO_POD
  // variable argument number 
  BVector(const vararg va, ...);
#endif
  inline void destroy ();
  ~BVector ();
 
  // basics
  BVector<T>& init_0 ()
        { if (UNLIKELY(dim)) TBCICLEAR (vec, T, dim); return *this; }
  BVector<T>& fill (const T&) HOT;
  // 2.97 was clever enough to detect this bug!
  //BVector<T>& fill (const Matrix<T>& m);
  friend NOINST /*inline*/ BVector<T>& FRIEND_TBCI2__ bvfillm FGD (BVector<T>&, const Matrix<T>& m);
        
  BVector<T>& resize (const BVector<T>&);
  BVector<T>& resize (const unsigned long);
  BVector<T>& resize (const T&, const unsigned long);
  BVector<T>& cheapdownsize (const unsigned long);
  // BVector<T>& free () { return this->resize (0); }
  BVector<T>& clear ()  { return this->fill((T)0); };
  BVector<T>& append (const T&);
  BVector<T>& push_back (const T& value) { append(value); return *this; }
  BVector<T>& append (const BVector<T>&);
  BVector<T>& remove (const unsigned long);
  BVector<T>& revert ();
  BVector<T>& bubble_sort ();
 
  BVector<T>& swap (BVector<T>& v);
 
  // Member access
  T& operator () (const unsigned long) HOT;
  typename tbci_traits<T>::const_refval_type
     operator () (const unsigned long) const HOT;
  // internal, no checking
  const T& getcref (const unsigned long idx) const
     { return vec[idx]; }
  typename tbci_traits<T>::const_refval_type
     get(const unsigned long idx) const HOT
     { return vec[idx]; }
  T& set(const unsigned long idx) HOT
     { return vec[idx]; }
  T& set(const T& val, unsigned long idx)
     { return this->set(idx) = val; }
 
  // valarray compatibility
  T& operator [] (const unsigned long i) { return this->operator() (i); }
  typename tbci_traits<T>::const_refval_type
     operator [] (const unsigned long i) const { return this->operator() (i); }
 
  inline unsigned long size () const HOT { return dim; }
 
  // Needed for emul et al.
  T* const & vecptr () const  { return vec; }

  typedef T* iterator;
  iterator begin () { return vec; }
  iterator   end () { return vec+dim; }
  typedef const T* const_iterator;
  const_iterator begin () const { return vec; }
  const_iterator   end () const { return vec+dim; }
 
  // You'd better know, what you are doing!
  BVector<T>& setptr (T* pointer)  { vec = pointer; return *this; }
  BVector<T>& setsize (const unsigned long size)  { dim = size; return *this; }
  // access to data pointer for e.g. fortran libs
  T* const & get_fortran_vector () const  { return vecptr (); }
  void set_fortran_vector (T* pointer)  { setptr (pointer); }
  
  // assignment operators
  BVector<T>& operator = (const T& a)  { return fill (a); }
  BVector<T>& operator = (const BVector<T>&) HOT;
  BVector<T>& alias (const BVector<T>& bv)
     { destroy (); vec = bv.vec; dim = bv.dim; keep = true; return *this; }

  BVector<T>& copy (const BVector<T>& bv);
  
  // basic comparison  
  bool operator == (const BVector<T>&) const HOT;
  bool operator != (const BVector<T>& bv) const  { return !(*this == bv); }
  bool operator <= (const BVector<T>& bv) const;
  bool operator >= (const BVector<T>& bv) const;
  bool operator <  (const BVector<T>& bv) const  { return !((*this) >= bv); }
  bool operator >  (const BVector<T>& bv) const  { return !((*this) <= bv); }
   
  // val in vector 
  bool contains (const T&, unsigned long * = 0) const;
   
  // io-streams
  friend STD__ ostream& operator << FGD (STD__ ostream&, const BVector<T>&);
  friend STD__ istream& operator >> FGD (STD__ istream&, BVector<T>&);
 
  // for Index class
  //for friend BVector<T> concat FGD (const BVector<T>&, const BVector<T>&);
  BVector<T> concat (const BVector<T>&) const;
 
  /*virtual*/ static const char* vec_info()  { return "BVector"; }
  
#ifndef HAVE_PROMOTION_BUG
# ifndef HAVE_GCC295_TMPLFRNDCLS_BUG
  // Promotion (only explicit)
  template <typename U> friend class BVector;
  template <typename U> explicit BVector (const BVector<U>& bv) 
    : dim (bv.dim), keep (bv.keep)
    { 
      if (LIKELY(dim)) vec = NEW (T, dim); else vec = 0;
      if (UNLIKELY(!vec)) dim = 0;
      for (unsigned long i = 0; i < dim; i++) vec[i] = bv.vec[i];
    }
# else
  template <typename U> explicit BVector (const BVector<U>& bv) 
    : keep (false)
    { 
      dim = bv.size(); //keep = bv.keep;
      if (LIKELY(dim)) vec = NEW (T, dim); else vec = 0;
      if (UNLIKELY(!vec)) dim = 0;
      for (unsigned long i = 0; i < dim; i++) vec[i] = bv(i);
    }
# endif
#endif  
 
#ifdef HAVE_BCXX_INHER_BREAKS_FRIEND_BUG
    /* Speed mat-vec-mul */
   friend NOINST void FRIEND_TBCI2__ do_mat_vec_mult FGD (const unsigned start, const unsigned end, \
                        TVector<T> * res, const Matrix<T> * mat, const Vector<T> * vec);
   friend NOINST void FRIEND_TBCI2__ do_mat_vec_transmult FGD (const unsigned start, const unsigned end, \
                        TVector<T> * res, const Matrix<T> * mat, const Vector<T> * vec);
#endif
 
};
 
 
//definitions
 
 
template <typename T>
INLINE BVector<T>::BVector (const unsigned long c)
  : vec ((T*)0), dim (c), keep (false)
{
        if (UNLIKELY(c)) {
                vec = NEW (T, c);
                if (UNLIKELY(!vec))
                        dim = 0;
        }
}
 
#ifndef NO_POD
template <typename T>
INLINE BVector<T>::BVector (const vararg va, ...)
  : vec ((T*)0), dim (va), keep (false)
{         
        if (LIKELY(dim))
                vec = NEW (T, dim);
        if (UNLIKELY(!vec))
                dim = 0;
        va_list vl;  
        va_start (vl, va);
#if !defined(__clang__) || !defined(CPLX)
        for (unsigned long i=0; i < dim; ++i)
                vec[i] = va_arg (vl, T);
#else
        throw VecErr("vararg not supported for cplx in clang");
#warning no vararg support with cplx numbers and clang
#endif
        va_end (vl);
}
#endif
 
template <typename T>
inline void BVector<T>::destroy ()
{
        if (LIKELY(dim)) 
                TBCIDELETE(T, vec, dim);
}
 
template <typename T>
INLINE BVector<T>::~BVector ()
{
        if (UNLIKELY(!keep))
                destroy ();
}
 
template <typename T>
INLINE BVector<T>::BVector (const BVector<T>& v) 
        : vec ((T*)0), dim(0), keep (false)
{
        if (LIKELY(v.dim)) {
                vec = NEW (T, v.dim);
                if (LIKELY(vec)) {
                        dim = v.dim;
                        TBCICOPY(vec, v.vec, T, dim);
                }
        }
}
 
 
template <typename T>
INLINE BVector<T>::BVector (const T& value, const unsigned long c)
  : vec ((T*)0), dim (c), keep (false)
{
        if (LIKELY(dim))
                vec = NEW (T, dim);
        if (LIKELY(vec)) 
                TBCIFILL (vec, value, T, dim);
        else 
                dim = 0;
}
 
template <typename T>
BVector<T>& BVector<T>::cheapdownsize (const unsigned long nd)
{
        BCHK (nd > dim, VecErr, cheapdownsize does not upsize, nd, *this);
        if (UNLIKELY(!nd))
                return this->resize (nd);
        dim = nd;
        return *this;
}
 
#ifndef LAPACK_INLINE
# define LAPACK_INLINE  
#endif
template <typename T>
LAPACK_INLINE BVector<T>& BVector<T>::resize (const unsigned long c)
{
        if (LIKELY(c == dim))
                return *this;
        if (UNLIKELY(c == 0)) {
                if (LIKELY(dim))
                        TBCIDELETE(T, vec, dim);
                dim = c; 
                return *this; 
        }
        T* tmp = vec;
        if (LIKELY(dim == 0))
                vec = NEW (T, c);
        else 
                REALLOC(vec,dim,T,c);
        BCHK (!vec, VecErr, Realloc failed in resize, 0, (vec=tmp, *this));
        dim = c; 
        return *this;
}
 
template <typename T>
BVector<T>& BVector<T>::resize (const T& value, const unsigned long c)
{
        if (UNLIKELY(dim != c)) {
                if (LIKELY(dim))
                        TBCIDELETE(T, vec, dim);
                dim = c;
                if (UNLIKELY(c))
                        vec = NEW (T, c); 
                else 
                        vec = (T*)0;
        }
        if (UNLIKELY(dim && vec))
                return this->fill (value);
        else 
                dim = 0;
        return *this;
}
 
template <typename T>
BVector<T>& BVector<T>::resize (const BVector<T>& v) 
{
        if (UNLIKELY(dim))
                TBCIDELETE (T, vec, dim);
        dim = v.dim; 
        if (UNLIKELY(v.dim)) {
                vec = NEW (T, v.dim);
                if (LIKELY(vec))
                        TBCICOPY(vec, v.vec, T, dim);
                else 
                        dim = 0;
        } else
                vec = (T*)0;
        return *this;
}
 
template <typename T>
INLINE BVector<T>& BVector<T>::fill (const T& val)
{
        TBCIFILL(vec, val, T, dim);
        return *this;
}
 
template <typename T>
INLINE BVector<T>& BVector<T>::revert ()
{
        BVector<T> buf(dim);
        for (REGISTER unsigned long t = 0; t < dim; t++) 
                buf[t] = vec[dim-t-1];
        return *this = buf;
}
 
template <typename T>
BVector<T>& BVector<T>::bubble_sort ()
{
        T temp;
        if(dim > 1) {
                for (REGISTER unsigned long i = 0; i < dim-1; i++)
                        for (REGISTER unsigned long j = dim-1; i < j; j--) {
                                if ( vec[j] < vec[j-1] ) {
                                        temp = vec[j];
                                        vec[j] = vec[j-1];
                                        vec[j-1] = temp;
                                }
                        }
        }
        return *this;
}
 
//append one element to vector
template <typename T>
INLINE BVector<T>& BVector<T>::append (const T& val)
{
        T* tmp = vec;
        REALLOC (vec,dim,T,dim+1);
        BCHK (!vec, VecErr, Realloc failed in append, 0, (vec=tmp, *this));
        vec[dim++] = val;
        return *this;
}
 
//append one vector to vector
template <typename T>
INLINE BVector<T>& BVector<T>::append (const BVector<T>& v)
{
        T* tmp = vec;
        REALLOC (vec, dim, T, dim+v.dim);
        BCHK (!vec, VecErr, Realloc failed in append, 0, (vec=tmp, *this));
        TBCICOPY(vec+dim, v.vec, T, v.dim);
        dim += v.dim;
        return *this;
}
 
//remove one element
template <typename T>
BVector<T>& BVector<T>::remove (const unsigned long i)
{
        BCHK (i>=dim, VecErr, Remove idx out of range, i, *this);
        T* tmp = NEW (T, dim-1);
        BCHK (!tmp, VecErr, Memory allocation failed in remove, (dim-1)*sizeof(T), *this);
        if (UNLIKELY(i))
                TBCICOPY (tmp, vec, T, i);
        dim -= 1;
        if (UNLIKELY(dim-i))
                TBCICOPY (&tmp[i], &vec[i+1], T, (dim-i));
        TBCIDELETE(T, vec, dim+1);
        vec = tmp;
        return *this;
}
 
template <typename T>
INLINE BVector<T>& BVector<T>::swap (BVector<T>& v)
{
        //BCHK(dim != v.dim, VecErr, swap with different size vecs, v.dim, *this);
        SWAP (v.dim, dim);
        SWAP (v.vec, vec);
        SWAP (v.keep, keep);
        return *this;
}
 
template <typename T>
inline typename tbci_traits<T>::const_refval_type
  BVector<T>::operator () (const unsigned long i) const
{ 
        EXPCHK(i>=dim, VecErr, Illegal index, i, vec[0]);
        return vec[i]; 
}
 
template <typename T>
inline T& BVector<T>::operator () (const unsigned long i)
{
        EXPCHK(i>=dim, VecErr, Illegal index, i, vec[0]);
        return vec[i]; 
}
 
template <typename T>
INLINE BVector<T>& BVector<T>::operator = (const BVector<T>& a)
{
        BCHK(dim != a.dim, VecErr, Assignment from wrong dim vector, a.dim, *this);
        TBCICOPY (vec, a.vec, T, dim);
        return *this;
}
 
template <typename T>
inline BVector<T>& BVector<T>::copy (const BVector<T>& a)
{
        //BCHK(dim != a.dim, VecErr, Assignment from wrong dim vector, a.dim, *this);
        return this->resize (a);
        //TBCICOPY(vec, a.vec, T, dim);
        //return *this;
}

template <typename T>
/*inline*/ bool BVector<T>::operator == (const BVector<T>& bv) const
{
        if (LIKELY(dim != bv.dim))
                return false;
        if (LIKELY(vec == bv.vec || dim == 0))
                return true;
        if (TBCICOMP (vec, bv.vec, T, dim)) 
                return false;
        else 
                return true;
}
 
 
template<typename T>
bool BVector<T>::operator <= (const BVector<T>& bv) const
{ 
        BCHK(dim != bv.dim, VecErr, comparison between wrong dim vectors, bv.dim, false);
        for (unsigned long i=0; i<dim; i++) 
                if (UNLIKELY(vec[i] > bv.vec[i]))
                        return false;
        return true;
}
 
 
template<typename T>
bool BVector<T>::operator >= (const BVector<T>& bv) const
{ 
        BCHK(dim != bv.dim, VecErr, comparison between wrong dim vectors, bv.dim, false);
        for (unsigned long i=0; i<dim; i++) 
                if( UNLIKELY(vec[i] < bv.vec[i]) )
                        return false;
        return true;
}
 
 
template <typename T>
inline STD__ ostream& operator << (STD__ ostream& os, const BVector<T>& v)
{
#ifdef VEC_COL
        for (unsigned long i = 0; i < v.dim; i++) 
                os << v.vec[i] << "\n";
#else
        for (unsigned long i = 0; i < v.dim; i++) 
                os << v.vec[i] << " ";
#endif
        return os.flush();
}
 
#ifdef PTR
template <typename T>
STD__ istream& operator >> (STD__ istream& in, BVector<T>& v)
{
        STD__ cerr << "BVector (" << __FILE__ << ":" << __LINE__ 
                   <<") Input on pointers not possible!" << STD__ endl;
        abort();
        return in;
}
#else
template <typename T>
STD__ istream& operator >> (STD__ istream& in, BVector<T>& v)
{
        T r ALIGN(MIN_ALIGN) = 0;
        char s = ',';
 
/*
        if( !( in >> s ) )
                return in;
        if (s != '(') 
                in.putback(s);
*/
        if( !( in >> r ) )
                return in;
        v.vec[0] = r;
        for (unsigned long i = 1; i < v.dim; i++) {
                in >> s;
                if (s != ',') in.putback(s);
                if( !( in >> r ) )
                        return in;
                v.vec[i] = r;
        } 
        return in;
}
#endif
 
template <typename T>
inline bool BVector<T>::contains (const T& val, unsigned long *ind) const
{
        for (REGISTER unsigned long i = 0; i < dim; ++i)
                if (UNLIKELY(val == vec[i])) {
                        if (ind) 
                                *ind = i;
                        return true;
                }
        return false;
}
 
 
template <typename T>
BVector<T> BVector<T>::concat (const BVector<T>& b2) const
{
        BVector<T> c (dim + b2.dim);
        TBCICOPY(c.vec, vec, T, dim);
        TBCICOPY(c.vec+dim, b2.vec, T, dim);
        return c;
}
 
 
INST(template <typename T> class BVector friend BVector<T> concat (const BVector<T>&, const BVector<T>&);)
template <typename T>
inline BVector<T> concat (const BVector<T>& b1, const BVector<T>& b2)
{ return (b1.concat (b2)); }
 

#if defined(SMP) && !defined(SMP_VECSLICE)
# define SMP_VECSLICE 262144
#endif
#if defined(SMP) && defined(__i386__) && !defined(SMP_VECSCALAR)
# define NOSMP_VECSCALAR 
# define NOSMP_VECFABS
#endif
#ifdef SMP
# define SMP_VECSLICE2 (SMP_VECSLICE/sizeof(T))
#endif
 
 
#if defined(SMP) && !defined(NO_SMP_VECVEC)
INST(template <typename T> class BVector friend void job_vv_comp (struct thr_ctrl*);)
template <typename T>
void job_vv_comp (struct thr_ctrl *tc)
{       
        tc->t_res_l = 0;
        do_vv_comp (tc->t_size, 
                    (const T*)(tc->t_par[0]), (const T*)(tc->t_par[1]),
                    tc->t_res_l);
}
 
HOTDECL(template <typename T> 
int _par_comp (const unsigned long sz, const T* v1, const T* v2))
{               
        volatile long res = 0;
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        //smp_barrier();
        
        if (LIKELY(n_thr < 2)) {
                do_vv_comp<T> (sz, v1, v2, res);
        } else {
                PREFETCH_R(v1, 3); PREFETCH_R(v2, 3);
                const unsigned first = slice_offset(1, n_thr, sz, v1);
                unsigned long st, en = first;
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, v1);
                        thread_start (t, (thr_job_t)job_vv_comp<T>, en - st,
                                      v1+st, v2+st, (void*)0);
                }
                /* The first slice is handled by the main thread */
                do_vv_comp<T> (first, v1, v2, res);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        job_output out;
                        thread_wait (t, &out);
                        res += (long) out.t_res_l;
                }
        }
        return res;
}
 
INST(template <typename T> class BVector friend void job_vv_copy (struct thr_ctrl*);)
template <typename T>
void job_vv_copy (struct thr_ctrl *tc)
{       
        _tbci_copy (tc->t_size,
                    (T*)(tc->t_par[0]), (const T*)(tc->t_par[1]));
}
 
HOTDECL(template <typename T> 
void _par_copy (const unsigned long sz, T* v1, const T* v2))
{               
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        //smp_barrier();
        
        if (LIKELY(n_thr < 2)) {
                _tbci_copy<T> (sz, v1, v2);
        } else {
                PREFETCH_W(v1, 3); PREFETCH_R(v2, 3);
                const unsigned first = slice_offset(1, n_thr, sz, v1);
                unsigned long st, en = first;
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, v1);
                        thread_start (t, (thr_job_t)job_vv_copy<T>, en - st,
                                      v1+st, v2+st, (void*)0);
                }
                /* The first slice is handled by the main thread */
                _tbci_copy (first, v1, v2);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t)
                        thread_wait (t);
        }
}
 
INST(template <typename T> class BVector friend void job_vec_fill (struct thr_ctrl*);)
template <typename T>
void job_vec_fill (struct thr_ctrl *tc)
{       
        _tbci_fill<T> (tc->t_size,
                    (T*)(tc->t_par[0]), *(const T*)(tc->t_par[1]));
}
 
 
HOTDECL(template <typename T> 
void _par_fill (const unsigned long sz, T* v1, typename tbci_traits<T>::loop_const_refval_type val))
{               
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        //smp_barrier();
        
        if (LIKELY(n_thr < 2)) {
                _tbci_fill<T> (sz, v1, val);
        } else {
                PREFETCH_W(v1, 3);
                const unsigned first = slice_offset(1, n_thr, sz, v1);
                unsigned long st, en = first;
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, v1);
                        thread_start (t, (thr_job_t)job_vec_fill<T>, en - st,
                                      v1+st, &val, (void*)0);
                }
                /* The first slice is handled by the main thread */
                _tbci_fill<T> (first, v1, val);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t)
                        thread_wait (t);
        }
}
 
INST(template <typename T> class BVector friend void job_vec_fill_fn (struct thr_ctrl*);)
template <typename T>
void job_vec_fill_fn (struct thr_ctrl *tc)
{       
        _tbci_fill_fn<T> (tc->t_size,
                    (T*)(tc->t_par[0]), 
                    *(vec_fill_fn<T>*)(tc->t_par[1]), 
                    tc->t_par[2]);
}
 
 
HOTDECL(template <typename T> 
void _par_fill_fn (const unsigned long sz, T* v1, vec_fill_fn<T> fn, void *par))
{               
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        //smp_barrier();
        
        if (LIKELY(n_thr < 2)) {
                _tbci_fill_fn<T> (sz, v1, fn, par);
        } else {
                PREFETCH_W(v1, 3);
                const unsigned first = slice_offset(1, n_thr, sz, v1);
                unsigned long st, en = first;
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, v1);
                        thread_start (t, (thr_job_t)job_vec_fill_fn<T>, en - st,
                                      v1+st, (void*)&fn, par, (void*)0);
                }
                /* The first slice is handled by the main thread */
                _tbci_fill_fn<T> (first, v1, fn, par);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t)
                        thread_wait (t);
        }
}
 
#else
template <typename T> 
int _par_comp (const unsigned long sz, const T* v1, const T* v2)
{               
        /* volatile */ long res = 0;
        do_vv_comp<T> (sz, v1, v2, res);
        return res;
}
 
template <typename T>
void _par_copy(const unsigned long sz, T* v1, const T* v2)
{
        _tbci_copy<T> (sz, v1, v2);
}
 
template <typename T>
void _par_fill (const unsigned long sz, T* v1, typename tbci_traits<T>::loop_const_refval_type val)
{
        _tbci_fill<T> (sz, v1, val);
}
 
template <typename T> 
void _par_fill_fn (const unsigned long sz, T* v1, vec_fill_fn<T> fn, void *par)
{
        _tbci_fill_fn<T> (sz, v1, fn, par);
}
#endif  
 
/* TODO: Explicit instatiation! */
#if defined(SMP) && defined(HAVE_LIBNUMA)
template <typename T>
int numa_optimize(const BVector<T>& bv, bool fault_in)
{
        if (!numa_avail && !fault_in)
                return 0;
        /* use some heuristic to decide for the num of threads */
        const unsigned long sz = bv.size();
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        //update_n_thr(n_thr);
        //smp_barrier();
        const T* v1 = &bv.getcref(0);
        
        if (LIKELY(n_thr < 2)) {
                if (fault_in)
                        do_numa_move_pages(main_numa_node, 1, 
                                (unsigned long)v1, (unsigned long)(v1+sz));
                return 0;
        } else {
                const unsigned first = slice_offset(1, n_thr, sz, v1);
                unsigned long st, en = first;
                /* Start threads */
                unsigned long res = 0;
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, v1);
                        thread_start_off (t, (thr_job_t)numa_move_pages_job, 
                                          threads[t].numa_node, fault_in,
                                          v1+st, v1+en, (void*)0);
                }
                /* The first slice is handled by the main thread */
                res = do_numa_move_pages(main_numa_node, fault_in,
                                (unsigned long)v1, (unsigned long)(v1+first));
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        job_output out; 
                        thread_wait (t, &out);
                        res += out.t_res_l;
                }
                //fprintf(stderr, "NUMA Optimize BVector %p: %li pages moved\n", v1, res); 
                return res;
        }
}
#else
template <typename T>
int numa_optimize(const BVector<T>& bv, bool fault_in)
{
        /* FIXME: We should emulate fault_in here, no? */
        return 0;
}
#endif
 
NAMESPACE_END
 
#endif  /* TBCI_BVECTOR_H */