vector.h

Go to the documentation of this file.

#ifndef TBCI_VECTOR_H
#define TBCI_VECTOR_H  

#include "tbci/bvector.h"
//#include "cplx.h"

// Include vector kernels
#include "tbci/vec_kern_unr_pref.h"

#ifndef TBCI_NO_SIMD
# if defined(__AVX__) && !defined(NO_AVX)
#  include "tbci/vec_kern_special2.h"
# else
#  include "tbci/vec_kern_special.h"
# endif
#endif

// Avoid -fguiding-decls
#if !defined(NO_GD) && !defined(AUTO_DECL)
# include "tbci/vector_gd.h"
#endif

NAMESPACE_TBCI

//template <typename T> class cplx;

template <typename T> class Vector;
template <typename T> class TVector;
template <typename T> class TSVector;
template <typename T> class Matrix;
template <typename T> class TMatrix;
template <typename T> class TSMatrix;
template <typename T> class F_Matrix;
template <typename T> class F_TMatrix;
template <typename T> class F_TSMatrix;
template <typename T> class Mat_Brack;
template <typename T> class CRMatrix;
template <typename T> class CSCMatrix;
template <typename T> class BdMatrix;
template <typename T> class Symm_BdMatrix;

#ifdef PRAGMA_I
# pragma interface "vector.h"
#endif

//------------------------------------------------------------------------

template <typename T>
class TVector : public BVector<T>, public Vector_Sig<T>
{
 public:
  friend class Vector<T>;
  friend class TSVector<T>;  // shouldn't be necessary ...
  friend class Matrix<T>;
  friend class BdMatrix<T>;

  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);

  explicit TVector (const unsigned long d = 0)  : BVector<T> (d) {}
  TVector (const T& val, const unsigned long d) : BVector<T> (val, d)  {}
  TVector (const BVector<T>& bv) : BVector<T> (bv)  {}
  TVector (const  Vector<T>&  v) : BVector<T> ( v)  {}
  //alias
  TVector (const  TVector<T>& tv) HOT
  { this->dim = tv.dim; this->vec = tv.vec; }
  explicit TVector (const TSVector<T>& ts) 
  { ts.eval (); this->vec = ts.vec; this->dim = ts.dim; }
  TVector (const Mat_Brack<T>&);
  // dtor: Don't destroy!
  ~TVector ()  { this->keep = true; }
  /* TODO: Check whether it would be more efficient to save time by
   * setting keep in c'tors or even only when aliasing takes place;
   * we'd be doing refcounting more explicitly this way. 
   * Question is whether we'd not end up with full ref counting, which
   * is what TBCI tries to avoid.       -KG, 2004-04-27.
   */

  // avoid ambiguity
  unsigned long size () const  { return this->dim; }

  TVector<T>& operator  = (const T& a) 
    { this->BVector<T>::operator = (a); return *this; }
  TVector<T>& operator  = (const Vector<T>& v) 
    { this->BVector<T>::operator = (v); return *this; }
  TVector<T>& operator  = (const  TVector<T>& tv);
  TVector<T>& operator  = (const TSVector<T>& ts);

  bool operator == (const TVector<T>& tv) const
    { Vector<T> v(*this); if (LIKELY(this->vec == tv.vec)) return true;
      /* else { */ Vector<T> v2(tv); return (v == v2); /*}*/ }
      
  bool operator != (const TVector<T>& tv) const  { return !(*this == tv); }
  bool operator == (const Vector<T>& v) const
    { Vector<T> tv(*this); return (tv == v); }
  bool operator != (const Vector<T>& v) const  { return !(*this == v); }
  bool operator == (const BVector<T>& bv) const
    { Vector<T> v(*this); return bv.operator == (v); }
  bool operator != (const BVector<T>& v) const  { return !(*this == v); }
  bool operator == (const TSVector<T>& tsv) const  { return (tsv == *this); }
  bool operator != (const TSVector<T>& tsv) const  { return !(*this == tsv); }

  TVector<T>& operator += (const T&);
  TVector<T>& operator -= (const T&);
  TVector<T>& operator *= (const T&);
  TVector<T>& operator /= (const T&);

  TVector<T>& operator += (const Vector<T>&) HOT;
  TVector<T>& operator -= (const Vector<T>&) HOT;

  TVector<T>& operator += (const TVector<T>& tv)
    { return this->operator += (Vector<T> (tv)); }
  TVector<T>& operator -= (const TVector<T>& tv)
    { return this->operator -= (Vector<T> (tv)); }

  TVector<T>& operator += (const TSVector<T>& tsv) HOT;
  TVector<T>& operator -= (const TSVector<T>& tsv) HOT;

  // addtinional arithmetics (valarray like)

  T min ()  { Vector<T> th (*this); return th.min (); }
  T max ()  { Vector<T> th (*this); return th.max (); }
  T sum ()  { Vector<T> th (*this); return th.sum (); }

  double fabssqr ()  { Vector<T> th (*this); return th.fabssqr (); }
  double fabs ()     { Vector<T> th (*this); return th.fabs ();    }
  T abs ()           { Vector<T> th (*this); return th.abs ();     }
  /* NOTE: The const is wrong! Just to make the compiler happy, when compiling 
   * operator / (const T&, const TVector<T>&) */
  
  // basic operations
  
  TVector<T>& operator + (const Vector<T>&) HOT;
  TVector<T>& operator - (const Vector<T>&) HOT;
  TVector<T>& operator + (const TSVector<T>& ts) HOT;
  TVector<T>& operator - (const TSVector<T>& ts) HOT;
  const TVector<T>& operator + (const TVector<T>& a) /*const*/ HOT;
  const TVector<T>& operator - (const TVector<T>& a) /*const*/ HOT;

  // fortunately, we don't need friendship here
#ifdef TV_OP_FRIENDS
  friend NOINST const TVector<T>& operator + FGD (const T&, const TVector<T>&);
  friend NOINST const TVector<T>& operator - FGD (const T&, const TVector<T>&);
  friend NOINST TSVector<T> operator * FGD (const T&, const TVector<T>&);
  friend NOINST TSVector<T> operator / FGD (const T&, const TVector<T>&);
  friend NOINST TVector<T>  operator + FGD (const T&, const  Vector<T>&);
  friend NOINST TVector<T>  operator - FGD (const T&, const  Vector<T>&);
  friend NOINST TSVector<T> operator * FGD (const T&, const  Vector<T>&);
  //friend NOINST TVector<T>  operator / FGD (const T&, const Vector<T>&);

  // More friendships ...
  friend NOINST TVector<T> operator + FGD (const T&, const TSVector<T>&);
  friend NOINST TVector<T> operator - FGD (const T&, const TSVector<T>&);
#endif

  // operator () destroys TVector !!! 
  T operator () (const unsigned long i)
    { T val = this->vec[i]; /*if (!keep)*/ this->destroy (); return val; }
  //T& operator () (const unsigned long i) { STD__ cerr << "TVectors may not be written to by the () op\n"; abort(); }
  T operator [] (const unsigned long i)  { return this->operator() (i); }
  //T& operator [] (const unsigned long i) { return this->operator() (i); }
  
  const T& getcref (const unsigned long i) const  { return this->vec[i]; }
  typename tbci_traits<T>::const_refval_type
     get (const unsigned long i) const   { return this->vec[i]; }

  T& setval (const unsigned long i) const { return this->vec[i]; }
  T& setval (const T& val, const unsigned long i) const { return this->vec[i] = val; }
  T& set (const T& val, const unsigned long i) const    { return this->vec[i] = val; }

  TVector<T> slice (const unsigned long, const unsigned long);
  TVector<T>& incr (const unsigned long, const T = (T)1);
 
  TVector<T>& operator + (const T&);
  TVector<T>& operator - (const T&);
  TSVector<T> operator * (const T&);
  TSVector<T> operator / (const T&);
  TSVector<T> operator - ();

  T operator * (const Vector<T>& v)
    { Vector<T> th (*this); return th * v; }
  T operator * (TVector<T>& tv)
    { Vector<T> th (*this); Vector<T> v (tv); return th * v; }


  friend TVector<T>  conj FGD (const Vector<T>&);
  friend TVector<T>  real FGD (const Vector<T>&);
  friend TVector<T>  imag FGD (const Vector<T>&);

  friend TVector<T>& conj FGD (TVector<T>&);
  friend TVector<T>& real FGD (TVector<T>&);
  friend TVector<T>& imag FGD (TVector<T>&);
  
  static const char* vec_info()  { return "TVector"; }
   
  friend STD__ ostream& operator << FGD (STD__ ostream&, const TVector<T>&);

#ifndef HAVE_PROMOTION_BUG
  // Promotion (only explicit)
//# ifndef HAVE_GCC295_TMPLFRNDCLS_BUG
//  template <typename U> friend class BVector;
//# endif
  template <typename U> explicit TVector (const BVector<U>& bv) : BVector<T> (bv) {}
#endif  
   /* 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_tsv_mult FGD (const unsigned start, const unsigned end, \
                        TVector<T> * res, const Matrix<T> * mat, const TSVector<T> * rsv);
   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);
} ;

/* ************************ definitions ********************************** */

template <typename T>
inline TVector<T>& TVector<T>::operator = (const TVector<T>& tv)
{ 
        BCHK (this->dim != tv.dim, VecErr, diff size in assignment, tv.dim, *this);
        if (LIKELY(this->vec != tv.vec)) {
                this->destroy (); this->vec = tv.vec; this->dim = tv.dim; 
        }
        return *this; 
}


template <typename T>
inline STD__ ostream& operator << (STD__ ostream& os, const TVector<T>& tv)
{ Vector<T> v(tv); return os << (BVector<T>)v; }

template <typename T>
inline TVector<T>& TVector<T>::operator = (const TSVector<T>& tsv)
{ 
        BCHK (this->dim != tsv.dim, VecErr, diff size in assignment, tsv.dim, *this);
        if (LIKELY(tsv.mut && tsv.fac == (T)1)) {
                this->destroy (); this->vec = tsv.vec; 
        }
        else 
                tsv.eval (this->vec);
        return *this; 
}

/*
template <typename T>
inline TVector<T>::TVector (const TSVector<T>& ts)
{
   ts.eval (); vec = ts.vec; dim = ts.dim; ts.mut = false;
}
 */

#if !defined(SMP) || defined(NOSMP_VECVEC)
/* nothing, actually */
#else   /* SMP */

INST(template <typename T> class Vector friend void job_vec_vec_add (struct thr_ctrl *);)
template <typename T>
void job_vec_vec_add (struct thr_ctrl *tc)
{
        do_vec_vec_add (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_vec_sub (struct thr_ctrl *);)
template <typename T>
void job_vec_vec_sub (struct thr_ctrl *tc)
{
        do_vec_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_add_vec (struct thr_ctrl *);)
template <typename T>
void job_vec_add_vec (struct thr_ctrl *tc)
{
        do_vec_add_vec (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_sub_vec (struct thr_ctrl *);)
template <typename T>
void job_vec_sub_vec (struct thr_ctrl *tc)
{
        do_vec_sub_vec (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_sub_vec_inv (struct thr_ctrl *);)
template <typename T>
void job_vec_sub_vec_inv (struct thr_ctrl *tc)
{
        do_vec_sub_vec_inv (tc->t_size, (T*)(tc->t_par[0]),
                            (const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_val_add (struct thr_ctrl *);)
template <typename T>
void job_vec_val_add (struct thr_ctrl *tc)
{
        do_vec_val_add (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_val_sub (struct thr_ctrl *);)
template <typename T>
void job_vec_val_sub (struct thr_ctrl *tc)
{
        do_vec_val_sub (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_val_mul (struct thr_ctrl *);)
template <typename T>
void job_vec_val_mul (struct thr_ctrl *tc)
{
        do_vec_val_mul (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_val_vec_mul (struct thr_ctrl *);)
template <typename T>
void job_val_vec_mul (struct thr_ctrl *tc)
{
        do_val_vec_mul (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_val_vec_add (struct thr_ctrl *);)
template <typename T>
void job_val_vec_add (struct thr_ctrl *tc)
{
        do_val_vec_add (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_val_vec_sub (struct thr_ctrl *);)
template <typename T>
void job_val_vec_sub (struct thr_ctrl *tc)
{
        do_val_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_val_vec_div (struct thr_ctrl *);)
template <typename T>
void job_val_vec_div (struct thr_ctrl *tc)
{
        do_val_vec_div (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_add_val (struct thr_ctrl *);)
template <typename T>
void job_vec_add_val (struct thr_ctrl *tc)
{
        do_vec_add_val (tc->t_size, (T*)(tc->t_par[0]),
                        *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_val_add_vec (struct thr_ctrl *);)
template <typename T>
void job_val_add_vec (struct thr_ctrl *tc)
{
        do_val_add_vec(tc->t_size, (T*)(tc->t_par[0]),
                       *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_sub_val (struct thr_ctrl *);)
template <typename T>
void job_vec_sub_val (struct thr_ctrl *tc)
{
        do_vec_sub_val (tc->t_size, (T*)(tc->t_par[0]),
                        *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_val_sub_vec (struct thr_ctrl *);)
template <typename T>
void job_val_sub_vec (struct thr_ctrl *tc)
{
        do_val_sub_vec(tc->t_size, (T*)(tc->t_par[0]),
                       *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_mul_val (struct thr_ctrl *);)
template <typename T>
void job_vec_mul_val (struct thr_ctrl *tc)
{
        do_vec_mul_val (tc->t_size, (T*)(tc->t_par[0]),
                        *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_vec_div_val (struct thr_ctrl *);)
template <typename T>
void job_vec_div_val (struct thr_ctrl *tc)
{
        do_vec_div_val (tc->t_size, (T*)(tc->t_par[0]),
                        *(const T*)(tc->t_par[1]));
}
INST(template <typename T> class Vector friend void job_val_div_vec (struct thr_ctrl *);)
template <typename T>
void job_val_div_vec (struct thr_ctrl *tc)
{
        do_val_div_vec(tc->t_size, (T*)(tc->t_par[0]),
                       *(const T*)(tc->t_par[1]));
}

/* TSV variants */
INST(template <typename T> class Vector friend void job_vec_svc_add (struct thr_ctrl *);)
template <typename T>
void job_vec_svc_add (struct thr_ctrl *tc)
{
        do_vec_svc_add (tc->t_size, (T*)(tc->t_par[0]),
                        (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_svc_vec_add (struct thr_ctrl *);)
template <typename T>
void job_svc_vec_add (struct thr_ctrl *tc)
{
        do_svc_vec_add (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_svc_svc_add (struct thr_ctrl *);)
template <typename T>
void job_svc_svc_add (struct thr_ctrl *tc)
{
        do_svc_svc_add (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),  (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]), *(const T*)(tc->t_par[4]));
}
INST(template <typename T> class Vector friend void job_vec_svc_sub (struct thr_ctrl *);)
template <typename T>
void job_vec_svc_sub (struct thr_ctrl *tc)
{
        do_vec_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_svc_vec_sub (struct thr_ctrl *);)
template <typename T>
void job_svc_vec_sub (struct thr_ctrl *tc)
{
        do_svc_vec_sub (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]), (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_svc_svc_sub (struct thr_ctrl *);)
template <typename T>
void job_svc_svc_sub (struct thr_ctrl *tc)
{
        do_svc_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),  (const T*)(tc->t_par[2]),
                        *(const T*)(tc->t_par[3]), *(const T*)(tc->t_par[4]));
}
INST(template <typename T> class Vector friend void job_vec_add_svc (struct thr_ctrl *);)
template <typename T>
void job_vec_add_svc (struct thr_ctrl *tc)
{
        do_vec_add_svc (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_sub_svc (struct thr_ctrl *);)
template <typename T>
void job_vec_sub_svc (struct thr_ctrl *tc)
{
        do_vec_sub_svc (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]));
}
INST(template <typename T> class Vector friend void job_vec_sub_svc_inv (struct thr_ctrl *);)
template <typename T>
void job_vec_sub_svc_inv (struct thr_ctrl *tc)
{
        do_vec_sub_svc_inv (tc->t_size, (T*)(tc->t_par[0]),
                             (const T*)(tc->t_par[1]),
                            *(const T*)(tc->t_par[2]));
}

INST(template <typename T> class Vector friend void job_svc_val_add (struct thr_ctrl *);)
template <typename T>
void job_svc_val_add (struct thr_ctrl *tc)
{
        do_svc_val_add (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3])); 
}
INST(template <typename T> class Vector friend void job_svc_val_sub (struct thr_ctrl *);)
template <typename T>
void job_svc_val_sub (struct thr_ctrl *tc)
{
        do_svc_val_sub (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3])); 
}
INST(template <typename T> class Vector friend void job_val_svc_add (struct thr_ctrl *);)
template <typename T>
void job_val_svc_add (struct thr_ctrl *tc)
{
        do_val_svc_add (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_val_svc_sub (struct thr_ctrl *);)
template <typename T>
void job_val_svc_sub (struct thr_ctrl *tc)
{
        do_val_svc_sub (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
}
INST(template <typename T> class Vector friend void job_val_svc_div (struct thr_ctrl *);)
template <typename T>
void job_val_svc_div (struct thr_ctrl *tc)
{
        do_val_svc_div (tc->t_size, (T*)(tc->t_par[0]),
                         (const T*)(tc->t_par[1]),
                        *(const T*)(tc->t_par[2]), *(const T*)(tc->t_par[3]));
}


/* TODO: conj, unary - */


#endif

/* Naming rules:
 *  1. first arg is always pointer
 *  2. we've max 3 pointers and two consts
 *  3. the last two args are designated by letters, C = const, V = vector
 *  4. letters that are clear because of rules 2 or 3 are omitted
 *  5. pointer(vector) args are always before consts
 */
#if !defined(SMP) || defined(NOSMP_VECVEC)
#define STD_SMP_TEMPLATE2V(oper, dm, a1, a2)                            \
        do_##oper <T> (dm, a1, a2)
#define STD_SMP_TEMPLATE2C(oper, dm, a1, a2)                            \
        do_##oper <T> (dm, a1, a2)
#define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3)                       \
        do_##oper <T> (dm, a1, a2, a3)
#define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3)                       \
        do_##oper <T> (dm, a1, a2, a3)
#define STD_SMP_TEMPLATE3CC(oper, dm, a1, a2, a3)                       \
        do_##oper <T> (dm, a1, a2, a3)
#define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4)                    \
        do_##oper <T> (dm, a1, a2, a3, a4)
#define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4)                    \
        do_##oper <T> (dm, a1, a2, a3, a4)
#define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5)                 \
        do_##oper <T> (dm, a1, a2, a3, a4, a5)

#else   /* SMP */

#define _SMP_TMPL2V(oper,st,en,a1,a2,a3,a4,a5)   do_##oper <T> (en-st, a1+st, a2+st)
#define _JOB_TMPL2V(oper,st,en,a1,a2,a3,a4,a5)  job_##oper <T>, en-st, a1+st, a2+st
#define _SMP_TMPL2C(oper,st,en,a1,a2,a3,a4,a5)   do_##oper <T> (en-st, a1+st,  a2)
#define _JOB_TMPL2C(oper,st,en,a1,a2,a3,a4,a5)  job_##oper <T>, en-st, a1+st, &a2
#define _SMP_TMPL3VV(oper,st,en,a1,a2,a3,a4,a5)  do_##oper <T> (en-st, a1+st, a2+st, a3+st)
#define _JOB_TMPL3VV(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, a3+st
#define _SMP_TMPL3VC(oper,st,en,a1,a2,a3,a4,a5)  do_##oper <T> (en-st, a1+st, a2+st,  a3)
#define _JOB_TMPL3VC(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, a2+st, &a3
#define _SMP_TMPL3CC(oper,st,en,a1,a2,a3,a4,a5)  do_##oper <T> (en-st, a1+st,  a2,  a3)
#define _JOB_TMPL3CC(oper,st,en,a1,a2,a3,a4,a5) job_##oper <T>, en-st, a1+st, &a2, &a3
#define _SMP_TMPL4V(oper,st,en,a1,a2,a3,a4,a5)   do_##oper <T> (en-st, a1+st, a2+st, a3+st,  a4)
#define _JOB_TMPL4V(oper,st,en,a1,a2,a3,a4,a5)  job_##oper <T>, en-st, a1+st, a2+st, a3+st, &a4
#define _SMP_TMPL4C(oper,st,en,a1,a2,a3,a4,a5)   do_##oper <T> (en-st, a1+st, a2+st,  a3,  a4)
#define _JOB_TMPL4C(oper,st,en,a1,a2,a3,a4,a5)  job_##oper <T>, en-st, a1+st, a2+st, &a3, &a4
#define _SMP_TMPL5(oper,st,en,a1,a2,a3,a4,a5)    do_##oper <T> (en-st, a1+st, a2+st, a3+st,  a4,  a5)
#define _JOB_TMPL5(oper,st,en,a1,a2,a3,a4,a5)   job_##oper <T>, en-st, a1+st, a2+st, a3+st, &a4, &a5

#define STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, a5, SMP_TMPL, JOB_TMPL, p2, p3)      \
        /* use some heuristic to decide for the num of threads */               \
        const unsigned n_thr = threads_avail (dm / SMP_VECSLICE2);              \
        update_n_thr(n_thr);                                                    \
        if (LIKELY(n_thr < 2)) {                                                \
                SMP_TMPL(oper, 0, dm, a1, a2, a3, a4, a5);                      \
        } else {                                                                \
                PREFETCH_W_MANY(a1, 3);                                         \
                p2                                                              \
                p3                                                              \
                const unsigned long first = slice_offset(2, n_thr, dm, a1);     \
                unsigned long st, en = first;                                   \
                unsigned _t;                                                    \
                /* smp_barrier(); */                                            \
                /* Start threads */                                             \
                for (_t = 0; _t < n_thr-1; ++_t) {                              \
                        st = en; en = slice_offset(_t+2, n_thr, dm, a1);        \
                        thread_start (_t, JOB_TMPL(oper, st, en, a1, a2, a3, a4, a5), (void*)0); \
                }                                                               \
                /* The first slice is handled by the main thread */             \
                SMP_TMPL(oper, 0UL, first, a1, a2, a3, a4, a5);                 \
                /* sched_yield (); */                                           \
                /* Wait for the end */                                          \
                for (_t = 0; _t < n_thr-1; ++_t)                                \
                        thread_wait (_t);                                       \
        }

#define STD_SMP_TEMPLATE2V(oper, dm, a1, a2)                            \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, NULL, NULL, NULL, _SMP_TMPL2V, _JOB_TMPL2V, PREFETCH_R_MANY(a2,2);,)
#define STD_SMP_TEMPLATE2C(oper, dm, a1, a2)                            \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, NULL, NULL, NULL, _SMP_TMPL2C, _JOB_TMPL2C,,)
#define STD_SMP_TEMPLATE3VV(oper, dm, a1, a2, a3)                       \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3VV, _JOB_TMPL3VV, PREFETCH_R_MANY(a2,2);, PREFETCH_R_MANY(a3,2);)
#define STD_SMP_TEMPLATE3VC(oper, dm, a1, a2, a3)                       \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3VC, _JOB_TMPL3VC, PREFETCH_R_MANY(a2,2);,)
#define STD_SMP_TEMPLATE3CC(oper, dm, a1, a2, a3)                       \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, NULL, NULL, _SMP_TMPL3CC, _JOB_TMPL3CC,,)
#define STD_SMP_TEMPLATE4V(oper, dm, a1, a2, a3, a4)                    \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, NULL, _SMP_TMPL4V, _JOB_TMPL4V, PREFETCH_R_MANY(a2,2);,)
#define STD_SMP_TEMPLATE4C(oper, dm, a1, a2, a3, a4)                    \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, NULL, _SMP_TMPL4C, _JOB_TMPL4C, PREFETCH_R_MANY(a2,2);,)
#define STD_SMP_TEMPLATE5(oper, dm, a1, a2, a3, a4, a5)                 \
        STD_SMP_TEMPLATE(oper, dm, a1, a2, a3, a4, a5, _SMP_TMPL5, _JOB_TMPL5, PREFETCH_R_MANY(a2,2);, PREFETCH_R_MANY(a3,2);)


#endif  /* SMP */

template <typename T>
inline TVector<T>& TVector<T>::operator += (const T& a)
{
        STD_SMP_TEMPLATE2C(vec_add_val, this->dim, this->vec, a);
        return *this;
}
template <typename T>
inline TVector<T>& TVector<T>::operator -= (const T& a)
{
        STD_SMP_TEMPLATE2C(vec_sub_val, this->dim, this->vec, a);
        return *this;
}
template <typename T>
inline TVector<T>& TVector<T>::operator *= (const T& a)
{
        STD_SMP_TEMPLATE2C(vec_mul_val, this->dim, this->vec, a);
        return *this;
}

template <typename T>
inline TVector<T>& TVector<T>::operator /= (const T& a)
{       
        BCHK (a == (T)0, VecErr, Divide by zero (TVector), 0, *this);
        return this->operator *= ((T)1/a);
}


template <typename T>
inline TVector<T>& TVector<T>::operator += (const Vector<T>& a)
{
        BCHK(this->dim!=a.dim,VecErr,"TV += V operator with diff. sizes", a.dim, *this);
        STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, this->vec, a.vec);
        return *this;
}
template <typename T>
inline TVector<T>& TVector<T>::operator -= (const Vector<T>& a)
{               
        BCHK(this->dim!=a.dim,VecErr,"TV -= V operator with diff. sizes", a.dim, *this);
        STD_SMP_TEMPLATE2V(vec_sub_vec, this->dim, this->vec, a.vec);
        return *this;
}

template <typename T>
inline TVector<T>& TVector<T>::operator += (const TSVector<T>& a)
{ 
        BCHK(this->dim!=a.dim,VecErr,"TV += TSV operator with diff. sizes", a.dim, *this);
        STD_SMP_TEMPLATE3VC(vec_add_svc, this->dim, this->vec, a.vec, a.fac);
        if (LIKELY(this->vec != a.vec))
                a.destroy ();
        return *this;
}
template <typename T>
inline TVector<T>& TVector<T>::operator -= (const TSVector<T>& a)
{ 
        BCHK(this->dim!=a.dim,VecErr,"TV -= TSV operator with diff. sizes", a.dim, *this);
        STD_SMP_TEMPLATE3VC(vec_sub_svc, this->dim, this->vec, a.vec, a.fac);
        if (LIKELY(this->vec != a.vec))
                a.destroy ();
        return *this;
}


// Binary operations (for TV id. to unary)

template <typename T>
inline TVector<T>& TVector<T>::operator + (const Vector<T>& a)
{
        //BCHK(dim != a.dim, VecErr, Different sizes in operator TV + V, a.dim, *this);
        return this->operator += (a);
}
template <typename T>
inline TVector<T>& TVector<T>::operator - (const Vector<T>& a)
{
        //BCHK(dim != a.dim, VecErr, Different sizes in operator TV - V, a.dim, *this);
        return this->operator -= (a);
}


template <typename T>
const TVector<T>& TVector<T>::operator + (const TVector<T>& tv) /*const*/
{
        BCHK (tv.dim != this->dim, VecErr, Diff. dim in TV + TV, tv.dim, (this->destroy(), tv));
        //for (REGISTER unsigned long i = 0; i < dim; ++i) tv.vec[i] = vec[i] + tv.vec[i];
        STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, tv.vec, this->vec);
        this->destroy (); return tv;
}

template <typename T>
const TVector<T>& TVector<T>::operator - (const TVector<T>& tv) /*const*/
{ 
        BCHK (tv.dim != this->dim, VecErr, Diff. dim in TV - TV, tv.dim, (this->destroy(), tv));
        //for (REGISTER unsigned long i = 0; i < dim; ++i) tv.vec[i] = vec[i] - tv.vec[i];
        STD_SMP_TEMPLATE2V(vec_sub_vec_inv, this->dim, tv.vec, this->vec);
        this->destroy (); return tv;
}


INST(template<typename T> class TVector friend const TVector<T>& operator + (const T&, const TVector<T>&);)
template <typename T>
inline const TVector<T>& operator + (const T& a, const TVector<T>& b)
{
        //for (REGISTER unsigned long i=0; i<b.size(); ++i) b.set(a + b.get(i), i);
        const unsigned long dim = b.size();
        T* vec = &b.setval(0);
        STD_SMP_TEMPLATE2C(val_add_vec, dim, vec, a);
        return b;
}
INST(template<typename T> class TVector friend const TVector<T>& operator - (const T&, const TVector<T>&);)
template <typename T>
inline const TVector<T>& operator - (const T& a, const TVector<T>& b)
{
        //for (REGISTER unsigned long i=0; i<b.size(); ++i) b.set(a - b.get(i), i);
        const unsigned long dim = b.size();
        T* vec = &b.setval(0);
        STD_SMP_TEMPLATE2C(val_sub_vec, dim, vec, a);
        return b;
}


INST(template <typename T> class TVector friend TSVector<T> operator * (const T&, const TVector<T>&);)
template <typename T>
inline TSVector<T> operator * (const T& a, const TVector<T>& b)
{
        return TSVector<T> (b, a);
}

// a/b = a*b/b^2
// We have to take care T not to be an int type!
INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
INST(template <typename T> class TVector friend TSVector<T> operator / (const T&, TVector<T>);)
INSTCTL(#endif)
template <typename T>
inline TSVector<T> operator / (const T& a, TVector<T> b)
{
        T dv = a / b.fabssqr();
        BCHK(dv==(T)0, VecErr, Divide by norm 0 Vector, 0, 0);
        return TSVector<T> (b, dv);
}


template <typename T>
inline TVector<T>& TVector<T>::operator + (const T& a)
{ return this->operator += (a); }
template <typename T>
inline TVector<T>& TVector<T>::operator - (const T& a)
{ return this->operator -= (a); }


template <typename T>
inline TSVector<T> TVector<T>::operator * (const T& a)
{
        return TSVector<T> (*this, a);
}

// We have to take care T not to be an int type!
template <typename T>
inline TSVector<T> TVector<T>::operator / (const T& a)
{
        return TSVector<T> (*this, (T)1/a);
}

#if 0
template <typename T>
inline TVector<T>& TVector<T>::operator - ()
{
        for (REGISTER unsigned long i=0; i<dim; ++i) 
                this->vec[i] = -this->vec[i];
        return *this;
}
#else
template <typename T>
inline TSVector<T> TVector<T>::operator - ()
{
        return TSVector<T> (*this, (T)-1);
}
#endif

//NAMESPACE_END
//NAMESPACE_CPLX

//INST(template <typename T> class TVector friend TVector<T>& conj (TVector<T>&);)
template <typename T>
inline TVector<T>& conj (TVector<T>& tv) 
{ 
        for (REGISTER unsigned long i=0; i<tv.size(); ++i) 
                tv.set (CPLX__ conj (tv.get(i)), i);
        return tv;
}

//INST(template <typename T> class TVector friend TVector<T>& real (TVector<T>&);)
template <typename T>
inline TVector<T>& real (TVector<T>& tv)
{
        for (REGISTER unsigned long i=0; i<tv.size(); ++i)
                tv.set(CPLX__ real (tv.get(i)), i);
        //tv.vec[i] = (tv.vec[i] + CPLX__ conj (tv.vec[i])) * 0.5;
        return tv;
}

//INST(template <typename T> class TVector friend TVector<T>& imag (TVector<T>&);)
template <typename T>
inline TVector<T>& imag (TVector<T>& tv)
{
        for (REGISTER unsigned long i=0; i<tv.size(); ++i)
                tv.set(CPLX__ imag (tv.get(i)), i);
        //tv.vec[i] = (tv.vec[i] - CPLX__ conj (tv.vec[i])) * 0.5;
        return tv;
}

//INST(template <typename T> class TVector friend TVector<T> conj (const Vector<T>&);)
template <typename T>
inline TVector<T> conj (const Vector<T>& v)
{ 
        TVector<T> tv(v.size());
        for (REGISTER unsigned long i=0; i<v.size(); ++i) 
                tv.set (CPLX__ conj (v(i)), i);
        return tv;
}

//INST(template <typename T> class TVector friend TVector<T> real (const Vector<T>&);)
template <typename T>
inline TVector<T> real (const Vector<T>& v)
{
        TVector<T> tv (v.size());
        for (REGISTER unsigned long i=0; i<v.size(); ++i) 
                tv.set (CPLX__ real (v(i)), i);
        //tv.vec[i] = (v.vec[i] + CPLX__ conj (v.vec[i])) * 0.5;
        return tv;
}

//INST(template <typename T> class TVector friend TVector<T> imag (const Vector<T>&);)
template <typename T>
inline TVector<T> imag (const Vector<T>& v)
{
        TVector<T> tv (v.size());
        for (REGISTER unsigned long i=0; i<v.size(); ++i)
                tv.set(CPLX__ imag (v(i)), i);
        //tv.vec[i] = (v.vec[i] - CPLX__ conj (v.vec[i])) * 0.5;
        return tv;
}

//NAMESPACE_CPLX_END
//NAMESPACE_TBCI        

template <typename T>
inline TVector<T>& TVector<T>::operator + (const TSVector<T>& ts)
{
        BCHK(this->dim!=ts.dim,VecErr,Wrong dim TV + TSV,ts.dim,*this);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  vec[t] op##= ts.fac * ts.vec[t];
        STD_SMP_TEMPLATE3VC(vec_add_svc, this->dim, this->vec, ts.vec, ts.fac);
        /* included by AA */
        if (LIKELY(this->vec != ts.vec))
                ts.destroy();
        return *this;
}

template <typename T>
inline TVector<T>& TVector<T>::operator - (const TSVector<T>& ts)
{
        BCHK(this->dim!=ts.dim,VecErr,Wrong dim TV - TSV,ts.dim,*this);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  vec[t] op##= ts.fac * ts.vec[t];
        STD_SMP_TEMPLATE3VC(vec_sub_svc, this->dim, this->vec, ts.vec, ts.fac);
        /* included by AA */
        if (LIKELY(this->vec != ts.vec))
                ts.destroy();
        return *this;
}


template <typename T>
inline TVector<T> TVector<T>::slice (const unsigned long i0, const unsigned long i1)
{
        BCHK (i1 < i0, VecErr, Slicing needs ordered arguments, i0, *this);
        BCHK (i0 >= this->dim, VecErr, Slice: arg1 out of range, i0, *this);
        BCHK (i1 > this->dim, VecErr, Slice: arg2 out of range, i1, *this);
        this->dim = i1 - i0;
        // BEWARE: This relies on the C library freeing the right chunk of
        //         memory, even if it moved within the allocated region!!!
        if (UNLIKELY(this->dim == 0)) 
                TBCIDELETE (T, this->vec, this->dim);
        else
                this->vec += i0;
        return *this;
}


template <typename T>
inline TVector<T>& TVector<T>::incr (const unsigned long wh, const T i)
{
        BCHK (wh >= this->dim, VecErr, Tried to change not-exitsting idx, wh, *this);
        this->vec[wh] += i;
        return *this;
}

//fwd-decl
//template <typename T> inline double fabssqr(TVector<T> tv);
template <typename T> inline double fabssqr(const Vector<T>& v);

NAMESPACE_END

NAMESPACE_CSTD

INST(template <typename T> class TBCI__ Vector friend double fabs (const TBCI__ Vector<T>&);)
template <typename T>
inline double fabs (const TBCI__ Vector<T>& v)
{
        return sqrt (TBCI__ fabssqr (v));
}

INST(template <typename T> class TBCI__ Vector friend T abs (const TBCI__ Vector<T>&);)
template <typename T>
inline T abs (const TBCI__ Vector<T>& v)
{
        return v.abs();
}

INST(template <typename T> class TBCI__ TVector friend double fabs (TBCI__ TVector<T>);)
template <typename T>
inline double fabs (TBCI__ TVector<T> tv)
{ 
        TBCI__ Vector<T> v(tv); return fabs (v); 
}

INST(template <typename T> class TBCI__ Vector friend T abs (TBCI__ TVector<T>);)
template <typename T>
inline T abs (TBCI__ TVector<T> tv)
{ 
        TBCI__ Vector<T> v(tv); return abs (v); 
}

NAMESPACE_CSTD_END

NAMESPACE_TBCI

template <typename T>
class TSVector : public Vector_Sig<T> //: public TVector<T>
{
 protected:
   mutable T* vec; mutable unsigned long dim; //bool keep;
   mutable T fac ALIGN(MIN_ALIGN2);
        
   void clone (const bool = false, const T* = 0) const;
 public:
   void detach (const T* = 0) const;
   
 public:
   mutable bool mut;  // mutability (ie private data)
   void destroy () const;

 public:
   friend class Vector<T>;
   friend class TVector<T>;
   friend class CRMatrix<T>;
   friend NOINST void FRIEND_TBCI2__ do_mat_tsv_mult FGD (const unsigned start, const unsigned end, \
                        TVector<T> * res, const Matrix<T> * mat, const TSVector<T> * rsv);
   
   typedef T value_type;
   typedef T element_type;
   typedef T aligned_value_type TALIGN(MIN_ALIGN2);
   
   const TSVector<T>& eval (const T* vv = 0) const;
   
   TSVector () : vec((T*)0), dim(0), fac((T)1), mut(true)  {}
   ~TSVector ()  { /* keep = true; */   }
   explicit TSVector (const unsigned long d) : vec(new T[d]), dim(d), fac((T)1), mut(true)  {}
   
   TSVector (const TSVector<T>& ts)
     : vec(ts.vec), dim(ts.dim), fac(ts.fac), mut(ts.mut)  {}

   TSVector (const TVector<T>& tv, const T& f = (T)1) 
     : vec(tv.vec), dim(tv.dim), fac(f), mut(true)  {}
   TSVector (const Vector<T>& v, const T& f = (T)1)
     : vec(v.vec), dim(v.dim), fac(f), mut(false)  {}
   
        
   T operator () (const unsigned long i) const
     { T val = vec[i] * fac; destroy (); return val; }
   //T& operator () (const unsigned long i) { STD__ cerr << "TSVectors may not be written to by () op.\n"; abort(); }
   T operator [] (const unsigned long i) const { return this->operator() (i); }
   //T& operator [] (const unsigned long i) { return this->operator() (i); }
        
   T get (const unsigned long i) const HOT { return vec[i] * fac; }
   const T& getfac () const  { return fac; }
   T& getfacref () const  { return fac; }
   const T& getcref (const unsigned long i) const HOT { return vec[i]; }

   // alias equ ops
   TSVector<T>& operator = (const TSVector<T>& ts)
      { destroy (); vec = ts.vec; dim = ts.dim;
        fac = ts.fac; mut = ts.mut; return *this; }
   TSVector<T>& operator = (const TVector<T>& tv)
      { destroy (); vec = tv.vec; dim = tv.dim;
        fac = (T)1; mut = true; return *this; }
   /*
   TSVector<T>& operator = (const Vector<T>& v)
      { destroy (); vec = v.vec; dim = v.dim;
        fac = (T)1; mut = false; return *this; }
    */

   // first the easy to implement operators
   TSVector<T>& operator *= (const T& f)  { fac *= f; return *this; }
   TSVector<T>& operator /= (const T& f)  { fac /= f; return *this; }
   const TSVector<T>& operator *  (const T& f) const  { fac *= f; return *this; }
   const TSVector<T>& operator /  (const T& f) const  { fac /= f; return *this; }

   const TSVector<T>& operator - () const  { fac = -fac; return *this;  }

   // now for the more complicated ones
   TVector<T>  operator + (const Vector<T>&) const HOT;
   TVector<T>  operator + (TVector<T>) const HOT;
   TVector<T>  operator + (const TSVector<T>&) const HOT;
   TVector<T>  operator + (const T&) const;
   TVector<T>  operator - (const Vector<T>&) const HOT;
   TVector<T>  operator - (TVector<T>) const HOT;
   TVector<T>  operator - (const TSVector<T>&) const HOT;
   TVector<T>  operator - (const T&) const;
   
#ifdef TSV_OP_FRIENDS
   friend const TSVector<T>& operator * FGD (const T& f, const TSVector<T>& ts);
   friend TVector<T> operator + FGD (const T&, const TSVector<T>&);
   friend TVector<T> operator - FGD (const T&, const TSVector<T>&);
   // for those "friends"
#else
   TVector<T> add_t_tsv (const T&) const;
   TVector<T> sub_t_tsv (const T&) const;
#endif
   
   TVector<T> incr (const unsigned long, const T = (T)1) const;
  
   inline unsigned long size () const { return dim; }

   T min ()  { Vector<T> th(*this); return th.min (); }
   T max ()  { Vector<T> th(*this); return th.max (); }
   T abs () const /* NOTE: The const IS wrong */
        { return ((TVector<T>)*this).abs () * GLBL__ CSTD__ abs(fac); }
   double fabs () const 
        { return ((TVector<T>)*this).fabs () * GLBL__ MATH__ fabs(fac); }
   double fabssqr () const
        { return ((TVector<T>)*this).fabssqr () * GLBL2__ TBCI__ fabssqr(fac); }

   T sum ()  { T val = ((TVector<T>)*this).sum () * fac; destroy (); return val; }
   
   /* Why are those declared but not defined ? */
#ifdef TSVEC_NEED_DOT
   T operator * (const TSVector<T>&);
   T operator * (const TVector<T>&);
   T operator * (const Vector<T>&);
#endif

#if 0   /* not implemented! */
   friend TVector<T> conj FGD (const TSVector<T>&);
   friend TVector<T> real FGD (const TSVector<T>&);
   friend TVector<T> imag FGD (const TSVector<T>&);
#endif
   
// They fortunately also don't need to be friends 
#ifndef HAVE_GCC295_FRIEND_BUG
   friend NOINST double FRIEND_TBCI__ fabssqr FGDT (const TSVector<T>&);
#endif
#if 0   
   friend NOINST double MATH__ fabs FGD (const TBCI__ TSVector<T>&);
   friend NOINST T CSTD__ abs FGD (const TBCI__ TSVector<T>&);
#endif

   bool operator == (const Vector<T>&) const;
   bool operator != (const Vector<T>& v) const  { return !(*this == v); }
   bool operator == (const TSVector<T>&) const;
   bool operator != (const TSVector<T>& ts) const  { return !(*this == ts); }
   bool operator == (const TVector<T>& tv) const  { Vector<T> v(tv); return (*this == v); }
   bool operator != (const TVector<T>& tv) const { return !(*this == tv); }

  // Needed for emul et al.
  T* const & vecptr () const  { return vec; }

  static const char* vec_info()  { return "TSVector"; }

  friend STD__  ostream& operator << FGD (STD__  ostream&, const TSVector<T>&); 
} ;


template <typename T>
inline void TSVector<T>::destroy () const
{
        if (LIKELY(!mut || !dim)) return;
        TBCIDELETE(T, vec, dim);
}
   
template <typename T>
inline void TSVector<T>::detach (const T* vv) const
{
        if (LIKELY(!dim || (mut && !vv)))
                return;
        mut = true;
        if (vv) 
                vec = const_cast<T*> (vv); 
        else {
                vec = NEW(T,dim);
                if (UNLIKELY(!vec)) {
                        dim = 0; mut = false;
                }
        }
}

INST(template <typename T> class TSVector<T> friend const TSVector<T>& operator * (const T&, const TSVector<T>&);)
template <typename T>
const TSVector<T>& operator * (const T& f, const TSVector<T>& ts)
{ ts.getfacref() = f * ts.getfac(); return ts; }

template <typename T>
inline void TSVector<T>::clone (const bool evl, const T* vv) const
{
        if (LIKELY(!dim || (mut && !vv)))
                return;
        const T* oldv = vec; bool oldm = mut; detach (vv);
        if (LIKELY(evl && fac != (T)1)) {
#if 0 //def __i386__
                for (REGISTER unsigned long t = 0; t < dim; t++)
                        vec[t] = oldv[t] * fac;
#else      
                STD_SMP_TEMPLATE3VC(vec_val_mul, dim, vec, oldv, fac);
#endif     
                fac = (T)1;
        }
        else 
                TBCICOPY(vec, oldv, T, dim);
        if (UNLIKELY(vv && oldm))
                TBCIDELETE(T, oldv, dim);
}

template <typename T>
const TSVector<T>& TSVector<T>::eval (const T* vv) const
{
        if (!vv && mut) {
                if (fac != (T)1) {
#if 0 //def __i386__
                        for (REGISTER unsigned long t = 0; t < dim; t++)
                                vec[t] *= fac;
#else
                        STD_SMP_TEMPLATE2C(vec_mul_val, dim, vec, fac);
#endif
                        fac = (T)1;
                }
                // else do nothing
        }
        else 
                clone (true, vv);
        return *this;
}

template <typename T>
bool TSVector<T>::operator == (const Vector<T>& v) const
{
        if (LIKELY(dim != v.dim)) {
                destroy (); return false; 
        }
        if (LIKELY(dim == 0)) { 
                /*destroy();*/ return true; 
        }
        if (LIKELY(fac == (T)1)) {
                if (LIKELY(vec == v.vec))
                        return true;
                if (TBCICOMP(vec, v.vec, T, dim)) { 
                        destroy (); return false; 
                }
        } else {
                if (LIKELY(vec == v.vec))
                        return false;
                for (REGISTER unsigned long i = 0; i < dim; ++i)
                        if (UNLIKELY(vec[i]*fac != v.vec[i])) {
                                destroy (); return false; 
                        }
        }
        destroy (); 
        return true;
}

template <typename T>
bool TSVector<T>::operator == (const TSVector<T>& ts) const
{
        if (LIKELY(dim != ts.dim)) {
                ts.destroy (); destroy (); return false; 
        }
        if (LIKELY(dim == 0)) {
                /*destroy(); ts.destroy();*/ return true; 
        }
        if (LIKELY(fac == ts.fac)) {
                if (LIKELY(vec == ts.vec)) {
                        destroy (); return true; 
                }
                if (TBCICOMP(vec, ts.vec, T, dim)) { 
                        ts.destroy (); destroy (); return false; 
                }
        } else {
                if (LIKELY(vec == ts.vec)) { 
                        destroy (); return false; 
                }
                for (REGISTER unsigned long i = 0; i < dim; ++i)
                        if (UNLIKELY(vec[i]*fac != ts.vec[i]*ts.fac)) {
                                ts.destroy (); destroy (); return false; 
                        }
        }
        ts.destroy (); destroy (); return true;
}

template <typename T>
inline TVector<T>  TSVector<T>::operator + (TVector<T> tv) const
{
        BCHK(dim!=tv.dim,VecErr,Wrong dims in TSV + TV,tv.dim,tv);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  tv.set (fac*vec[t] + tv.vec[t], t);
        STD_SMP_TEMPLATE3VC(vec_add_svc, dim, tv.vec, vec, fac);
        destroy (); return tv;
}
template <typename T>
inline TVector<T>  TSVector<T>::operator - (TVector<T> tv) const
{
        BCHK(dim!=tv.dim,VecErr,Wrong dims in TSV - TV,tv.dim,tv);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  tv.set (fac*vec[t] - tv.vec[t], t);
        STD_SMP_TEMPLATE3VC(vec_sub_svc_inv, dim, tv.vec, vec, fac);
        destroy (); return tv;
}

template <typename T>
inline TVector<T> TSVector<T>::operator + (const Vector<T>& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else 
                tv.resize (dim);
        BCHK(dim!=v.dim,VecErr,Wrong dims in TSV + V,v.dim,tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = (vec[i]*fac + v.vec[i]); }
        STD_SMP_TEMPLATE4V(svc_vec_add, dim, tv.vec, vec, v.vec, fac);
        return (tv);
}
template <typename T>
inline TVector<T> TSVector<T>::operator - (const Vector<T>& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else 
                tv.resize (dim);
        BCHK(dim!=v.dim,VecErr,Wrong dims in TSV - V,v.dim,tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = (vec[i]*fac - v.vec[i]); }
        STD_SMP_TEMPLATE4V(svc_vec_sub, dim, tv.vec, vec, v.vec, fac);
        return (tv);
}

template <typename T>
inline TVector<T> TSVector<T>::operator + (const TSVector<T>& tsv) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else if (LIKELY(tsv.mut)) {
                tv.vec = tsv.vec; tv.dim = tsv.dim; 
        } else 
                tv.resize (dim);
        BCHK(dim!=tsv.dim,VecErr,Wrong dims in TSV + TSV,tsv.dim,tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = vec[i]*fac + tsv.vec[i]*tsv.fac; }
        STD_SMP_TEMPLATE5(svc_svc_add, dim, tv.vec, vec, tsv.vec, fac, tsv.fac);
        if (UNLIKELY(mut && tsv.mut)) 
                tsv.destroy (); 
        return tv;
}
template <typename T>
inline TVector<T> TSVector<T>::operator - (const TSVector<T>& tsv) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else if (LIKELY(tsv.mut)) { 
                tv.vec = tsv.vec; tv.dim = tsv.dim; 
        } else 
                tv.resize (dim);
        BCHK(dim!=tsv.dim,VecErr,Wrong dims in TSV - TSV,tsv.dim,tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = vec[i]*fac - tsv.vec[i]*tsv.fac; }
        STD_SMP_TEMPLATE5(svc_svc_sub, dim, tv.vec, vec, tsv.vec, fac, tsv.fac);
        if (UNLIKELY(mut && tsv.mut)) 
                tsv.destroy ();
        return tv;
}

template <typename T>
inline TVector<T> TSVector<T>::operator + (const T& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else 
                tv.resize (dim);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = vec[i]*fac op v; }
        STD_SMP_TEMPLATE4C(svc_val_add, dim, tv.vec, vec, fac, v);
        return tv;
}
template <typename T>
inline TVector<T> TSVector<T>::operator - (const T& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        } else 
                tv.resize (dim);
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = vec[i]*fac op v; }
        STD_SMP_TEMPLATE4C(svc_val_sub, dim, tv.vec, vec, fac, v);
        return tv;
}

template <typename T>
inline TVector<T> TSVector<T>::add_t_tsv (const T& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        }
        else 
                tv.resize (size());
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = (v op vec[i]*fac); }
        STD_SMP_TEMPLATE4C(val_svc_add, dim, tv.vec, vec, v, fac);
        return tv;
}
INST(template <typename T> class TSVector<T> friend TVector<T> operator + (const T&, const TSVector<T>&);)
template <typename T>
inline TVector<T> operator + (const T& v, const TSVector<T>& tsv)
{ return tsv.add_t_tsv (v); }

template <typename T>
inline TVector<T> TSVector<T>::sub_t_tsv (const T& v) const
{
        TVector<T> tv;
        if (LIKELY(mut)) { 
                tv.vec = vec; tv.dim = dim; 
        }
        else 
                tv.resize (size());
        //for (REGISTER unsigned long i=0; i<dim; ++i)
        //  { tv.vec[i] = (v op vec[i]*fac); }
        STD_SMP_TEMPLATE4C(val_svc_sub, dim, tv.vec, vec, v, fac);
        return tv;
}
INST(template <typename T> class TSVector<T> friend TVector<T> operator - (const T&, const TSVector<T>&);)
template <typename T>
inline TVector<T> operator - (const T& v, const TSVector<T>& tsv)
{ return tsv.sub_t_tsv (v); }

template <typename T>
STD__  ostream& operator << (STD__  ostream& os, const TSVector<T>& ts)
{
        for (unsigned long i = 0; i < ts.dim; ++i) 
                os << ts.vec[i]*ts.fac << " ";
        ts.destroy (); return os;
}


template <typename T>
inline TVector<T> TSVector<T>::incr (const unsigned long wh, const T i) const
{
        TVector<T> tv (*this);
        BCHK (wh >= dim, VecErr, Tried to change not-exitsting idx, wh, tv);
        tv.vec[wh] += i; /* destroy(); */
        return tv;
}


INST(template <typename T> class TSVector friend double fabssqr (const TSVector<T>&);)
template <typename T>
double fabssqr (const TSVector<T>& ts)
{ return ts.fabssqr (); }

NAMESPACE_END

NAMESPACE_CSTD

INST(template <typename T> class TBCI__ TSVector friend double fabs (const TBCI__ TSVector<T>&);)
template <typename T>
double inline fabs (const TBCI__ TSVector<T>& ts)
{ return ts.fabs(); }

INST(template <typename T> class TBCI__ TSVector friend T abs (const TBCI__ TSVector<T>&);)
template <typename T>
inline T abs (const TBCI__ TSVector<T>& ts)
{ return ts.abs(); }

NAMESPACE_CSTD_END

NAMESPACE_TBCI

//------------------------------------------------------------

template <typename T>
class Vector : public TVector<T> //, public Vector_Sig<T>
{
 public:
  //typedef T TALIGN(sizeof(T)) vec_t;
  typedef T value_type;
  typedef T element_type;
  typedef T aligned_value_type TALIGN(MIN_ALIGN2);
  
  friend class TSVector<T>;
  // constructors are NEVER inherited, 
  // so we inherit them explicitly 
  // (or do the default constructors do exactly the same ?) 
  
  explicit Vector (const unsigned long d = 0) : TVector<T> (d)  {}
  Vector (const T& val, const unsigned long d) : TVector<T> (val, d)  {}
  Vector (const BVector<T>& bv) : TVector<T> (bv)  {}
  Vector (const Vector<T>& v) : TVector<T> (v)  {}
  //alias
  Vector (const TSVector<T>& ts)
     { ts.eval (this->vec); this->vec = ts.vec; this->dim = ts.dim; }
  Vector (const TVector<T>& tv) : TVector<T> (tv)  {} //remove const ?
  Vector (const Mat_Brack<T>&);
#ifndef NO_POD
  Vector (const vararg va, ... );
#endif
  //disallow automatic conversion to detect errors in the code
  // KG, 99/05/19: 
  // Danger of recursion: int -> cplx<int> -> cplx<cplx<int> > ... compiler death
  // Strange enough, egcs can handle it, if no namespaces are used (!)
  #if defined(__GNUC__) && defined VEC_CPLX_CONSTR
    explicit Vector (const Vector<cplx<T> >& cv);
  #endif

  ~Vector() { this->destroy (); }

   /* Read-only access */
   typename tbci_traits<T>::const_refval_type
      operator () (const unsigned long i) const
        { return BVector<T>::operator () (i); }
   typename tbci_traits<T>::const_refval_type
      operator [] (const unsigned long i) const 
        { return this->operator() (i); }
   // Read-write access
   T& operator () (const unsigned long i)
        { return BVector<T>::operator () (i); }
   T& operator [] (const unsigned long i) 
        { return this->operator() (i); }
   
   //T& operator () (const unsigned long i) const
     //{ return BVector<T>::operator () (i); }

   Vector<T>& fill (const T& v) { BVector<T>::fill (v); return *this; }
   Vector<T>& clear ()          { return this->fill((T)0); }
        
  // Ahhh! Slicing!
   TVector<T> slice (unsigned long, unsigned long) const;

   TVector<T> incr (const unsigned long, const T = (T)1) const;
    
  // assignment operators

  Vector<T>& operator = (const T& a) 
    { this->TVector<T>::operator = (a); return *this; }
  Vector<T>& operator = (const Vector<T>& v) HOT
    { this->TVector<T>::operator = (v); return *this; }
  Vector<T>& operator = (const TSVector<T>& ts) HOT
     { TVector<T>::operator = (ts); return *this; }
  Vector<T>& operator = (const TVector<T>& tv) HOT
    { this->TVector<T>::operator = (tv); return *this; }
#ifdef TBCI_NEW_BRACKET
  Vector<T>& operator = (const T* ptr)
    { TBCICOPY(this->vec,ptr,T,this->dim); return *this; }
#endif
   
  // The ones from TVector
  bool operator == (const Vector<T>& v) const
    { return BVector<T>::operator == ((BVector<T>)v); }
  bool operator != (const Vector<T>& v) const  { return !(*this == v); }
  bool operator == (const TSVector<T>&) const;
  bool operator != (const TSVector<T>& ts) const  { return !(*this == ts); }
  bool operator == (const TVector<T>& tv) const  { return (tv == *this); }
  bool operator != (const TVector<T>& tv) const  { return !(*this == tv); }
   
  bool operator <= (const Vector<T>& v) const  { return BVector<T>::operator <= (v); }
  bool operator >= (const Vector<T>& v) const  { return BVector<T>::operator >= (v); }
  bool operator <  (const Vector<T>& v) const  { return !((*this) >= v); }
  bool operator >  (const Vector<T>& v) const  { return !((*this) <= v); }
        
  TVector<T>  operator + (const Vector<T>&) const HOT;
  TVector<T>  operator - (const Vector<T>&) const HOT;
  TVector<T>  operator + (/*const*/ TVector<T>/*&*/ tv) const HOT;
  TVector<T>  operator - (/*const*/ TVector<T>/*&*/ tv) const HOT;
  TVector<T>  operator + (const TSVector<T>&) const HOT;
  TVector<T>  operator - (const TSVector<T>&) const HOT; 

  TVector<T>  operator + (const T&) const;
  TVector<T>  operator - (const T&) const;
  TSVector<T> operator * (const T&) const;
  TSVector<T> operator / (const T&) const;

  TVector<T> operator - () const;

  // Additional

  friend T dot FGD (const Vector<T>&, const Vector<T>&) HOT;

  T operator * (const Vector<T>&) const;
  T operator * (TVector<T>& tv) const
    { Vector<T> v (tv); return this->operator * (v); }
  T operator / (const Vector<T>&) const; //NB: a/b := a*b/b^2

  T min () const;
  T max () const;
  T min (unsigned long& pos) const;
  T max (unsigned long& pos) const;
  T sum () const;
  double fabssqr () const HOT;
  double fabs () const  { return MATH__ sqrt (this->fabssqr()); }
  T abs () const { return (T) MATH__ sqrt (this->fabssqr()); }
  /* { T res(0.0); do_vec_sumsqr<T> (dim, vec, res); return res; } */

  bool contains (const T& v) const { return BVector<T>::contains (v); }

  //friend TVector<T> operator + FGD (const T&, const Vector<T>&);
  //friend TVector<T> operator - FGD (const T&, const Vector<T>&);
  //friend TSVector<T> operator * FGD (const T&, const Vector<T>&);

  static const char* vec_info()  { return "Vector"; }
  // useful?
  //Vector<T> apply(T op(T)); 
  
#if 0
  friend STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v);
#endif
  /* An inheritance problem in GCC */
#ifdef OLD_GCC
  friend STD__ istream& operator >> (STD__ istream& is, Vector<T>& v)
    {  return is >> (BVector<T>&)v;  }
#endif

#ifndef HAVE_PROMOTION_BUG
  // Promotion (only explicit)
//# ifndef HAVE_GCC295_TMPLFRNDCLS_BUG
//  template <typename U> friend class BVector;
//# endif
  template <typename U> explicit Vector (const BVector<U>& bv) : TVector<T> (bv) {}
#endif  

} ;

#if 1
INST(template <typename T> friend inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v);)
template <typename T>
inline STD__ ostream& operator << (STD__ ostream& os, const Vector<T>& v)
{  return os << (BVector<T>&)v;  }
#endif

#ifndef NO_POD
template <typename T>
inline Vector<T>::Vector (const vararg va, ...)
{         
        this->vec = (T*)0; this->dim = (unsigned)va; this->keep = false;
        if (UNLIKELY(this->dim)) {
                this->vec = NEW (T, this->dim);
                if (UNLIKELY(!this->vec))
                        this->dim = 0;
        }
        va_list vl;  
        va_start (vl, va);
#if !defined(__clang__) || !defined(CPLX)
        for (unsigned long i=0; i < this->dim; ++i) 
                this->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

#if defined(__GNUC__) && defined VEC_CPLX_CONSTR
template <typename T>
inline Vector<T>::Vector (const Vector<cplx<T> >& vc) : TVector<T> (vc.size())
{
        for (unsigned long i = 0; i < this->dim; ++i)
                this->vec[i] = CPLX__ real (vc(i));
}
#endif

template <typename T>
bool Vector<T>::operator == (const TSVector<T>& tsv) const
{
        if (LIKELY(this->dim != tsv.dim)) { 
                tsv.destroy (); return false; 
        }
        if (LIKELY(this->dim == 0)) { 
                /*tsv.destroy();*/ return true; 
        }
        if (LIKELY(tsv.fac == (T)1)) {
                if (LIKELY(this->vec == tsv.vec)) 
                        return true;
                if (TBCICOMP(this->vec, tsv.vec, T, this->dim))  { 
                        tsv.destroy (); return false; 
                }
        } else {
                if (LIKELY(this->vec == tsv.vec))
                        return false;
                for (REGISTER unsigned long i = 0; i < this->dim; ++i)
                        if (UNLIKELY(tsv.vec[i]*tsv.fac != this->vec[i])) { 
                                tsv.destroy (); return false; 
                        }
        }
        tsv.destroy (); return true;
}


#if 0
template <typename T>
inline Vector<T>& Vector<T>::operator *= (const T& f)
{
        for (unsigned long i = 0; i < this->dim; ++i) 
                this->vec[i] *= f;
        return *this;
}

template <typename T>
inline Vector<T>& Vector<T>::operator /= (const T& f)
{
        for (unsigned long i = 0; i < this->dim; ++i) 
                this->vec[i] /= f;
        return *this;
}
#endif

#if !defined(SMP) || defined(NOSMP_VECSCALAR)

// Scalar product with conj
template <typename T>
inline T dot (const Vector<T>& a, const Vector<T>& b)
{
        BCHK(a.dim != b.dim, VecErr, Scalar dot product with diff. size Vectors, b.size(), T(0));
        T s = T(0.0);
        if (do_exactsum())
                do_vec_dot_exact <T> (a.dim, a.vec, b.vec, s);
        else
                do_vec_dot_quick <T> (a.dim, a.vec, b.vec, s);
        //STD__ cout << "dot: " << s << STD__ endl;
# ifdef SMPVEC_DEBUG
        STD__ cerr << "dot: " << STD__ setprecision (16) << s << STD__ endl;
# endif
        return s;
}

// Scalar product withOUT conj
template <typename T>
inline T Vector<T>::operator * (const Vector<T>& a) const
{
        BCHK(this->dim != a.dim, VecErr, Scalar product with diff. size Vectors, a.size(), 0);
        T s = T(0.0);
        if (do_exactsum())
                do_vec_mult_exact <T> (this->dim, this->vec, a.vec, s);
        else
                do_vec_mult_quick <T> (this->dim, this->vec, a.vec, s);
# ifdef SMPVEC_DEBUG
        STD__ cerr << "op *: " << STD__ setprecision (16) << s << STD__ endl;
# endif
        return s;
}


#else /* SMP */

INST(template <typename T> class TVector friend void job_vec_dot (struct thr_ctrl*);)
template <typename T>
void job_vec_dot (struct thr_ctrl *tc)
{
        if (sizeof(T) <= 16) /* tc->t_par[2] == 0) */
                tc->t_par[2] = &tc->t_res_ptr;
        *(T*)(tc->t_par[2]) = T(0.0);
        if (do_exactsum())
                do_vec_dot_exact (tc->t_size, (const T*)(tc->t_par[0]),
                    (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
        else
                do_vec_dot_quick (tc->t_size, (const T*)(tc->t_par[0]),
                    (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
}

// Scalar product with conj
template <typename T>
inline T dot (const Vector<T>& a, const Vector<T>& b)
{
        BCHK(a.size() != b.size(), VecErr, Scalar dot product with diff. size Vectors, b.size(), 0);
        T res = T(0); unsigned long sz = a.size();
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        if (LIKELY(n_thr < 2)) {
                if (do_exactsum())
                        do_vec_dot_exact<T> (sz, a.vec, b.vec, res);
                else
                        do_vec_dot_quick<T> (sz, a.vec, b.vec, res);
                //STD__ cout << "dot: " << res << STD__ endl;
        }
        else
        {
                PREFETCH_R(a.vec, 2); PREFETCH_R(b.vec, 2);
                const unsigned long first = slice_offset(1, n_thr, sz, a.vec);
                unsigned long st, en = first;
                T comp = 0.0;
                T* results = (T*)0;
                if (sizeof(T) > 16)
                        results = new T[n_thr];
                //_tbci_fill<T>(n_thr, results, T(0.0));
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, a.vec);
                        thread_start (t, (thr_job_t)job_vec_dot<T>, en - st,
                                      a.vec+st, b.vec+st, results+t, (void*)0);
                }
                /* The first slice is handled by the main thread */
                if (do_exactsum())
                        do_vec_dot_exact<T> (first, a.vec, b.vec, res);
                else
                        do_vec_dot_quick<T> (first, a.vec, b.vec, res);
                //sched_yield ();
                /* Wait for the end */
                //STD__ cout << "dot: " << res;
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        job_output out;
                        thread_wait(t, &out);
                        T y;
                        if (sizeof(T) > 16)
                                y = results[t];
                        else {
                                const T* thrres = (const T*)(&out.t_res_ptr);
                                y = *thrres;
                        }
                        T tmp = res + y;
                        if (do_exactsum())
                                comp += (tmp - res) - y;
                        res = tmp;
#ifdef SMPVEC_DEBUG
                        STD__ cerr << "+" 
                                   << (sizeof(T) > THREAD_MAX_RES_LN? 
                                       results[t]: *((T*)(threads[t].t_res_dummy)));
#endif
                }
                //STD__ cout << ":" << res << STD__ endl;
                if (results)
                        delete[] results;
                if (do_exactsum())
                        res -= comp;
        }
#ifdef SMPVEC_DEBUG
        STD__ cerr << "dot: " << STD__ setprecision (16) << res << STD__ endl;
#endif
        return res;
}

INST(template <typename T> class TVector friend void job_vec_mult (struct thr_ctrl*);)
template <typename T>
void job_vec_mult (struct thr_ctrl *tc)
{
        *(T*)(tc->t_par[2]) = T(0.0);
        if (do_exactsum())
                do_vec_mult_exact (tc->t_size, (const T*)(tc->t_par[0]),
                             (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
        else
                do_vec_mult_quick (tc->t_size, (const T*)(tc->t_par[0]),
                             (const T*)(tc->t_par[1]), *(T*)tc->t_par[2]);
}

// Scalar product withOUT conj
template <typename T>
inline T Vector<T>::operator * (const Vector<T>& a) const
{
        BCHK(this->dim != a.size(), VecErr, Scalar product with diff. size Vectors, a.size(), 0);
        T res = T(0); unsigned long sz = this->dim;
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        if (LIKELY(n_thr < 2)) {
                if (do_exactsum())
                        do_vec_mult_exact<T> (sz, this->vec, a.vec, res);
                else
                        do_vec_mult_quick<T> (sz, this->vec, a.vec, res);
                //STD__ cout << "mult: " << res << STD__ endl;
        }
        else
        {
                PREFETCH_R(this->vec, 2); PREFETCH_R(a.vec, 2);
                const unsigned long first = sz/n_thr - (sz/n_thr)%4;
                unsigned long st, en = first;
                T comp(0.0);
                T* results = new T[n_thr];
                //_tbci_fill<T>(n_thr, results, T(0.0));
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = (t+2)*sz / n_thr;
                        if (en < n_thr-2)
                                en -= en%4;
                        thread_start (t, (thr_job_t)job_vec_mult<T>, en - st,
                                      this->vec+st, a.vec+st, results+t, (void*)0);
                }
                /* The last slice is handled by the main thread */
                if (do_exactsum())
                        do_vec_dot_exact<T> (first, this->vec, a.vec, res);
                else
                        do_vec_dot_quick<T> (first, this->vec, a.vec, res);
                //sched_yield ();
                /* Wait for the end */
                //STD__ cout << "dot: " << res;
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        thread_wait (t);
                        T tmp = res + results[t];
                        if (do_exactsum())
                                comp += (tmp - res) - results[t];
                        res = tmp;
#ifdef SMPVEC_DEBUG
                        STD__ cerr << "+"
                                   << (sizeof(T) > THREAD_MAX_RES_LN? 
                                       results[t]: *((T*)(threads[t].t_res_dummy)));
#endif
                }
                //STD__ cout << ":" << res << STD__ endl;
                delete[] results;
                if (do_exactsum())
                        res -= comp;
        }
#ifdef SMPVEC_DEBUG
        STD__ cerr << "mlt: " << STD__ setprecision (16) << res << STD__ endl;
#endif
        return res;
}

#endif /* SMP */



HOTDECL(template <typename T> 
bool par_comp (const Vector<T>& v1, const Vector<T>& v2))
{               
        const unsigned long dim = v1.size();
        if (LIKELY(dim != v2.size()))
                return false;
        if (LIKELY(!dim))
                return true;
        if (LIKELY(v1(0) != v2(0)))
                return false;
        return !(_par_comp<T>(dim, &v1.getcref(0), &v2.getcref(0)));
}

HOTDECL(template <typename T> 
void par_copy (Vector<T>& v1, const Vector<T>& v2))
{               
        const unsigned long dim = v1.size();
        BCHK(dim != v2.size(), VecErr, Copy with wrong size, dim, v2.size());
        if (LIKELY(!dim))
                return;
        _par_copy<T>(dim, &v1.setval(0), &v2.getcref(0));
}

INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, const Vector<T>&);)
INST(template <typename T> class TVector friend bool par_comp (const Vector<T>&, TVector<T>);)
INST(template <typename T> class TVector friend bool par_comp (TVector<T>, const Vector<T>&);)
INST(template <typename T> class TVector friend bool par_comp (TVector<T>, TVector<T>);)

template <typename T> 
inline bool par_comp (const Vector<T>& v1, TVector<T> v2)
{               
        return par_comp (v1, (const Vector<T>)v2);
        // Note: The cast aliases the TV and the mem is freed ;-)
}

template <typename T> 
inline bool par_comp (TVector<T> v1, const Vector<T>& v2)
{               
        return par_comp ((const Vector<T>)v1, v2);
}

template <typename T> 
inline bool par_comp (TVector<T> v1, TVector<T> v2)
{               
        return par_comp ((const Vector<T>)v1, (const Vector<T>)v2);
}

template <typename T>
inline void par_fill(Vector<T>& vec, vec_fill_fn<T> fn, void* par)
{
        const unsigned long dim = vec.size();
        if (!dim)
                return;
        _par_fill_fn(dim, &vec.setval(0), fn, par);
}

template <typename T> 
T Vector<T>::operator / (const Vector<T>& a) const
{
        double div ALIGN(MIN_ALIGN) = a.fabssqr();
        BCHK(div == 0, VecErr, Division by Vector with norm 0, 0, div);
        return (T)((*this * a) / div); 
}

template <typename T>
T Vector<T>::min () const
{
        if (LIKELY(this->dim == 0))
                return 0;
        T min ALIGN(MIN_ALIGN) = this->vec[0];
        for (REGISTER unsigned long t = 1; t < this->dim; t++)
                if (UNLIKELY(this->vec[t] < min))  
                        min = this->vec[t];
        return min;
}

template <typename T>
T Vector<T>::min (unsigned long& pos) const
{
        if (LIKELY(this->dim == 0))
                return 0;
        BCHK (pos >= this->dim, VecErr, pos outside range, pos, 0);
        T min ALIGN(MIN_ALIGN) = this->vec[pos];
        for (REGISTER unsigned long t = pos+1; t < this->dim; t++)
                if (UNLIKELY(this->vec[t] < min)) {
                        min = this->vec[t]; pos = t; 
                }
        return min;
}

template <typename T>
T Vector<T>::max () const
{
        if (LIKELY(this->dim == 0))
                return 0;
        T max ALIGN(MIN_ALIGN) = this->vec[0];
        for (REGISTER unsigned long t = 1; t < this->dim; t++)
                if (UNLIKELY(this->vec[t] > max))
                        max = this->vec[t]; 
        return max;
}

template <typename T>
T Vector<T>::max (unsigned long& pos) const
{
        if (LIKELY(this->dim == 0))
                return 0;
        BCHK (pos >= this->dim, VecErr, pos outside range, pos, 0);
        T max ALIGN(MIN_ALIGN) = this->vec[pos];
        for (REGISTER unsigned long t = pos+1; t < this->dim; t++)
                if (UNLIKELY(this->vec[t] > max)) { 
                        max = this->vec[t]; pos = t;
                }
        return max;
}

#ifdef SMP
INST(template <typename T> class Vector friend void job_vec_sum(struct thr_ctrl*);)
template <typename T>
void job_vec_sum(struct thr_ctrl *tc)
{
        *(T*)tc->t_par[1] = 0.0;
        if (do_exactsum())
                do_vec_sum_exact<T>(tc->t_size, (const T*)tc->t_par[0], *(T*)tc->t_par[1]);
        else
                do_vec_sum_quick<T>(tc->t_size, (const T*)tc->t_par[0], *(T*)tc->t_par[1]);
}

template <typename T>
T Vector<T>::sum() const
{
        T res = 0.0; const unsigned long sz = this->dim;
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        if (LIKELY(n_thr < 2)) {
                if (do_exactsum())
                        do_vec_sum_exact<T> (sz, this->vec, res);
                else
                        do_vec_sum_quick<T> (sz, this->vec, res);
        }
        else {
                PREFETCH_R(this->vec, 2);
                const unsigned long first = slice_offset(1, n_thr, sz, this->vec);
                unsigned long st, en = first;
                T comp = 0.0;
                T* results = new T[n_thr];
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, this->vec);
                        thread_start (t, (thr_job_t)job_vec_sum<T>, en-st,
                                      this->vec+st, results+t, (void*)0);
                }
                /* The first slice is handled by the main thread */
                if (do_exactsum())
                        do_vec_sum_exact<T> (first, this->vec, res);
                else
                        do_vec_sum_quick<T> (first, this->vec, res);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        thread_wait (t);
                        T tmp = res + results[t];
                        if (do_exactsum())
                                comp += (tmp - res) - results[t];
                        res = tmp;
                }
                delete[] results;
                if (do_exactsum())
                        res -= comp;
        }
        return res;
}

#else
template <typename T>
T Vector<T>::sum () const
{
        T sum ALIGN(MIN_ALIGN) = T(0);
        if (do_exactsum())
                do_vec_sum_exact<T> (this->dim, this->vec, sum);
        else
                do_vec_sum_quick<T> (this->dim, this->vec, sum);
        return sum;
}
#endif

#if 0
{
  if (!this->dim) return 0.0;
  double s = CPLX__ real (vec[0] * CPLX__ conj(vec[0]));
  for (unsigned long i = 1; i < this->dim; ++i)
    s += CPLX__ real ((this->vec[i] * CPLX__ conj(this->vec[i])));
  return (sqrt(s));
//  STD__ cout << s << STD__ endl;
//  STD__ cout << (*this) * (*this) << STD__ endl;
//  return MATH__ fabs( MATH__ sqrt (CPLX__ real(((*this) * (*this))))); 
}
#endif

template <typename T>
inline TVector<T> Vector<T>::operator + (const Vector<T>& v) const
{
        TVector<T> tv (this->dim);
        BCHK(this->dim!=v.dim, VecErr, Unequal size in V + V, v.dim, tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] + v.vec[i];
        STD_SMP_TEMPLATE3VV(vec_vec_add, this->dim, tv.vec, this->vec, v.vec);
        return tv;
}
template <typename T>
inline TVector<T> Vector<T>::operator - (const Vector<T>& v) const
{
        TVector<T> tv (this->dim);
        BCHK(this->dim!=v.dim, VecErr, Unequal size in V - V, v.dim, tv);
        //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] - v.vec[i];
        STD_SMP_TEMPLATE3VV(vec_vec_sub, this->dim, tv.vec, this->vec, v.vec);
        return tv;
}

template <typename T>
inline TVector<T> Vector<T>::operator + (/*const*/ TVector<T>/*&*/ tv) const
{
        BCHK(this->dim!=tv.dim, VecErr, Unequal size in V + TV, tv.dim, tv);
        //for (REGISTER unsigned long t=0; t<dim; t++) tv.vec[t] = vec[t] + tv.vec[t];
        STD_SMP_TEMPLATE2V(vec_add_vec, this->dim, tv.vec, this->vec);
        return tv;
} 
template <typename T>
inline TVector<T> Vector<T>::operator - (/*const*/ TVector<T>/*&*/ tv) const
{
        BCHK(this->dim!=tv.dim, VecErr, Unequal size in V - TV, tv.dim, tv);
        //for (REGISTER unsigned long t=0; t<dim; t++) tv.vec[t] = vec[t] - tv.vec[t];
        STD_SMP_TEMPLATE2V(vec_sub_vec_inv, this->dim, tv.vec, this->vec);
        return tv;
}

// No need to destroy ts, as it's either an alias to a Vector (=> non-mutable)
// or an alias to the output TVector
template <typename T>
inline TVector<T> Vector<T>::operator + (const TSVector<T>& ts) const
{
        TVector<T> tv;
        BCHK(this->dim!=ts.dim, VecErr, Wrong dim in V + TSV, ts.dim, tv);
        if (LIKELY(ts.mut)) { 
                tv.vec = ts.vec; tv.dim = ts.dim; 
        } else 
                tv.resize (this->dim);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  tv.vec[t] = vec[t] + ts.fac * ts.vec[t];
        STD_SMP_TEMPLATE4V(vec_svc_add, this->dim, tv.vec, this->vec, ts.vec, ts.fac);
        return tv;
}
template <typename T>
inline TVector<T> Vector<T>::operator - (const TSVector<T>& ts) const
{
        TVector<T> tv;
        BCHK(this->dim!=ts.dim, VecErr, Wrong dim in V - TSV, ts.dim, tv);
        if (LIKELY(ts.mut)) { 
                tv.vec = ts.vec; tv.dim = ts.dim; 
        } else 
                tv.resize (this->dim);
        //for (REGISTER unsigned long t = 0; t < dim; t++)
        //  tv.vec[t] = vec[t] - ts.fac * ts.vec[t];
        STD_SMP_TEMPLATE4V(vec_svc_sub, this->dim, tv.vec, this->vec, ts.vec, ts.fac);
        return tv;
}

template <typename T>
inline TVector<T> Vector<T>::operator + (const T& a) const
{
        TVector<T> tv (this->dim);
        //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] + a;
        STD_SMP_TEMPLATE3VC(vec_val_add, this->dim, tv.vec, this->vec, a);
        return tv;
}
template <typename T>
inline TVector<T> Vector<T>::operator - (const T& a) const
{
        TVector<T> tv (this->dim);
        //for (REGISTER unsigned long i=0; i<dim; ++i) tv.vec[i] = vec[i] - a;
        STD_SMP_TEMPLATE3VC(vec_val_sub, this->dim, tv.vec, this->vec, a);
        return tv;
}


template <typename T>
inline TSVector<T> Vector<T>::operator * (const T& a) const
{
        return TSVector<T> (*this, a);
}

template <typename T>
inline TSVector<T> Vector<T>::operator / (const T& a) const
{
        BCHK (a==(T)0, VecErr, Division by zero in TSVec / T, 0, TSVector<T> (*this, 1));
        return TSVector<T> (*this, (T)1 / a);
}


INST(template <typename T> class Vector friend TVector<T> operator + (const T&, const Vector<T>&);)
template <typename T>
inline TVector<T> operator + (const T& a, const Vector<T>& v)
{
        const unsigned long dim = v.size();
        TVector<T> tv (dim);
        const T* const vvec(&v.getcref(0));
        T* const tvvec(&tv.setval(0));
        //for (REGISTER unsigned long i=0; i<v.size(); ++i)
        //      tv.set (a + v.get(i), i);
        STD_SMP_TEMPLATE3VC(val_vec_add, dim, tvvec, vvec, a);
        return tv;
}
INST(template <typename T> class Vector friend TVector<T> operator - (const T&, const Vector<T>&);)
template <typename T>
inline TVector<T> operator - (const T& a, const Vector<T>& v)
{
        const unsigned long dim = v.size();
        TVector<T> tv (dim);
        const T* const vvec(&v.getcref(0));
        T* const tvvec(&tv.setval(0));
        //for (REGISTER unsigned long i=0; i<v.size(); ++i)
        //      tv.set (a - v.get(i), i);
        STD_SMP_TEMPLATE3VC(val_vec_sub, dim, tvvec, vvec, a);
        return tv;
}


INST(template <typename T> class Vector friend TSVector<T> operator * (const T&, const Vector<T>&);)
template <typename T>
inline TSVector<T> operator * (const T& a, const Vector<T>& b)
{
        return TSVector<T> (b, a);
}

// a/b = a*b/b^2
INSTCTL(#ifndef VEC_INT_DIV_SUPPORT)
INST(template <typename T> class Vector friend TSVector<T> operator / (const T&, const Vector<T>&);)
INSTCTL(#endif)
template <typename T>
TSVector<T> operator / (const T& a, const Vector<T>& b)
{
        T dv = a / fabssqr (b);
        BCHK(dv==(T)0, VecErr, Divide by norm 0 Vector, 0, 0);
        return TSVector<T> (b, dv);
}

template <typename T>
inline TVector<T> Vector<T>::operator - () const
{
        TVector<T> tv (this->dim);
        do_vec_neg_vec<T> (this->dim, tv.vec, this->vec);
        return tv;
}


template <typename T>
inline TVector<T> Vector<T>::slice (unsigned long i0, unsigned long i1) const
{
        BCHK (i1 < i0, VecErr, Slicing needs ordered arguments, i0, *this);
        BCHK (i0 >= this->dim, VecErr, Slice: arg1 out of range, i0, *this);
        BCHK (i1 > this->dim, VecErr, Slice: arg2 out of range, i1, *this);
        TVector<T> tv (i1 - i0);
        if (UNLIKELY(tv.dim))
                TBCICOPY(tv.vec, this->vec+i0, T, tv.dim);
        return tv;
}


template <typename T>
inline TVector<T> Vector<T>::incr (const unsigned long wh, const T i) const
{
        TVector<T> tv (*this);
        BCHK (wh >= this->dim, VecErr, Tried to change not-exitsting idx, wh, tv);
        tv.vec[wh] += i;
        return tv;
}


/* Cost for these ops with Vecs of dim d */
#define COST_VECFABSSQR(d) (d*(COST_UNIT_LOAD+COST_MULT+COST_ADD+COST_LOOP))
#define COST_VECSCALAR(d) (d*(2*COST_UNIT_LOAD+COST_MULT+COST_ADD+COST_LOOP))

//INST(template <typename T> class Vector friend void do_vec_fabssqr (const unsigned long, const T*, double &);)

INST(template <typename T> class Vector friend double fabssqr (const Vector<T>&);)
template <typename T>
inline double fabssqr (const Vector<T>& v)
{ return v.fabssqr(); ;}


#if !defined(SMP) || defined(NOSMP_VECFABS)

template <typename T>
inline double Vector<T>::fabssqr () const
{
        double res = 0.0;
        if (do_exactsum())
                do_vec_fabssqr_exact (this->dim, this->vec, res);
        else
                do_vec_fabssqr_quick (this->dim, this->vec, res);
#ifdef SMPVEC_DEBUG
        STD__ cerr << "fabssqr: " << STD__ setprecision (16) << res << STD__ endl;
#endif
        return res;
}

#else /* SMP */

INST(template <typename T> class Vector friend void job_vec_fabssqr (struct thr_ctrl*);)
template <typename T>
void job_vec_fabssqr (struct thr_ctrl *tc)
{
        tc->t_res_d = 0.0;
        if (do_exactsum())
                do_vec_fabssqr_exact (tc->t_size, (const T*)(tc->t_par[0]), tc->t_res_d);
        else
                do_vec_fabssqr_quick (tc->t_size, (const T*)(tc->t_par[0]), tc->t_res_d);
}

template <typename T>
inline double Vector<T>::fabssqr () const
{
        double res = 0.0; unsigned long sz = this->dim;
        /* use some heuristic to decide for the num of threads */
        const unsigned n_thr = threads_avail (sz / SMP_VECSLICE2);
        update_n_thr(n_thr);
        if (LIKELY(n_thr < 2)) {
                if (do_exactsum())
                        do_vec_fabssqr_exact<T> (sz, this->vec, res);
                else
                        do_vec_fabssqr_quick<T> (sz, this->vec, res);
        }
        else {
                PREFETCH_R(this->vec, 2);
                const unsigned long first = slice_offset(1, n_thr, sz, this->vec);
                unsigned long st, en = first;
                double comp(0.0);
                /* Start threads */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        st = en; en = slice_offset(t+2, n_thr, sz, this->vec);
                        thread_start (t, (thr_job_t)job_vec_fabssqr<T>, en-st,
                                      this->vec+st, (void*)0);
                }
                /* The first slice is handled by the main thread */
                if (do_exactsum())
                        do_vec_fabssqr_exact<T> (first, this->vec, res);
                else
                        do_vec_fabssqr_quick<T> (first, this->vec, res);
                //sched_yield ();
                /* Wait for the end */
                for (unsigned t = 0; t < n_thr-1; ++t) {
                        double y = thread_wait_result (t);
                        double tmp = res + y;
                        if (do_exactsum())
                                comp += (tmp - res) - y;
                        res = tmp;
                }
                if (do_exactsum())
                        res -= comp;
        }
#ifdef SMPVEC_DEBUG
        STD__ cerr << "fabssqr: " << STD__ setprecision (16) << res << STD__ endl;
#endif
        return res;
}

#endif /* SMP */

INST(template <typename T> class TVector friend double fabssqr (TVector<T>);)
template <typename T>
inline double fabssqr (TVector<T> tv)
{ return tv.fabssqr (); }


/*
template <typename T>
Vector<T>::operator Vector<double> () const
{ 
  Vector<double> v(dim); for (unsigned long i = 0; i < dim; ++i)
    v.vec[i] = CPLX__ real(vec[i]); 
  return v;
}
*/

// Here are the emul, cemul, ediv and cediv variations:
#ifdef INC_VEC_EXTRA
# include "vector_extra.h"
#endif

NAMESPACE_END

// ---------------------------------------------------------
// Specializations for int types
//----------------------------------------------------------

#ifdef VEC_INT_DIV_SUPPORT

NAMESPACE_TBCI
   
//INST(template <> class TSVector friend TSVector<int> operator / (const int&, const TVector<int>&);)
template <>
inline TSVector<int> operator / (const int& a, TVector<int> b)
{
  double dv = (double)a / b.fabssqr ();
  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
  for (unsigned long i = 0; i < b.size (); ++i)
     b.set((int)(b.get(i) * dv), i);
  return TSVector<int> (b, (int)1);
}

template <>
inline TVector<int>& TVector<int>::operator /= (const int& a)
{
  BCHK(a==0, VecErr, Divide TVector<int> by zero, 0, *this);
  for (REGISTER unsigned long i = 0; i < this->dim; ++i) 
        this->vec[i] /= a;
  return *this;
}

template <> /* Was commented out. Why ??? */
inline TSVector<int> operator / (const int& a, const Vector<int>& b)
{
  double dv = (double)a / b.fabssqr ();
  TVector<int> r (b.size ());
  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, r);
  for (unsigned long i = 0; i < b.size (); ++i)
     r.set ((int)(dv * b.get(i)), i);
  return TSVector<int>(r, 1);
}

template <>
inline TSVector<int> Vector<int>::operator / (const int& a) const
{
  TVector<int> r (this->dim);
  BCHK(a==0, VecErr, Divide Vector<int> by 0, 0, r);
  for (unsigned long i = 0; i < this->dim; ++i)
        r.vec[i] = this->vec[i] / a;
  return TSVector<int> (r, 1);
}

// ---------------------------------------------------------
// Specializations for unsigned types
//----------------------------------------------------------

template <>
inline TVector<unsigned>& TVector<unsigned>::operator /= (const unsigned& a)
{
  BCHK(a==0, VecErr, Divide TVector<unsigned> by zero, 0, *this);
  for (REGISTER unsigned long i = 0; i < this->dim; ++i) 
        this->vec[i] /= a;
  return *this;
}


#ifndef HAVE_WIN_32
//does not compile under VC 6.0 because it has been already used
template <>
inline TSVector<unsigned> operator / (const unsigned& a, TVector<unsigned> b)
{
  double dv = (double)a / b.fabssqr ();
  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
  for (unsigned long i = 0; i < b.size (); ++i)
     b.set((unsigned)(b.get(i) * dv), i);
  return TSVector<unsigned> (b, (unsigned)1);
}
#endif /* HAVE_WIN_32 */
   
template <>
inline TSVector<unsigned> operator / (const unsigned& a, const Vector<unsigned>& b)
{
  double dv = (double)a / b.fabssqr ();
  BCHK(!dv, VecErr, Divide by norm 0 Vector, 0, b);
  TVector<unsigned> r (b.size ());
  for (unsigned long i = 0; i < b.size (); ++i)
     r.set ((unsigned)(dv * b.get(i)), i);
  return TSVector<unsigned>(r, 1);
}

template <>
inline TSVector<unsigned> Vector<unsigned>::operator / (const unsigned& a) const
{
  TVector<unsigned> r (this->dim);
  BCHK(a==0, VecErr, Divide Vector<unsigned> by 0, 0, r);
  for (unsigned long i = 0; i < this->dim; ++i)
        r.vec[i] = this->vec[i] / a;
  return TSVector<unsigned> (r, 1);
}

NAMESPACE_END

#endif /* VEC_INT_DIV_SUPPORT */

#endif  /* TBCI_VECTOR_H */