vec.h

Go to the documentation of this file.

#ifndef VEC_H
#define VEC_H

#include <itpp/base/itassert.h>
#include <itpp/base/math/misc.h>
#include <itpp/base/copy_vector.h>
#include <itpp/base/factory.h>
#include <vector>
#include <itpp/itexports.h>

namespace itpp
{

// Declaration of Vec
template<class Num_T> class Vec;
// Declaration of Mat
template<class Num_T> class Mat;
// Declaration of bin
class bin;

//-----------------------------------------------------------------------------------
// Declaration of Vec Friends
//-----------------------------------------------------------------------------------

template<class Num_T>
Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
template<class Num_T>
Vec<Num_T> operator+(const Vec<Num_T> &v, Num_T t);
template<class Num_T>
Vec<Num_T> operator+(Num_T t, const Vec<Num_T> &v);

template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v, Num_T t);
template<class Num_T>
Vec<Num_T> operator-(Num_T t, const Vec<Num_T> &v);
template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v);

template<class Num_T>
Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
template<class Num_T>
Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
template<class Num_T>
Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                         bool hermitian = false);
template<class Num_T>
Vec<Num_T> operator*(const Vec<Num_T> &v, Num_T t);
template<class Num_T>
Vec<Num_T> operator*(Num_T t, const Vec<Num_T> &v);

template<class Num_T>
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b);
template<class Num_T>
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b,
                     const Vec<Num_T> &c);
template<class Num_T>
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b,
                     const Vec<Num_T> &c, const Vec<Num_T> &d);

template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b,
                   Vec<Num_T> &out);
template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b,
                   const Vec<Num_T> &c, Vec<Num_T> &out);
template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b,
                   const Vec<Num_T> &c, const Vec<Num_T> &d,
                   Vec<Num_T> &out);

template<class Num_T>
void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b);
template<class Num_T>
Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);

template<class Num_T>
Vec<Num_T> operator/(const Vec<Num_T> &v, Num_T t);
template<class Num_T>
Vec<Num_T> operator/(Num_T t, const Vec<Num_T> &v);

template<class Num_T>
Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b);
template<class Num_T>
Vec<Num_T> elem_div(Num_T t, const Vec<Num_T> &v);
template<class Num_T>
void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
template<class Num_T>
Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v, Num_T a);
template<class Num_T>
Vec<Num_T> concat(Num_T a, const Vec<Num_T> &v);
template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3);
template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3, const Vec<Num_T> &v4);
template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3, const Vec<Num_T> &v4,
                  const Vec<Num_T> &v5);

//-----------------------------------------------------------------------------------
// Declaration of Vec
//-----------------------------------------------------------------------------------

template<class Num_T>
class Vec
{
public:
  typedef Num_T value_type;

  explicit Vec(const Factory &f = DEFAULT_FACTORY);
  explicit Vec(int size, const Factory &f = DEFAULT_FACTORY);
  Vec(const Vec<Num_T> &v);
  Vec(const Vec<Num_T> &v, const Factory &f);
  Vec(const char *str, const Factory &f = DEFAULT_FACTORY);
  Vec(const std::string &str, const Factory &f = DEFAULT_FACTORY);
  Vec(const Num_T *c_array, int size, const Factory &f = DEFAULT_FACTORY);

  ~Vec();

  int length() const { return datasize; }
  int size() const { return datasize; }

  void set_size(int size, bool copy = false);
  void set_length(int size, bool copy = false) { set_size(size, copy); }
  void zeros();
  void clear() { zeros(); }
  void ones();
  void set(const char *str);
  void set(const std::string &str);

  const Num_T &operator[](int i) const;
  const Num_T &operator()(int i) const;
  Num_T &operator[](int i);
  Num_T &operator()(int i);
  Vec<Num_T> operator()(int i1, int i2) const;
  Vec<Num_T> operator()(const Vec<int> &indexlist) const;
  Vec<Num_T> operator()(const Vec<bin> &binlist) const;

  const Num_T &get(int i) const;
  Vec<Num_T> get(int i1, int i2) const;
  Vec<Num_T> get(const Vec<int> &indexlist) const;
  Vec<Num_T> get(const Vec<bin> &binlist) const;

  void set(int i, Num_T t);

  Mat<Num_T> transpose() const;
  Mat<Num_T> T() const { return this->transpose(); }
  Mat<Num_T> hermitian_transpose() const;
  Mat<Num_T> H() const { return this->hermitian_transpose(); }

  Vec<Num_T>& operator+=(const Vec<Num_T> &v);
  Vec<Num_T>& operator+=(Num_T t);
  friend Vec<Num_T> operator+<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
  friend Vec<Num_T> operator+<>(const Vec<Num_T> &v, Num_T t);
  friend Vec<Num_T> operator+<>(Num_T t, const Vec<Num_T> &v);

  Vec<Num_T>& operator-=(const Vec<Num_T> &v);
  Vec<Num_T>& operator-=(Num_T t);
  friend Vec<Num_T> operator-<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
  friend Vec<Num_T> operator-<>(const Vec<Num_T> &v, Num_T t);
  friend Vec<Num_T> operator-<>(Num_T t, const Vec<Num_T> &v);
  friend Vec<Num_T> operator-<>(const Vec<Num_T> &v);

  Vec<Num_T>& operator*=(Num_T t);
  friend Vec<Num_T> operator*<>(const Vec<Num_T> &v, Num_T t);
  friend Vec<Num_T> operator*<>(Num_T t, const Vec<Num_T> &v);

  friend Vec<Num_T> elem_mult<>(const Vec<Num_T> &a, const Vec<Num_T> &b);
  friend Vec<Num_T> elem_mult<>(const Vec<Num_T> &a, const Vec<Num_T> &b,
                                const Vec<Num_T> &c);
  friend Vec<Num_T> elem_mult<>(const Vec<Num_T> &a, const Vec<Num_T> &b,
                                const Vec<Num_T> &c, const Vec<Num_T> &d);

  friend void elem_mult_out<>(const Vec<Num_T> &a, const Vec<Num_T> &b,
                              Vec<Num_T> &out);
  friend void elem_mult_out<>(const Vec<Num_T> &a, const Vec<Num_T> &b,
                              const Vec<Num_T> &c, Vec<Num_T> &out);
  friend void elem_mult_out<>(const Vec<Num_T> &a, const Vec<Num_T> &b,
                              const Vec<Num_T> &c, const Vec<Num_T> &d,
                              Vec<Num_T> &out);

  friend void elem_mult_inplace<>(const Vec<Num_T> &a, Vec<Num_T> &b);
  friend Num_T elem_mult_sum<>(const Vec<Num_T> &a, const Vec<Num_T> &b);

  Vec<Num_T>& operator/=(Num_T t);
  Vec<Num_T>& operator/=(const Vec<Num_T> &v);

  friend Vec<Num_T> operator/<>(const Vec<Num_T> &v, Num_T t);
  friend Vec<Num_T> operator/<>(Num_T t, const Vec<Num_T> &v);

  friend Vec<Num_T> elem_div<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
  friend Vec<Num_T> elem_div<>(Num_T t, const Vec<Num_T> &v);
  friend void elem_div_out<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                             Vec<Num_T> &out);
  friend Num_T elem_div_sum<>(const Vec<Num_T> &a, const Vec<Num_T> &b);

  Vec<Num_T> right(int nr) const;
  Vec<Num_T> left(int nr) const;
  Vec<Num_T> mid(int start, int nr) const;
  Vec<Num_T> split(int pos);
  void shift_right(Num_T t, int n = 1);
  void shift_right(const Vec<Num_T> &v);
  void shift_left(Num_T t, int n = 1);
  void shift_left(const Vec<Num_T> &v);

  friend Vec<Num_T> concat<>(const Vec<Num_T> &v, Num_T t);
  friend Vec<Num_T> concat<>(Num_T t, const Vec<Num_T> &v);
  friend Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
  friend Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                             const Vec<Num_T> &v3);
  friend Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                             const Vec<Num_T> &v3, const Vec<Num_T> &v4);
  friend Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                             const Vec<Num_T> &v3, const Vec<Num_T> &v4,
                             const Vec<Num_T> &v5);

  void set_subvector(int i1, int i2, const Vec<Num_T> &v);
  void set_subvector(int i, const Vec<Num_T> &v);
  void set_subvector(int i1, int i2, Num_T t);
  void replace_mid(int i, const Vec<Num_T> &v);
  void del(int i);
  void del(int i1, int i2);
  void ins(int i, Num_T t);
  void ins(int i, const Vec<Num_T> &v);

  Vec<Num_T>& operator=(Num_T t);
  Vec<Num_T>& operator=(const Vec<Num_T> &v);
  Vec<Num_T>& operator=(const Mat<Num_T> &m);
  Vec<Num_T>& operator=(const char *str);
  Vec<Num_T>& operator=(const std::string &str);

  Vec<bin> operator==(Num_T t) const;
  Vec<bin> operator!=(Num_T t) const;
  Vec<bin> operator<(Num_T t) const;
  Vec<bin> operator<=(Num_T t) const;
  Vec<bin> operator>(Num_T t) const;
  Vec<bin> operator>=(Num_T t) const;

  bool operator==(const Vec<Num_T> &v) const;
  bool operator!=(const Vec<Num_T> &v) const;

  Num_T &_elem(int i) { return data[i]; }
  const Num_T &_elem(int i) const { return data[i]; }

  Num_T *_data() { return data; }
  const Num_T *_data() const { return data; }

protected:
  void alloc(int size);
  void free();

  int datasize;
  Num_T *data;
  const Factory &factory;
private:
  // Clean up and tokenize input initialisation string
  std::vector<std::string> tokenize(const std::string &str,
                                    bool &abc_format) const;
  // Parse double and integer values from string tokens
  Num_T parse_token(const std::string &s) const;
  // Parse \c a, \c b and \c c values from "a:b:c" format
  void parse_abc_token(const std::string &s, Num_T &a, Num_T &b,
                       Num_T &c) const;
  bool in_range(int i) const { return ((i < datasize) && (i >= 0)); }
};

//-----------------------------------------------------------------------------------
// Type definitions of vec, cvec, ivec, svec, and bvec
//-----------------------------------------------------------------------------------

typedef Vec<double> vec;

typedef Vec<std::complex<double> > cvec;

typedef Vec<int> ivec;

typedef Vec<short int> svec;

typedef Vec<bin> bvec;

} //namespace itpp


#include <itpp/base/mat.h>

namespace itpp
{

//-----------------------------------------------------------------------------------
// Declaration of input and output streams for Vec
//-----------------------------------------------------------------------------------

template<class Num_T>
std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v);

template<class Num_T>
std::istream &operator>>(std::istream &is, Vec<Num_T> &v);

//-----------------------------------------------------------------------------------
// Implementation of templated Vec members and friends
//-----------------------------------------------------------------------------------

template<class Num_T> inline
void Vec<Num_T>::alloc(int size)
{
  if (size > 0) {
    create_elements(data, size, factory);
    datasize = size;
  }
  else {
    data = 0;
    datasize = 0;
  }
}

template<class Num_T> inline
void Vec<Num_T>::free()
{
  destroy_elements(data, datasize);
  datasize = 0;
}


template<class Num_T> inline
Vec<Num_T>::Vec(const Factory &f) : datasize(0), data(0), factory(f) {}

template<class Num_T> inline
Vec<Num_T>::Vec(int size, const Factory &f) : datasize(0), data(0), factory(f)
{
  it_assert_debug(size >= 0, "Negative size in Vec::Vec(int)");
  alloc(size);
}

template<class Num_T> inline
Vec<Num_T>::Vec(const Vec<Num_T> &v) : datasize(0), data(0), factory(v.factory)
{
  alloc(v.datasize);
  copy_vector(datasize, v.data, data);
}

template<class Num_T> inline
Vec<Num_T>::Vec(const Vec<Num_T> &v, const Factory &f) : datasize(0), data(0), factory(f)
{
  alloc(v.datasize);
  copy_vector(datasize, v.data, data);
}

template<class Num_T> inline
Vec<Num_T>::Vec(const char *str, const Factory &f) : datasize(0), data(0), factory(f)
{
  set(std::string(str));
}

template<class Num_T> inline
Vec<Num_T>::Vec(const std::string &str, const Factory &f) : datasize(0), data(0), factory(f)
{
  set(str);
}

template<class Num_T> inline
Vec<Num_T>::Vec(const Num_T *c_array, int size, const Factory &f) : datasize(0), data(0), factory(f)
{
  alloc(size);
  copy_vector(size, c_array, data);
}

template<class Num_T> inline
Vec<Num_T>::~Vec()
{
  free();
}

template<class Num_T>
void Vec<Num_T>::set_size(int size, bool copy)
{
  it_assert_debug(size >= 0, "Vec::set_size(): New size must not be negative");
  if (datasize == size)
    return;
  if (copy) {
    // create a temporary pointer to the allocated data
    Num_T* tmp = data;
    // store the current number of elements
    int old_datasize = datasize;
    // check how many elements we need to copy
    int min = datasize < size ? datasize : size;
    // allocate new memory
    alloc(size);
    // copy old elements into a new memory region
    copy_vector(min, tmp, data);
    // initialize the rest of resized vector
    for (int i = min; i < size; ++i)
      data[i] = Num_T(0);
    // delete old elements
    destroy_elements(tmp, old_datasize);
  }
  else {
    free();
    alloc(size);
  }
}

template<class Num_T> inline
const Num_T& Vec<Num_T>::operator[](int i) const
{
  it_assert_debug(in_range(i), "Vec<>::operator[]: Index out of range");
  return data[i];
}

template<class Num_T> inline
const Num_T& Vec<Num_T>::operator()(int i) const
{
  return (*this)[i];
}

template<class Num_T> inline
Num_T& Vec<Num_T>::operator[](int i)
{
  it_assert_debug(in_range(i), "Vec<>::operator[]: Index out of range");
  return data[i];
}

template<class Num_T> inline
Num_T& Vec<Num_T>::operator()(int i)
{
  return (*this)[i];
}

template<class Num_T> inline
Vec<Num_T> Vec<Num_T>::operator()(int i1, int i2) const
{
  if (i1 == -1) i1 = datasize - 1;
  if (i2 == -1) i2 = datasize - 1;

  it_assert_debug((i1 >= 0) && (i1 <= i2) && (i2 < datasize),
                  "Vec<>::operator()(i1, i2): Indexing out of range");

  Vec<Num_T> s(i2 - i1 + 1);
  copy_vector(s.datasize, data + i1, s.data);

  return s;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::operator()(const Vec<int> &indexlist) const
{
  int size = indexlist.size();
  Vec<Num_T> temp(size);
  for (int i = 0; i < size; ++i) {
    it_assert_debug(in_range(indexlist(i)), "Vec<>::operator()(ivec &): "
                    "Index i=" << i << " out of range");
    temp(i) = data[indexlist(i)];
  }
  return temp;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::operator()(const Vec<bin> &binlist) const
{
  int size = binlist.size();
  it_assert_debug(datasize == size, "Vec<>::operator()(bvec &): "
                  "Wrong size of binlist vector");
  Vec<Num_T> temp(size);
  int j = 0;
  for (int i = 0; i < size; ++i)
    if (binlist(i) == bin(1))
      temp(j++) = data[i];
  temp.set_size(j, true);
  return temp;
}


template<class Num_T> inline
const Num_T& Vec<Num_T>::get(int i) const
{
  return (*this)[i];
}

template<class Num_T> inline
Vec<Num_T> Vec<Num_T>::get(int i1, int i2) const
{
  return (*this)(i1, i2);
}

template<class Num_T> inline
Vec<Num_T> Vec<Num_T>::get(const Vec<int> &indexlist) const
{
  return (*this)(indexlist);
}

template<class Num_T> inline
Vec<Num_T> Vec<Num_T>::get(const Vec<bin> &binlist) const
{
  return (*this)(binlist);
}

template<class Num_T> inline
void Vec<Num_T>::zeros()
{
  for (int i = 0; i < datasize; i++)
    data[i] = Num_T(0);
}

template<class Num_T> inline
void Vec<Num_T>::ones()
{
  for (int i = 0; i < datasize; i++)
    data[i] = Num_T(1);
}

template<class Num_T> inline
void Vec<Num_T>::set(int i, Num_T t)
{
  it_assert_debug(in_range(i), "Vec<>::set(i, t): Index out of range");
  data[i] = t;
}

template<class Num_T> inline
void Vec<Num_T>::set(const std::string &str)
{
  it_error("Vec::set(): Only `double', `complex<double>', `int', "
           "`short int' and `bin' types supported");
}

template<class Num_T> inline
void Vec<Num_T>::set(const char *str)
{
  set(std::string(str));
}

template<>
ITPP_EXPORT void Vec<double>::set(const std::string &str);
template<>
ITPP_EXPORT void Vec<std::complex<double> >::set(const std::string &str);
template<>
ITPP_EXPORT void Vec<int>::set(const std::string &str);
template<>
ITPP_EXPORT void Vec<short int>::set(const std::string &str);
template<>
ITPP_EXPORT void Vec<bin>::set(const std::string &str);

template<class Num_T>
Mat<Num_T> Vec<Num_T>::transpose() const
{
  Mat<Num_T> temp(1, datasize);
  copy_vector(datasize, data, temp._data());
  return temp;
}

template<class Num_T>
Mat<Num_T> Vec<Num_T>::hermitian_transpose() const
{
  Mat<Num_T> temp(1, datasize);
  copy_vector(datasize, data, temp._data());
  return temp;
}

template<>
ITPP_EXPORT Mat<std::complex<double> > Vec<std::complex<double> >::hermitian_transpose() const;


template<class Num_T>
Vec<Num_T>& Vec<Num_T>::operator+=(const Vec<Num_T> &v)
{
  if (datasize == 0) { // if not assigned a size.
    if (this != &v) { // check for self addition
      alloc(v.datasize);
      copy_vector(datasize, v.data, data);
    }
  }
  else {
    it_assert_debug(datasize == v.datasize, "Vec::operator+=: Wrong sizes");
    for (int i = 0; i < datasize; i++)
      data[i] += v.data[i];
  }
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator+=(Num_T t)
{
  for (int i = 0;i < datasize;i++)
    data[i] += t;
  return *this;
}

template<class Num_T>
Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
{
  int i;
  Vec<Num_T> r(v1.datasize);

  it_assert_debug(v1.datasize == v2.datasize, "Vec::operator+: wrong sizes");
  for (i = 0; i < v1.datasize; i++)
    r.data[i] = v1.data[i] + v2.data[i];

  return r;
}

template<class Num_T>
Vec<Num_T> operator+(const Vec<Num_T> &v, Num_T t)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = v.data[i] + t;

  return r;
}

template<class Num_T>
Vec<Num_T> operator+(Num_T t, const Vec<Num_T> &v)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = t + v.data[i];

  return r;
}

template<class Num_T>
Vec<Num_T>& Vec<Num_T>::operator-=(const Vec<Num_T> &v)
{
  if (datasize == 0) { // if not assigned a size.
    if (this != &v) { // check for self decrementation
      alloc(v.datasize);
      for (int i = 0; i < v.datasize; i++)
        data[i] = -v.data[i];
    }
  }
  else {
    it_assert_debug(datasize == v.datasize, "Vec::operator-=: Wrong sizes");
    for (int i = 0; i < datasize; i++)
      data[i] -= v.data[i];
  }
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator-=(Num_T t)
{
  for (int i = 0;i < datasize;i++)
    data[i] -= t;
  return *this;
}

template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
{
  int i;
  Vec<Num_T> r(v1.datasize);

  it_assert_debug(v1.datasize == v2.datasize, "Vec::operator-: wrong sizes");
  for (i = 0; i < v1.datasize; i++)
    r.data[i] = v1.data[i] - v2.data[i];

  return r;
}

template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v, Num_T t)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = v.data[i] - t;

  return r;
}

template<class Num_T>
Vec<Num_T> operator-(Num_T t, const Vec<Num_T> &v)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = t - v.data[i];

  return r;
}

template<class Num_T>
Vec<Num_T> operator-(const Vec<Num_T> &v)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = -v.data[i];

  return r;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator*=(Num_T t)
{
  scal_vector(datasize, t, data);
  return *this;
}

template<class Num_T> inline
Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
{
  return dot(v1, v2);
}

template<class Num_T>
Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
{
  it_assert_debug(v1.length() == v2.length(), "Vec::dot(): Wrong sizes");
  Num_T r = Num_T(0);
  for (int i = 0; i < v1.length(); ++i)
    r += v1._data()[i] * v2._data()[i];
  return r;
}

template<>
ITPP_EXPORT double dot(const vec &v1, const vec &v2);


template<class Num_T>
Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2, bool)
{
  it_assert_debug((v1.length() > 0) && (v2.length() > 0),
                  "Vec::outer_product:: Input vector of zero size");

  Mat<Num_T> r(v1.length(), v2.length());
  for (int i = 0; i < v1.length(); ++i) {
    for (int j = 0; j < v2.length(); ++j) {
      r(i, j) = v1._data()[i] * v2._data()[j];
    }
  }
  return r;
}

template<>
ITPP_EXPORT mat outer_product(const vec &v1, const vec &v2, bool);

template<>
ITPP_EXPORT cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian);

template<class Num_T>
Vec<Num_T> operator*(const Vec<Num_T> &v, Num_T t)
{
  int i;
  Vec<Num_T> r(v.datasize);
  for (i = 0; i < v.datasize; i++)
    r.data[i] = v.data[i] * t;

  return r;
}

template<class Num_T> inline
Vec<Num_T> operator*(Num_T t, const Vec<Num_T> &v)
{
  return operator*(v, t);
}

template<class Num_T> inline
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b)
{
  Vec<Num_T> out;
  elem_mult_out(a, b, out);
  return out;
}

template<class Num_T> inline
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b,
                     const Vec<Num_T> &c)
{
  Vec<Num_T> out;
  elem_mult_out(a, b, c, out);
  return out;
}

template<class Num_T> inline
Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b,
                     const Vec<Num_T> &c, const Vec<Num_T> &d)
{
  Vec<Num_T> out;
  elem_mult_out(a, b, c, d, out);
  return out;
}

template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
{
  it_assert_debug(a.datasize == b.datasize,
                  "Vec<>::elem_mult_out(): Wrong sizes");
  out.set_size(a.datasize);
  for (int i = 0; i < a.datasize; i++)
    out.data[i] = a.data[i] * b.data[i];
}

template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b,
                   const Vec<Num_T> &c, Vec<Num_T> &out)
{
  it_assert_debug((a.datasize == b.datasize) && (a.datasize == c.datasize),
                  "Vec<>::elem_mult_out(): Wrong sizes");
  out.set_size(a.datasize);
  for (int i = 0; i < a.datasize; i++)
    out.data[i] = a.data[i] * b.data[i] * c.data[i];
}

template<class Num_T>
void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b,
                   const Vec<Num_T> &c, const Vec<Num_T> &d, Vec<Num_T> &out)
{
  it_assert_debug((a.datasize == b.datasize) && (a.datasize == c.datasize)
                  && (a.datasize == d.datasize),
                  "Vec<>::elem_mult_out(): Wrong sizes");
  out.set_size(a.datasize);
  for (int i = 0; i < a.datasize; i++)
    out.data[i] = a.data[i] * b.data[i] * c.data[i] * d.data[i];
}

template<class Num_T>
#ifndef _MSC_VER
inline
#endif
void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b)
{
  it_assert_debug(a.datasize == b.datasize,
                  "Vec<>::elem_mult_inplace(): Wrong sizes");
  for (int i = 0; i < a.datasize; i++)
    b.data[i] *= a.data[i];
}

template<class Num_T> inline
Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
{
  it_assert_debug(a.datasize == b.datasize,
                  "Vec<>::elem_mult_sum(): Wrong sizes");
  Num_T acc = 0;
  for (int i = 0; i < a.datasize; i++)
    acc += a.data[i] * b.data[i];
  return acc;
}

template<class Num_T>
Vec<Num_T> operator/(const Vec<Num_T> &v, Num_T t)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = v.data[i] / t;

  return r;
}

template<class Num_T>
Vec<Num_T> operator/(Num_T t, const Vec<Num_T> &v)
{
  int i;
  Vec<Num_T> r(v.datasize);

  for (i = 0; i < v.datasize; i++)
    r.data[i] = t / v.data[i];

  return r;
}

template<class Num_T>
Vec<Num_T> elem_div(Num_T t, const Vec<Num_T> &v)
{
  it_warning("Vec<>::elem_div(Num_T, const Vec<Num_T> &): This function is "
             "deprecated and might be removed from future IT++ releases. "
             "Please use Vec<>::operator/(Num_T, const Vec<Num_T> &) "
             "instead.");
  return operator/(t, v);
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator/=(Num_T t)
{
  for (int i = 0; i < datasize; ++i) {
    data[i] /= t;
  }
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator/=(const Vec<Num_T> &v)
{
  it_assert_debug(datasize == v.datasize, "Vec::operator/=(): wrong sizes");
  for (int i = 0; i < datasize; ++i) {
    data[i] /= v.data[i];
  }
  return *this;
}

template<class Num_T> inline
Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b)
{
  Vec<Num_T> out;
  elem_div_out(a, b, out);
  return out;
}

template<class Num_T>
void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
{
  it_assert_debug(a.datasize == b.datasize, "Vecelem_div_out: wrong sizes");

  out.set_size(a.size());

  for (int i = 0; i < a.datasize; i++)
    out.data[i] = a.data[i] / b.data[i];
}

template<class Num_T> inline
Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
{
  it_assert_debug(a.datasize == b.datasize, "Vec::elem_div_sum: wrong sizes");

  Num_T acc = 0;

  for (int i = 0; i < a.datasize; i++)
    acc += a.data[i] / b.data[i];

  return acc;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::right(int nr) const
{
  it_assert_debug(nr <= datasize, "Vec::right(): index out of range");
  Vec<Num_T> temp(nr);
  if (nr > 0) {
    copy_vector(nr, &data[datasize-nr], temp.data);
  }
  return temp;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::left(int nr) const
{
  it_assert_debug(nr <= datasize, "Vec::left(): index out of range");
  Vec<Num_T> temp(nr);
  if (nr > 0) {
    copy_vector(nr, data, temp.data);
  }
  return temp;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::mid(int start, int nr) const
{
  it_assert_debug((start >= 0) && ((start + nr) <= datasize),
                  "Vec::mid(): indexing out of range");
  Vec<Num_T> temp(nr);
  if (nr > 0) {
    copy_vector(nr, &data[start], temp.data);
  }
  return temp;
}

template<class Num_T>
Vec<Num_T> Vec<Num_T>::split(int pos)
{
  it_assert_debug((pos >= 0) && (pos <= datasize),
                  "Vec<>::split(): Index out of range");
  Vec<Num_T> temp1(pos);
  if (pos > 0) {
    copy_vector(pos, data, temp1.data);
    if (pos < datasize) {
      Vec<Num_T> temp2(datasize - pos);
      copy_vector(datasize - pos, &data[pos], temp2.data);
      (*this) = temp2;
    }
    else {
      set_size(0);
    }
  }
  return temp1;
}

template<class Num_T>
void Vec<Num_T>::shift_right(Num_T t, int n)
{
  int i = datasize;

  it_assert_debug(n >= 0, "Vec::shift_right: index out of range");
  while (--i >= n)
    data[i] = data[i-n];
  while (i >= 0)
    data[i--] = t;
}

template<class Num_T>
void Vec<Num_T>::shift_right(const Vec<Num_T> &v)
{
  for (int i = datasize - 1; i >= v.datasize; i--)
    data[i] = data[i-v.datasize];
  for (int i = 0; i < v.datasize; i++)
    data[i] = v[i];
}

template<class Num_T>
void Vec<Num_T>::shift_left(Num_T t, int n)
{
  int i;

  it_assert_debug(n >= 0, "Vec::shift_left: index out of range");
  for (i = 0; i < datasize - n; i++)
    data[i] = data[i+n];
  while (i < datasize)
    data[i++] = t;
}

template<class Num_T>
void Vec<Num_T>::shift_left(const Vec<Num_T> &v)
{
  for (int i = 0; i < datasize - v.datasize; i++)
    data[i] = data[i+v.datasize];
  for (int i = datasize - v.datasize; i < datasize; i++)
    data[i] = v[i-datasize+v.datasize];
}

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v, Num_T t)
{
  int size = v.size();
  Vec<Num_T> temp(size + 1);
  copy_vector(size, v.data, temp.data);
  temp(size) = t;
  return temp;
}

template<class Num_T>
Vec<Num_T> concat(Num_T t, const Vec<Num_T> &v)
{
  int size = v.size();
  Vec<Num_T> temp(size + 1);
  temp(0) = t;
  copy_vector(size, v.data, &temp.data[1]);
  return temp;
}

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
{
  int size1 = v1.size();
  int size2 = v2.size();
  Vec<Num_T> temp(size1 + size2);
  copy_vector(size1, v1.data, temp.data);
  copy_vector(size2, v2.data, &temp.data[size1]);
  return temp;
}

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3)
{
  int size1 = v1.size();
  int size2 = v2.size();
  int size3 = v3.size();
  Vec<Num_T> temp(size1 + size2 + size3);
  copy_vector(size1, v1.data, temp.data);
  copy_vector(size2, v2.data, &temp.data[size1]);
  copy_vector(size3, v3.data, &temp.data[size1+size2]);
  return temp;
}

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3, const Vec<Num_T> &v4)
{
  int size1 = v1.size();
  int size2 = v2.size();
  int size3 = v3.size();
  int size4 = v4.size();
  Vec<Num_T> temp(size1 + size2 + size3 + size4);
  copy_vector(size1, v1.data, temp.data);
  copy_vector(size2, v2.data, &temp.data[size1]);
  copy_vector(size3, v3.data, &temp.data[size1+size2]);
  copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
  return temp;
}

template<class Num_T>
Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
                  const Vec<Num_T> &v3, const Vec<Num_T> &v4,
                  const Vec<Num_T> &v5)
{
  int size1 = v1.size();
  int size2 = v2.size();
  int size3 = v3.size();
  int size4 = v4.size();
  int size5 = v5.size();
  Vec<Num_T> temp(size1 + size2 + size3 + size4 + size5);
  copy_vector(size1, v1.data, temp.data);
  copy_vector(size2, v2.data, &temp.data[size1]);
  copy_vector(size3, v3.data, &temp.data[size1+size2]);
  copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
  copy_vector(size5, v5.data, &temp.data[size1+size2+size3+size4]);
  return temp;
}

template<class Num_T>
void Vec<Num_T>::set_subvector(int i1, int, const Vec<Num_T> &v)
{
  it_warning("Vec<>::set_subvector(int, int, const Vec<> &): This function "
             "is deprecated and might be removed from future IT++ releases. "
             "Please use Vec<>::set_subvector(int, const Vec<> &) instead.");
  set_subvector(i1, v);
}

template<class Num_T> inline
void Vec<Num_T>:: set_subvector(int i, const Vec<Num_T> &v)
{
  it_assert_debug((i >= 0) && (i + v.datasize <= datasize),
                  "Vec<>::set_subvector(int, const Vec<> &): "
                  "Index out of range or too long input vector");
  copy_vector(v.datasize, v.data, data + i);
}

template<class Num_T> inline
void Vec<Num_T>::set_subvector(int i1, int i2, Num_T t)
{
  if (i1 == -1) i1 = datasize - 1;
  if (i2 == -1) i2 = datasize - 1;
  it_assert_debug((i1 >= 0) && (i1 <= i2) && (i2 < datasize),
                  "Vec<>::set_subvector(int, int, Num_T): Indexing out "
                  "of range");
  for (int i = i1; i <= i2; i++)
    data[i] = t;
}

template<class Num_T> inline
void Vec<Num_T>::replace_mid(int i, const Vec<Num_T> &v)
{
  set_subvector(i, v);
}

template<class Num_T>
void Vec<Num_T>::del(int index)
{
  it_assert_debug(in_range(index), "Vec<>::del(int): Index out of range");
  Vec<Num_T> temp(*this);
  set_size(datasize - 1, false);
  copy_vector(index, temp.data, data);
  copy_vector(datasize - index, &temp.data[index+1], &data[index]);
}

template<class Num_T>
void Vec<Num_T>::del(int i1, int i2)
{
  if (i1 == -1) i1 = datasize - 1;
  if (i2 == -1) i2 = datasize - 1;
  it_assert_debug((i1 >= 0) && (i1 <= i2) && (i2 < datasize),
                  "Vec<>::del(int, int): Indexing out of range");
  Vec<Num_T> temp(*this);
  int new_size = datasize - (i2 - i1 + 1);
  set_size(new_size, false);
  copy_vector(i1, temp.data, data);
  copy_vector(datasize - i1, &temp.data[i2+1], &data[i1]);
}

template<class Num_T>
void Vec<Num_T>::ins(int index, const Num_T t)
{
  it_assert_debug((index >= 0) && (index <= datasize),
                  "Vec<>::ins(): Index out of range");
  Vec<Num_T> Temp(*this);

  set_size(datasize + 1, false);
  copy_vector(index, Temp.data, data);
  data[index] = t;
  copy_vector(Temp.datasize - index, Temp.data + index, data + index + 1);
}

template<class Num_T>
void Vec<Num_T>::ins(int index, const Vec<Num_T> &v)
{
  it_assert_debug((index >= 0) && (index <= datasize),
                  "Vec<>::ins(): Index out of range");
  Vec<Num_T> Temp(*this);

  set_size(datasize + v.length(), false);
  copy_vector(index, Temp.data, data);
  copy_vector(v.size(), v.data, &data[index]);
  copy_vector(Temp.datasize - index, Temp.data + index, data + index + v.size());
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator=(Num_T t)
{
  for (int i = 0;i < datasize;i++)
    data[i] = t;
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator=(const Vec<Num_T> &v)
{
  if (this != &v) {
    set_size(v.datasize, false);
    copy_vector(datasize, v.data, data);
  }
  return *this;
}

template<class Num_T>
Vec<Num_T>& Vec<Num_T>::operator=(const Mat<Num_T> &m)
{
  if (m.cols() == 1) {
    set_size(m.rows(), false);
    copy_vector(m.rows(), m._data(), data);
  }
  else if (m.rows() == 1) {
    set_size(m.cols(), false);
    copy_vector(m.cols(), m._data(), m.rows(), data, 1);
  }
  else
    it_error("Vec<>::operator=(Mat<Num_T> &): Wrong size of input matrix");
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator=(const char *str)
{
  set(std::string(str));
  return *this;
}

template<class Num_T> inline
Vec<Num_T>& Vec<Num_T>::operator=(const std::string &str)
{
  set(str);
  return *this;
}

template<class Num_T>
bvec Vec<Num_T>::operator==(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator==(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] == t);
  return temp;
}

template<class Num_T>
bvec Vec<Num_T>::operator!=(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator!=(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] != t);
  return temp;
}

template<>
bvec Vec<std::complex<double> >::operator<(std::complex<double>) const;

template<class Num_T>
bvec Vec<Num_T>::operator<(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator<(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] < t);
  return temp;
}

template<class Num_T>
bvec Vec<Num_T>::operator<=(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator<=(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] <= t);
  return temp;
}

template<>
ITPP_EXPORT bvec Vec<std::complex<double> >::operator<=(std::complex<double>) const;

template<class Num_T>
bvec Vec<Num_T>::operator>(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator>(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] > t);
  return temp;
}

template<>
ITPP_EXPORT bvec Vec<std::complex<double> >::operator>(std::complex<double>) const;

template<class Num_T>
bvec Vec<Num_T>::operator>=(Num_T t) const
{
  it_assert_debug(datasize > 0, "Vec<>::operator>=(): Wrong size");
  bvec temp(datasize);
  for (int i = 0; i < datasize; i++)
    temp(i) = (data[i] >= t);
  return temp;
}

template<>
ITPP_EXPORT bvec Vec<std::complex<double> >::operator>=(std::complex<double>) const;

template<class Num_T>
bool Vec<Num_T>::operator==(const Vec<Num_T> &invector) const
{
  // OBS ! if wrong size, return false
  if (datasize != invector.datasize) return false;
  for (int i = 0;i < datasize;i++) {
    if (data[i] != invector.data[i]) return false;
  }
  return true;
}

template<class Num_T>
bool Vec<Num_T>::operator!=(const Vec<Num_T> &invector) const
{
  if (datasize != invector.datasize) return true;
  for (int i = 0;i < datasize;i++) {
    if (data[i] != invector.data[i]) return true;
  }
  return false;
}

template<class Num_T>
std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v)
{
  int i, sz = v.length();

  os << "[" ;
  for (i = 0; i < sz; i++) {
    os << v(i) ;
    if (i < sz - 1)
      os << " ";
  }
  os << "]" ;

  return os;
}

template<class Num_T>
std::istream &operator>>(std::istream &is, Vec<Num_T> &v)
{
  std::ostringstream buffer;
  bool started = false;
  bool finished = false;
  bool brackets = false;
  char c = 0;

  while (!finished) {
    if (is.eof()) {
      finished = true;
    }
    else {
      is.get(c);

      if (is.eof() || (c == '\n')) {
        if (brackets) {
          // Right bracket missing
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else if (!((c == '\n') && !started)) {
          finished = true;
        }
      }
      else if ((c == ' ') || (c == '\t')) {
        if (started) {
          buffer << ' ';
        }
      }
      else if (c == '[') {
        if (started) {
          // Unexpected left bracket
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else {
          started = true;
          brackets = true;
        }
      }
      else if (c == ']') {
        if (!started || !brackets) {
          // Unexpected right bracket
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else {
          finished = true;
        }
        while (!is.eof() && (((c = static_cast<char>(is.peek())) == ' ')
                             || (c == '\t'))) {
          is.get();
        }
        if (!is.eof() && (c == '\n')) {
          is.get();
        }
      }
      else {
        started = true;
        buffer << c;
      }
    }
  }

  if (!started) {
    v.set_size(0, false);
  }
  else {
    v.set(buffer.str());
  }

  return is;
}


// ----------------------------------------------------------------------
// Private functions
// ----------------------------------------------------------------------

template<class Num_T>
void Vec<Num_T>::parse_abc_token(const std::string &s, Num_T &a, Num_T &b,
                                 Num_T &c) const
{
  std::string::size_type beg = 0;
  std::string::size_type end = s.find(':', 0);
  a = parse_token(s.substr(beg, end-beg));
  beg = end + 1;
  end = s.find(':', beg);
  if (end != std::string::npos) {
    b = parse_token(s.substr(beg, end-beg));
    c = parse_token(s.substr(end+1, s.size()-end));
  }
  else {
    b = Num_T(1);
    c = parse_token(s.substr(beg, end-beg-1));
  }
}

template<class Num_T>
Num_T Vec<Num_T>::parse_token(const std::string &s) const
{
  it_error("Vec::parse_token(): Only `double' and `int' types are supported");
  return 0;
}

template<>
ITPP_EXPORT double Vec<double>::parse_token(const std::string &s) const;
template<>
ITPP_EXPORT int Vec<int>::parse_token(const std::string &s) const;

template<class Num_T>
std::vector<std::string> Vec<Num_T>::tokenize(const std::string &str_in,
                                              bool &abc_format) const
{
  std::vector<std::string> vs;  // vector for storing parsed tokens
  std::string s;                // currently processed token string
  bool start = true;
  bool space = false;
  bool colon = false;
  bool comma = false;
  bool lparen = false;
  abc_format = false;
  for (std::string::size_type i = 0; i < str_in.size(); ++i) {
    char c = str_in[i];
    switch (c) {
    case ' ': case '\t':
      space = true;             // set flag for whitespaces
      break;
    case ',':
      if (lparen)
        comma = true;           // set flag for comma in "(re,im)" format
      else
        space = true;           // otherwise treat comma as separator
     break;
    case ')':
      s.push_back('i');         // replace right paren in "(re,im)" with 'i'
      break;
    case ':':
      colon = true;             // set flag for "a:b[:c]" format string
      space = false;            // reset flag for whitespaces
      abc_format = true;        // set external flag for "a:b[:c]" format
      s.push_back(c);
      break;
    case '(':
      lparen = true;            // set flag for complex "(re,im)" format
      break;
    default:
      if (colon) {              // reset colon and space flags
        colon = false;          // to get rid of whitespaces around ":"
        space = false;
      }
      else if (lparen && comma) { // support for "(re,im)" format
        lparen = false;
        comma = false;
        space = false;
        if ((c != '-') && (c != '+')) // if needed
          s.push_back('+');           // insert '+' between "re" and "im"
      }
      else if (space) {         // new token detected
        space = false;
        if (!start) {           // if not at the beginning of the string
          vs.push_back(s);      // store already parsed token
          s.clear();            // and start parsing the next token
        }
      }
      s.push_back(c);          // append next character to the current token
      start = false;           // reset the "beginning of the string" flag
      break;
    }
  }
  if (!s.empty())               // if the final token is not an empty string
    vs.push_back(s);            // store it in the output vector
  return vs;
}

// ----------------------------------------------------------------------
// Instantiations
// ---------------------------------------------------------------------- 
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Vec<double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Vec<int>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Vec<short int>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Vec<std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Vec<bin>;

// addition operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator+(const vec &v1, const vec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator+(const cvec &v1, const cvec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator+(const ivec &v1, const ivec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator+(const svec &v1, const svec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator+(const bvec &v1, const bvec &v2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator+(const vec &v1, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator+(const cvec &v1, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator+(const ivec &v1, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator+(const svec &v1, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator+(const bvec &v1, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator+(double t, const vec &v1);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator+(std::complex<double> t, const cvec &v1);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator+(int t, const ivec &v1);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator+(short t, const svec &v1);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator+(bin t, const bvec &v1);

// subtraction operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator-(const vec &v1, const vec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator-(const cvec &v1, const cvec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator-(const ivec &v1, const ivec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator-(const svec &v1, const svec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator-(const bvec &v1, const bvec &v2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator-(const vec &v, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator-(const cvec &v, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator-(const ivec &v, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator-(const svec &v, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator-(const bvec &v, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator-(double t, const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator-(std::complex<double> t, const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator-(int t, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator-(short t, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator-(bin t, const bvec &v);

// unary minus

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator-(const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator-(const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator-(const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator-(const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator-(const bvec &v);

// multiplication operator
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::complex<double> dot(const cvec &v1, const cvec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT int dot(const ivec &v1, const ivec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT short dot(const svec &v1, const svec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bin dot(const bvec &v1, const bvec &v2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT double operator*(const vec &v1, const vec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::complex<double> operator*(const cvec &v1, const cvec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT int operator*(const ivec &v1, const ivec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT short operator*(const svec &v1, const svec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bin operator*(const bvec &v1, const bvec &v2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT imat outer_product(const ivec &v1, const ivec &v2,
                                     bool hermitian);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT smat outer_product(const svec &v1, const svec &v2,
                                     bool hermitian);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bmat outer_product(const bvec &v1, const bvec &v2,
                                     bool hermitian);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator*(const vec &v, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator*(const cvec &v, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator*(const ivec &v, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator*(const svec &v, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator*(const bvec &v, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator*(double t, const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator*(std::complex<double> t, const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator*(int t, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator*(short t, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator*(bin t, const bvec &v);

// elementwise multiplication

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec elem_mult(const vec &a, const vec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec elem_mult(const cvec &a, const cvec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec elem_mult(const ivec &a, const ivec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec elem_mult(const svec &a, const svec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec elem_mult(const bvec &a, const bvec &b);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const vec &a, const vec &b, vec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const cvec &a, const cvec &b, cvec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const ivec &a, const ivec &b, ivec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const svec &a, const svec &b, svec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const bvec &a, const bvec &b, bvec &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec elem_mult(const vec &a, const vec &b, const vec &c);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec elem_mult(const cvec &a, const cvec &b, const cvec &c);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec elem_mult(const ivec &a, const ivec &b, const ivec &c);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec elem_mult(const svec &a, const svec &b, const svec &c);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec elem_mult(const bvec &a, const bvec &b, const bvec &c);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const vec &a, const vec &b,
                                     const vec &c, vec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const cvec &a, const cvec &b,
                                     const cvec &c, cvec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const ivec &a, const ivec &b,
                                     const ivec &c, ivec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const svec &a, const svec &b,
                                     const svec &c, svec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const bvec &a, const bvec &b,
                                     const bvec &c, bvec &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec elem_mult(const vec &a, const vec &b,
                                const vec &c, const vec &d);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec elem_mult(const cvec &a, const cvec &b,
                                 const cvec &c, const cvec &d);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec elem_mult(const ivec &a, const ivec &b,
                                 const ivec &c, const ivec &d);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec elem_mult(const svec &a, const svec &b,
                                 const svec &c, const svec &d);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec elem_mult(const bvec &a, const bvec &b,
                                 const bvec &c, const bvec &d);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const vec &a, const vec &b, const vec &c,
                                     const vec &d, vec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const cvec &a, const cvec &b,
                                     const cvec &c, const cvec &d, cvec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const ivec &a, const ivec &b,
                                     const ivec &c, const ivec &d, ivec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const svec &a, const svec &b,
                                     const svec &c, const svec &d, svec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_out(const bvec &a, const bvec &b,
                                     const bvec &c, const bvec &d, bvec &out);

// in-place elementwise multiplication

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_inplace(const vec &a, vec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_inplace(const cvec &a, cvec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_inplace(const ivec &a, ivec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_inplace(const svec &a, svec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_mult_inplace(const bvec &a, bvec &b);

// elementwise multiplication followed by summation

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT double elem_mult_sum(const vec &a, const vec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::complex<double> elem_mult_sum(const cvec &a,
    const cvec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT int elem_mult_sum(const ivec &a, const ivec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT short elem_mult_sum(const svec &a, const svec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bin elem_mult_sum(const bvec &a, const bvec &b);

// division operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator/(const vec &v, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator/(const cvec &v, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator/(const ivec &v, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator/(const svec &v, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator/(const bvec &v, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec operator/(double t, const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec operator/(std::complex<double> t, const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec operator/(int t, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec operator/(short t, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec operator/(bin t, const bvec &v);

// elementwise division operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec elem_div(const vec &a, const vec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec elem_div(const cvec &a, const cvec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec elem_div(const ivec &a, const ivec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec elem_div(const svec &a, const svec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec elem_div(const bvec &a, const bvec &b);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec elem_div(double t, const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec elem_div(std::complex<double> t, const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec elem_div(int t, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec elem_div(short t, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec elem_div(bin t, const bvec &v);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_div_out(const vec &a, const vec &b, vec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_div_out(const cvec &a, const cvec &b, cvec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_div_out(const ivec &a, const ivec &b, ivec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_div_out(const svec &a, const svec &b, svec &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT void elem_div_out(const bvec &a, const bvec &b, bvec &out);

// elementwise division followed by summation

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT double elem_div_sum(const vec &a, const vec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::complex<double> elem_div_sum(const cvec &a,
    const cvec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT int elem_div_sum(const ivec &a, const ivec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT short elem_div_sum(const svec &a, const svec &b);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bin elem_div_sum(const bvec &a, const bvec &b);

// concat operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(const vec &v, double a);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(const cvec &v, std::complex<double> a);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(const ivec &v, int a);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(const svec &v, short a);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(const bvec &v, bin a);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(double a, const vec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(std::complex<double> a, const cvec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(int a, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(short a, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(bin a, const bvec &v);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(const vec &v1, const vec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(const cvec &v1, const cvec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(const ivec &v1, const ivec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(const svec &v1, const svec &v2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(const bvec &v1, const bvec &v2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(const vec &v1, const vec &v2, const vec &v3);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(const cvec &v1, const cvec &v2, const cvec &v3);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(const ivec &v1, const ivec &v2, const ivec &v3);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(const svec &v1, const svec &v2, const svec &v3);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(const bvec &v1, const bvec &v2, const bvec &v3);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(const vec &v1, const vec &v2,
                             const vec &v3, const vec &v4);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(const cvec &v1, const cvec &v2,
                              const cvec &v3, const cvec &v4);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(const ivec &v1, const ivec &v2,
                              const ivec &v3, const ivec &v4);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(const svec &v1, const svec &v2,
                              const svec &v3, const svec &v4);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(const bvec &v1, const bvec &v2,
                              const bvec &v3, const bvec &v4);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT vec concat(const vec &v1, const vec &v2, const vec &v3,
                             const vec &v4, const vec &v5);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT cvec concat(const cvec &v1, const cvec &v2, const cvec &v3,
                              const cvec &v4, const cvec &v5);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT ivec concat(const ivec &v1, const ivec &v2, const ivec &v3,
                              const ivec &v4, const ivec &v5);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT svec concat(const svec &v1, const svec &v2, const svec &v3,
                              const svec &v4, const svec &v5);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT bvec concat(const bvec &v1, const bvec &v2, const bvec &v3,
                              const bvec &v4, const bvec &v5);

// I/O streams

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::ostream &operator<<(std::ostream& os, const vec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::ostream &operator<<(std::ostream& os, const cvec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::ostream &operator<<(std::ostream& os, const svec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::ostream &operator<<(std::ostream& os, const ivec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::ostream &operator<<(std::ostream& os, const bvec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::istream &operator>>(std::istream& is, vec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::istream &operator>>(std::istream& is, cvec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::istream &operator>>(std::istream& is, svec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::istream &operator>>(std::istream& is, ivec &vect);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT std::istream &operator>>(std::istream& is, bvec &vect);


} // namespace itpp

#endif // #ifndef VEC_H