IT++ Logo

vec.h

Go to the documentation of this file.
00001 
00030 #ifndef VEC_H
00031 #define VEC_H
00032 
00033 #ifndef _MSC_VER
00034 #  include <itpp/config.h>
00035 #else
00036 #  include <itpp/config_msvc.h>
00037 #endif
00038 
00039 #include <itpp/base/itassert.h>
00040 #include <itpp/base/math/misc.h>
00041 #include <itpp/base/copy_vector.h>
00042 #include <itpp/base/factory.h>
00043 
00044 
00045 namespace itpp {
00046 
00047   // Declaration of Vec
00048   template<class Num_T> class Vec;
00049   // Declaration of Mat
00050   template<class Num_T> class Mat;
00051   // Declaration of bin
00052   class bin;
00053 
00054   //-----------------------------------------------------------------------------------
00055   // Declaration of Vec Friends
00056   //-----------------------------------------------------------------------------------
00057 
00059   template<class Num_T> const Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00061   template<class Num_T> const Vec<Num_T> operator+(const Vec<Num_T> &v, const Num_T t);
00063   template<class Num_T> const Vec<Num_T> operator+(const Num_T t, const Vec<Num_T> &v);
00064 
00066   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00068   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v, const Num_T t);
00070   template<class Num_T> const Vec<Num_T> operator-(const Num_T t, const Vec<Num_T> &v);
00072   template<class Num_T> const Vec<Num_T> operator-(const Vec<Num_T> &v);
00073 
00075   template<class Num_T> Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00077   template<class Num_T> Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00086   template<class Num_T>
00087   const Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00088                                  bool hermitian = false);
00090   template<class Num_T> const Vec<Num_T> operator*(const Vec<Num_T> &v, const Num_T t);
00092   template<class Num_T> const Vec<Num_T> operator*(const Num_T t, const Vec<Num_T> &v);
00093 
00095   template<class Num_T> const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b);
00097   template<class Num_T> const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c);
00099   template<class Num_T> const 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);
00100 
00102   template<class Num_T> void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00104   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);
00106   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);
00107 
00109   template<class Num_T> void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b);
00111   template<class Num_T> Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);
00112 
00114   template<class Num_T> const Vec<Num_T> operator/(const Vec<Num_T> &v, const Num_T t);
00116   template<class Num_T> const Vec<Num_T> operator/(const Num_T t, const Vec<Num_T> &v);
00117 
00119   template<class Num_T> const Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b);
00121   template<class Num_T> const Vec<Num_T> elem_div(const Num_T t, const Vec<Num_T> &v);
00123   template<class Num_T> void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00125   template<class Num_T> Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b);
00126 
00128   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v, const Num_T a);
00130   template<class Num_T> const Vec<Num_T> concat(const Num_T a, const Vec<Num_T> &v);
00132   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1,const Vec<Num_T> &v2);
00134   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, const Vec<Num_T> &v3);
00136   template<class Num_T> const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, const Vec<Num_T> &v3, const Vec<Num_T> &v4);
00138   template<class Num_T> const 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);
00139 
00140   //-----------------------------------------------------------------------------------
00141   // Declaration of Vec
00142   //-----------------------------------------------------------------------------------
00143 
00207   template<class Num_T>
00208   class Vec {
00209   public:
00211     typedef Num_T value_type;
00212 
00214     explicit Vec(const Factory &f = DEFAULT_FACTORY);
00216     explicit Vec(int size, const Factory &f = DEFAULT_FACTORY);
00218     Vec(const Vec<Num_T> &v);
00220     Vec(const Vec<Num_T> &v, const Factory &f);
00222     Vec(const char *values, const Factory &f = DEFAULT_FACTORY);
00224     Vec(const std::string &values, const Factory &f = DEFAULT_FACTORY);
00226     Vec(const Num_T *c_array, int size, const Factory &f = DEFAULT_FACTORY);
00227 
00229     ~Vec();
00230 
00232     int length() const { return datasize; }
00234     int size() const { return datasize; }
00235 
00237     void set_size(int size, const bool copy=false);
00239     void set_length(int size, const bool copy=false) { set_size(size, copy); }
00241     void zeros();
00243     void clear() { zeros(); }
00245     void ones();
00247     void set(const char *str);
00249     void set(const std::string &str);
00250 
00252     const Num_T &operator[](int i) const;
00254     const Num_T &operator()(int i) const;
00256     Num_T &operator[](int i);
00258     Num_T &operator()(int i);
00260     const Vec<Num_T> operator()(int i1, int i2) const;
00262     const Vec<Num_T> operator()(const Vec<int> &indexlist) const;
00263 
00265     const Num_T &get(int i) const;
00267     const Vec<Num_T> get(int i1, int i2) const;
00269     void set(int i, const Num_T &v);
00270 
00272     Mat<Num_T> transpose() const;
00274     Mat<Num_T> T() const { return this->transpose(); }
00276     Mat<Num_T> hermitian_transpose() const;
00278     Mat<Num_T> H() const { return this->hermitian_transpose(); }
00279 
00281     Vec<Num_T>& operator+=(const Vec<Num_T> &v);
00283     Vec<Num_T>& operator+=(const Num_T t);
00285     friend const Vec<Num_T> operator+<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00287     friend const Vec<Num_T> operator+<>(const Vec<Num_T> &v, const Num_T t);
00289     friend const Vec<Num_T> operator+<>(const Num_T t, const Vec<Num_T> &v);
00290 
00292     Vec<Num_T>& operator-=(const Vec<Num_T> &v);
00294     Vec<Num_T>& operator-=(const Num_T t);
00296     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00298     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v, const Num_T t);
00300     friend const Vec<Num_T> operator-<>(const Num_T t, const Vec<Num_T> &v);
00302     friend const Vec<Num_T> operator-<>(const Vec<Num_T> &v);
00303 
00305     Vec<Num_T>& operator*=(const Num_T t);
00307     friend Num_T operator*<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00309     friend Num_T dot <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00311     friend const Mat<Num_T> outer_product <>(const Vec<Num_T> &v1,
00312                                              const Vec<Num_T> &v2,
00313                                              bool hermitian);
00315     friend const Vec<Num_T> operator*<>(const Vec<Num_T> &v, const Num_T t);
00317     friend const Vec<Num_T> operator*<>(const Num_T t, const Vec<Num_T> &v);
00318 
00320     friend const Vec<Num_T> elem_mult <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00322     friend const Vec<Num_T> elem_mult <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c);
00324     friend const 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);
00325 
00327     friend void elem_mult_out <>(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out);
00329     friend void elem_mult_out <>(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out);
00331     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);
00332 
00334     friend void elem_mult_inplace <>(const Vec<Num_T> &a, Vec<Num_T> &b);
00336     friend Num_T elem_mult_sum <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00337 
00339     Vec<Num_T>& operator/=(const Num_T t);
00341     Vec<Num_T>& operator/=(const Vec<Num_T> &v);
00342 
00344     friend const Vec<Num_T> operator/<>(const Vec<Num_T> &v, const Num_T t);
00346     friend const Vec<Num_T> operator/<>(const Num_T t, const Vec<Num_T> &v);
00347 
00349     friend const Vec<Num_T> elem_div <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2);
00351     friend const Vec<Num_T> elem_div <>(const Num_T t, const Vec<Num_T> &v);
00353     friend void elem_div_out <>(const Vec<Num_T> &v1, const Vec<Num_T> &v2, Vec<Num_T> &out);
00355     friend Num_T elem_div_sum <>(const Vec<Num_T> &a, const Vec<Num_T> &b);
00356 
00358     Vec<Num_T> get(const Vec<bin> &binlist) const;
00360     Vec<Num_T> right(int nr) const;
00362     Vec<Num_T> left(int nr) const;
00364     Vec<Num_T> mid(int start, int nr) const;
00366     Vec<Num_T> split(int pos);
00368     void shift_right(const Num_T In, int n=1);
00370     void shift_right(const Vec<Num_T> &In);
00372     void shift_left(const Num_T In, int n=1);
00374     void shift_left(const Vec<Num_T> &In);
00375 
00377     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v, const Num_T a);
00379     friend const Vec<Num_T> concat<>(const Num_T a, const Vec<Num_T> &v);
00381     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1,const Vec<Num_T> &v2);
00383     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00384                                      const Vec<Num_T> &v3);
00386     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00387                                      const Vec<Num_T> &v3, const Vec<Num_T> &v4);
00389     friend const Vec<Num_T> concat<>(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
00390                                      const Vec<Num_T> &v3, const Vec<Num_T> &v4,
00391                                      const Vec<Num_T> &v5);
00392 
00394     void set_subvector(int i1, int i2, const Vec<Num_T> &v);
00396     void set_subvector(int i, const Vec<Num_T> &v);
00398     void set_subvector(int i1, int i2, const Num_T t);
00400     void replace_mid(int pos, const Vec<Num_T> &v);
00402     void del(int index);
00404     void del(int i1, int i2);
00406     void ins(int index, const Num_T in);
00408     void ins(int index, const Vec<Num_T> &in);
00409 
00411     Vec<Num_T>& operator=(const Num_T t);
00413     Vec<Num_T>& operator=(const Vec<Num_T> &v);
00415     Vec<Num_T>& operator=(const Mat<Num_T> &m);
00417     Vec<Num_T>& operator=(const char *values);
00418 
00420     Vec<bin> operator==(const Num_T value) const;
00422     Vec<bin> operator!=(const Num_T value) const;
00424     Vec<bin> operator<(const Num_T value) const;
00426     Vec<bin> operator<=(const Num_T value) const;
00428     Vec<bin> operator>(const Num_T value) const;
00430     Vec<bin> operator>=(const Num_T value) const;
00431 
00433     bool operator==(const Vec<Num_T> &v) const;
00435     bool operator!=(const Vec<Num_T> &v) const;
00436 
00438     Num_T &_elem(int i) { return data[i]; }
00440     const Num_T &_elem(int i) const { return data[i]; }
00441 
00443     Num_T *_data() { return data; }
00445     const Num_T *_data() const { return data; }
00446 
00447   protected:
00449     void alloc(int size);
00451     void free();
00452 
00454     int datasize;
00456     Num_T *data;
00458     const Factory &factory;
00459   private:
00460     // This function is used in set() methods to replace commas with spaces
00461     std::string replace_commas(const std::string &str);
00462   };
00463 
00464   //-----------------------------------------------------------------------------------
00465   // Type definitions of vec, cvec, ivec, svec, and bvec
00466   //-----------------------------------------------------------------------------------
00467 
00472   typedef Vec<double> vec;
00473 
00478   typedef Vec<std::complex<double> > cvec;
00479 
00484   typedef Vec<int> ivec;
00485 
00490   typedef Vec<short int> svec;
00491 
00496   typedef Vec<bin> bvec;
00497 
00498 } //namespace itpp
00499 
00500 
00501 #include <itpp/base/mat.h>
00502 
00503 namespace itpp {
00504 
00505   //-----------------------------------------------------------------------------------
00506   // Declaration of input and output streams for Vec
00507   //-----------------------------------------------------------------------------------
00508 
00513   template<class Num_T>
00514   std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v);
00515 
00527   template<class Num_T>
00528   std::istream &operator>>(std::istream &is, Vec<Num_T> &v);
00529 
00530   //-----------------------------------------------------------------------------------
00531   // Implementation of templated Vec members and friends
00532   //-----------------------------------------------------------------------------------
00533 
00534   template<class Num_T> inline
00535   void Vec<Num_T>::alloc(int size)
00536   {
00537     if (size > 0) {
00538       create_elements(data, size, factory);
00539       datasize = size;
00540     }
00541     else {
00542       data = 0;
00543       datasize = 0;
00544     }
00545   }
00546 
00547   template<class Num_T> inline
00548   void Vec<Num_T>::free()
00549   {
00550     destroy_elements(data, datasize);
00551     datasize = 0;
00552   }
00553 
00554 
00555   template<class Num_T> inline
00556   Vec<Num_T>::Vec(const Factory &f) : datasize(0), data(0), factory(f) {}
00557 
00558   template<class Num_T> inline
00559   Vec<Num_T>::Vec(int size, const Factory &f) : datasize(0), data(0), factory(f)
00560   {
00561     it_assert_debug(size>=0, "Negative size in Vec::Vec(int)");
00562     alloc(size);
00563   }
00564 
00565   template<class Num_T> inline
00566   Vec<Num_T>::Vec(const Vec<Num_T> &v) : datasize(0), data(0), factory(v.factory)
00567   {
00568     alloc(v.datasize);
00569     copy_vector(datasize, v.data, data);
00570   }
00571 
00572   template<class Num_T> inline
00573   Vec<Num_T>::Vec(const Vec<Num_T> &v, const Factory &f) : datasize(0), data(0), factory(f)
00574   {
00575     alloc(v.datasize);
00576     copy_vector(datasize, v.data, data);
00577   }
00578 
00579   template<class Num_T> inline
00580   Vec<Num_T>::Vec(const char *values, const Factory &f) : datasize(0), data(0), factory(f)
00581   {
00582     set(values);
00583   }
00584 
00585   template<class Num_T> inline
00586   Vec<Num_T>::Vec(const std::string &values, const Factory &f) : datasize(0), data(0), factory(f)
00587   {
00588     set(values);
00589   }
00590 
00591   template<class Num_T> inline
00592   Vec<Num_T>::Vec(const Num_T *c_array, int size, const Factory &f) : datasize(0), data(0), factory(f)
00593   {
00594     alloc(size);
00595     copy_vector(size, c_array, data);
00596   }
00597 
00598   template<class Num_T> inline
00599   Vec<Num_T>::~Vec()
00600   {
00601     free();
00602   }
00603 
00604   template<class Num_T>
00605   void Vec<Num_T>::set_size(int size, bool copy)
00606   {
00607     it_assert_debug(size >= 0, "Vec::set_size(): New size must not be negative");
00608     if (datasize == size)
00609       return;
00610     if (copy) {
00611       // create a temporary pointer to the allocated data
00612       Num_T* tmp = data;
00613       // store the current number of elements
00614       int old_datasize = datasize;
00615       // check how many elements we need to copy
00616       int min = datasize < size ? datasize : size;
00617       // allocate new memory
00618       alloc(size);
00619       // copy old elements into a new memory region
00620       copy_vector(min, tmp, data);
00621       // initialize the rest of resized vector
00622       for (int i = min; i < size; ++i)
00623         data[i] = Num_T(0);
00624       // delete old elements
00625       destroy_elements(tmp, old_datasize);
00626     }
00627     else {
00628       free();
00629       alloc(size);
00630     }
00631   }
00632 
00633   template<class Num_T> inline
00634   const Num_T& Vec<Num_T>::operator[](int i) const
00635   {
00636     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator[]: Index out of range");
00637     return data[i];
00638   }
00639 
00640   template<class Num_T> inline
00641   const Num_T& Vec<Num_T>::operator()(int i) const
00642   {
00643     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator(): Index out of range");
00644     return data[i];
00645   }
00646 
00647   template<class Num_T> inline
00648   Num_T& Vec<Num_T>::operator[](int i)
00649   {
00650     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator[]: Index out of range");
00651     return data[i];
00652   }
00653 
00654   template<class Num_T> inline
00655   Num_T& Vec<Num_T>::operator()(int i)
00656   {
00657     it_assert_debug((i >= 0) && (i < datasize), "Vec::operator(): Index out of range");
00658     return data[i];
00659   }
00660 
00661   template<class Num_T> inline
00662   const Vec<Num_T> Vec<Num_T>::operator()(int i1, int i2) const
00663   {
00664     int ii1=i1, ii2=i2;
00665 
00666     if (ii1 == -1) ii1 = datasize-1;
00667     if (ii2 == -1) ii2 = datasize-1;
00668 
00669     it_assert_debug(ii1>=0 && ii2>=0 && ii1<datasize && ii2<datasize, "Vec::operator()(i1,i2): indicies out of range");
00670     it_assert_debug(ii2>=ii1, "Vec::op(i1,i2): i2 >= i1 necessary");
00671 
00672     Vec<Num_T> s(ii2-ii1+1);
00673     copy_vector(s.datasize, data+ii1, s.data);
00674 
00675     return s;
00676   }
00677 
00678   template<class Num_T>
00679   const Vec<Num_T> Vec<Num_T>::operator()(const Vec<int> &indexlist) const
00680   {
00681     Vec<Num_T> temp(indexlist.length());
00682     for (int i=0;i<indexlist.length();i++) {
00683       it_assert((indexlist(i)>=0) && (indexlist(i) < datasize), "Vec::operator()(ivec &): index outside range");
00684       temp(i)=data[indexlist(i)];
00685     }
00686     return temp;
00687   }
00688 
00689 
00690   template<class Num_T> inline
00691   const Num_T& Vec<Num_T>::get(int i) const
00692   {
00693     it_assert_debug((i >= 0) && (i < datasize), "method get()");
00694     return data[i];
00695   }
00696 
00697   template<class Num_T> inline
00698   const Vec<Num_T> Vec<Num_T>::get(int i1, int i2) const
00699   {
00700     return (*this)(i1, i2);
00701   }
00702 
00703 
00704   template<class Num_T> inline
00705   void Vec<Num_T>::zeros()
00706   {
00707     for (int i = 0; i < datasize; i++)
00708       data[i] = Num_T(0);
00709   }
00710 
00711   template<class Num_T> inline
00712   void Vec<Num_T>::ones()
00713   {
00714     for (int i = 0; i < datasize; i++)
00715       data[i] = Num_T(1);
00716   }
00717 
00718   template<class Num_T> inline
00719   void Vec<Num_T>::set(int i, const Num_T &v)
00720   {
00721     it_assert_debug((i >= 0) && (i < datasize), "method set()");
00722     data[i] = v;
00723   }
00724 
00726   template<>
00727   void Vec<double>::set(const std::string &str);
00728   template<>
00729   void Vec<std::complex<double> >::set(const std::string &str);
00730   template<>
00731   void Vec<int>::set(const std::string &str);
00732   template<>
00733   void Vec<short int>::set(const std::string &str);
00734   template<>
00735   void Vec<bin>::set(const std::string &str);
00737 
00738   template<class Num_T>
00739   void Vec<Num_T>::set(const std::string &str)
00740   {
00741     it_error("Vec::set(): Only `double', `complex<double>', `int', "
00742              "`short int' and `bin' types supported");
00743   }
00744 
00745   template<class Num_T>
00746   void Vec<Num_T>::set(const char *str)
00747   {
00748     set(std::string(str));
00749   }
00750 
00751 
00752   template<class Num_T>
00753   Mat<Num_T> Vec<Num_T>::transpose() const
00754   {
00755     Mat<Num_T> temp(1, datasize);
00756     copy_vector(datasize, data, temp._data());
00757     return temp;
00758   }
00759 
00760   template<>
00761   Mat<std::complex<double> > Vec<std::complex<double> >::hermitian_transpose() const;
00762 
00763   template<class Num_T>
00764   Mat<Num_T> Vec<Num_T>::hermitian_transpose() const
00765   {
00766     Mat<Num_T> temp(1, datasize);
00767     copy_vector(datasize, data, temp._data());
00768     return temp;
00769   }
00770 
00771   template<class Num_T> inline
00772   Vec<Num_T>& Vec<Num_T>::operator+=(const Vec<Num_T> &v)
00773   {
00774     if (datasize == 0) { // if not assigned a size.
00775       if (this != &v) { // check for self addition
00776         alloc(v.datasize);
00777         copy_vector(datasize, v.data, data);
00778       }
00779     } else {
00780       it_assert_debug(datasize == v.datasize, "Vec::operator+=: Wrong sizes");
00781       for (int i = 0; i < datasize; i++)
00782         data[i] += v.data[i];
00783     }
00784     return *this;
00785   }
00786 
00787   template<class Num_T> inline
00788   Vec<Num_T>& Vec<Num_T>::operator+=(const Num_T t)
00789   {
00790     for (int i=0;i<datasize;i++)
00791       data[i]+=t;
00792     return *this;
00793   }
00794 
00795   template<class Num_T> inline
00796   const Vec<Num_T> operator+(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00797   {
00798     int i;
00799     Vec<Num_T> r(v1.datasize);
00800 
00801     it_assert_debug(v1.datasize==v2.datasize, "Vec::operator+: wrong sizes");
00802     for (i=0; i<v1.datasize; i++)
00803       r.data[i] = v1.data[i] + v2.data[i];
00804 
00805     return r;
00806   }
00807 
00808   template<class Num_T> inline
00809   const Vec<Num_T> operator+(const Vec<Num_T> &v, const Num_T t)
00810   {
00811     int i;
00812     Vec<Num_T> r(v.datasize);
00813 
00814     for (i=0; i<v.datasize; i++)
00815       r.data[i] = v.data[i] + t;
00816 
00817     return r;
00818   }
00819 
00820   template<class Num_T> inline
00821   const Vec<Num_T> operator+(const Num_T t, const Vec<Num_T> &v)
00822   {
00823     int i;
00824     Vec<Num_T> r(v.datasize);
00825 
00826     for (i=0; i<v.datasize; i++)
00827       r.data[i] = t + v.data[i];
00828 
00829     return r;
00830   }
00831 
00832   template<class Num_T> inline
00833   Vec<Num_T>& Vec<Num_T>::operator-=(const Vec<Num_T> &v)
00834   {
00835     if (datasize == 0) { // if not assigned a size.
00836       if (this != &v) { // check for self decrementation
00837         alloc(v.datasize);
00838         for (int i = 0; i < v.datasize; i++)
00839           data[i] = -v.data[i];
00840       }
00841     } else {
00842       it_assert_debug(datasize == v.datasize, "Vec::operator-=: Wrong sizes");
00843       for (int i = 0; i < datasize; i++)
00844         data[i] -= v.data[i];
00845     }
00846     return *this;
00847   }
00848 
00849   template<class Num_T> inline
00850   Vec<Num_T>& Vec<Num_T>::operator-=(const Num_T t)
00851   {
00852     for (int i=0;i<datasize;i++)
00853       data[i]-=t;
00854     return *this;
00855   }
00856 
00857   template<class Num_T> inline
00858   const Vec<Num_T> operator-(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00859   {
00860     int i;
00861     Vec<Num_T> r(v1.datasize);
00862 
00863     it_assert_debug(v1.datasize==v2.datasize, "Vec::operator-: wrong sizes");
00864     for (i=0; i<v1.datasize; i++)
00865       r.data[i] = v1.data[i] - v2.data[i];
00866 
00867     return r;
00868   }
00869 
00870   template<class Num_T> inline
00871   const Vec<Num_T> operator-(const Vec<Num_T> &v, const Num_T t)
00872   {
00873     int i;
00874     Vec<Num_T> r(v.datasize);
00875 
00876     for (i=0; i<v.datasize; i++)
00877       r.data[i] = v.data[i] - t;
00878 
00879     return r;
00880   }
00881 
00882   template<class Num_T> inline
00883   const Vec<Num_T> operator-(const Num_T t, const Vec<Num_T> &v)
00884   {
00885     int i;
00886     Vec<Num_T> r(v.datasize);
00887 
00888     for (i=0; i<v.datasize; i++)
00889       r.data[i] = t - v.data[i];
00890 
00891     return r;
00892   }
00893 
00894   template<class Num_T> inline
00895   const Vec<Num_T> operator-(const Vec<Num_T> &v)
00896   {
00897     int i;
00898     Vec<Num_T> r(v.datasize);
00899 
00900     for (i=0; i<v.datasize; i++)
00901       r.data[i] = -v.data[i];
00902 
00903     return r;
00904   }
00905 
00906   template<class Num_T> inline
00907   Vec<Num_T>& Vec<Num_T>::operator*=(const Num_T t)
00908   {
00909     scal_vector(datasize, t, data);
00910     return *this;
00911   }
00912 
00913 #if defined(HAVE_BLAS)
00914   template<> inline
00915   double dot(const vec &v1, const vec &v2)
00916   {
00917     it_assert_debug(v1.datasize == v2.datasize, "vec::dot: wrong sizes");
00918     int incr = 1;
00919     return blas::ddot_(&v1.datasize, v1.data, &incr, v2.data, &incr);
00920   }
00921 
00922   template<> inline
00923   std::complex<double> dot(const cvec &v1, const cvec &v2)
00924   {
00925     it_assert_debug(v1.datasize == v2.datasize, "cvec::dot: wrong sizes");
00926     int incr = 1;
00927     std::complex<double> output;
00928 #if defined(HAVE_ZDOTU_VOID)
00929     blas::zdotu_(&output, &v1.datasize, v1.data, &incr, v2.data, &incr);
00930 #elif defined(HAVE_ZDOTU_RETURN)
00931     output = blas::zdotu_(&v1.datasize, v1.data, &incr, v2.data, &incr);
00932 #else
00933     blas::zdotusub_(&output, &v1.datasize, v1.data, &incr, v2.data, &incr);
00934 #endif // HAVE_ZDOTU_VOID
00935     return output;
00936   }
00937 #endif // HAVE_BLAS
00938 
00939   template<class Num_T> inline
00940   Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00941   {
00942     int i;
00943     Num_T r=Num_T(0);
00944 
00945     it_assert_debug(v1.datasize==v2.datasize, "Vec::dot: wrong sizes");
00946     for (i=0; i<v1.datasize; i++)
00947       r += v1.data[i] * v2.data[i];
00948 
00949     return r;
00950   }
00951 
00952   template<class Num_T> inline
00953   Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
00954   {
00955     return dot(v1, v2);
00956   }
00957 
00958 
00959 #if defined(HAVE_BLAS)
00960   template<> inline
00961   const mat outer_product(const vec &v1, const vec &v2, bool hermitian)
00962   {
00963     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
00964                     "Vec::outer_product():: Input vector of zero size");
00965 
00966     mat out(v1.datasize, v2.datasize);
00967     out.zeros();
00968     double alpha = 1.0;
00969     int incr = 1;
00970     blas::dger_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00971                 v2.data, &incr, out._data(), &v1.datasize);
00972     return out;
00973   }
00974 
00975   template<> inline
00976   const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian)
00977   {
00978     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
00979                     "Vec::outer_product():: Input vector of zero size");
00980 
00981     cmat out(v1.datasize, v2.datasize);
00982     out.zeros();
00983     std::complex<double> alpha(1.0);
00984     int incr = 1;
00985     if (hermitian) {
00986       blas::zgerc_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00987                    v2.data, &incr, out._data(), &v1.datasize);
00988     }
00989     else {
00990       blas::zgeru_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr,
00991                    v2.data, &incr, out._data(), &v1.datasize);
00992     }
00993     return out;
00994   }
00995 #else
00997   template<> inline
00998   const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian)
00999   {
01000     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
01001                     "Vec::outer_product():: Input vector of zero size");
01002 
01003     cmat out(v1.datasize, v2.datasize);
01004     if (hermitian) {
01005       for (int i = 0; i < v1.datasize; ++i) {
01006         for (int j = 0; j < v2.datasize; ++j) {
01007           out(i, j) = v1.data[i] * conj(v2.data[j]);
01008         }
01009       }
01010     }
01011     else {
01012       for (int i = 0; i < v1.datasize; ++i) {
01013         for (int j = 0; j < v2.datasize; ++j) {
01014           out(i, j) = v1.data[i] * v2.data[j];
01015         }
01016       }
01017     }
01018     return out;
01019   }
01020 #endif // HAVE_BLAS
01021 
01022   template<class Num_T> inline
01023   const Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01024                                  bool hermitian)
01025   {
01026     int i, j;
01027 
01028     it_assert_debug((v1.datasize > 0) && (v2.datasize > 0),
01029                     "Vec::outer_product:: Input vector of zero size");
01030 
01031     Mat<Num_T> r(v1.datasize, v2.datasize);
01032     for (i=0; i<v1.datasize; i++) {
01033       for (j=0; j<v2.datasize; j++) {
01034         r(i,j) = v1.data[i] * v2.data[j];
01035       }
01036     }
01037 
01038     return r;
01039   }
01040 
01041   template<class Num_T> inline
01042   const Vec<Num_T> operator*(const Vec<Num_T> &v, const Num_T t)
01043   {
01044     int i;
01045     Vec<Num_T> r(v.datasize);
01046     for (i=0; i<v.datasize; i++)
01047       r.data[i] = v.data[i] * t;
01048 
01049     return r;
01050   }
01051 
01052   template<class Num_T> inline
01053   const Vec<Num_T> operator*(const Num_T t, const Vec<Num_T> &v)
01054   {
01055     return operator*(v, t);
01056   }
01057 
01058   template<class Num_T> inline
01059   const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b)
01060   {
01061     Vec<Num_T> out;
01062     elem_mult_out(a,b,out);
01063     return out;
01064   }
01065 
01066   template<class Num_T> inline
01067   const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c)
01068   {
01069     Vec<Num_T> out;
01070     elem_mult_out(a,b,c,out);
01071     return out;
01072   }
01073 
01074   template<class Num_T> inline
01075   const 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)
01076   {
01077     Vec<Num_T> out;
01078     elem_mult_out(a,b,c,d,out);
01079     return out;
01080   }
01081 
01082   template<class Num_T> inline
01083   void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
01084   {
01085     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_out: wrong sizes");
01086 
01087     if(out.datasize != a.datasize)
01088       out.set_size(a.size());
01089 
01090     for(int i=0; i<a.datasize; i++)
01091       out.data[i] = a.data[i] * b.data[i];
01092   }
01093 
01094   template<class Num_T> inline
01095   void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out)
01096   {
01097     it_assert_debug(a.datasize==b.datasize==c.datasize, "Vec::elem_mult_out: wrong sizes");
01098 
01099     if(out.datasize != a.datasize)
01100       out.set_size(a.size());
01101 
01102     for(int i=0; i<a.datasize; i++)
01103       out.data[i] = a.data[i] * b.data[i] * c.data[i];
01104   }
01105 
01106   template<class Num_T> inline
01107   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)
01108   {
01109     it_assert_debug(a.datasize==b.datasize==c.datasize==d.datasize, "Vec::elem_mult_out: wrong sizes");
01110 
01111     if(out.datasize != a.datasize)
01112       out.set_size(a.size());
01113 
01114     for(int i=0; i<a.datasize; i++)
01115       out.data[i] = a.data[i] * b.data[i] * c.data[i] * d.data[i];
01116   }
01117 
01118   template<class Num_T> inline
01119   void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b)
01120   {
01121     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_inplace: wrong sizes");
01122 
01123     for(int i=0; i<a.datasize; i++)
01124       b.data[i] *= a.data[i];
01125   }
01126 
01127   template<class Num_T> inline
01128   Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
01129   {
01130     it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_sum: wrong sizes");
01131 
01132     Num_T acc = 0;
01133 
01134     for(int i=0; i<a.datasize; i++)
01135       acc += a.data[i] * b.data[i];
01136 
01137     return acc;
01138   }
01139 
01140   template<class Num_T> inline
01141   const Vec<Num_T> operator/(const Vec<Num_T> &v, const Num_T t)
01142   {
01143     int i;
01144     Vec<Num_T> r(v.datasize);
01145 
01146     for (i=0; i<v.datasize; i++)
01147       r.data[i] = v.data[i] / t;
01148 
01149     return r;
01150   }
01151 
01152   template<class Num_T> inline
01153   const Vec<Num_T> operator/(const Num_T t, const Vec<Num_T> &v)
01154   {
01155     int i;
01156     Vec<Num_T> r(v.datasize);
01157 
01158     for (i=0; i<v.datasize; i++)
01159       r.data[i] = t / v.data[i];
01160 
01161     return r;
01162   }
01163 
01164   template<class Num_T> inline
01165   Vec<Num_T>& Vec<Num_T>::operator/=(const Num_T t)
01166   {
01167     for (int i = 0; i < datasize; ++i) {
01168       data[i] /= t;
01169     }
01170     return *this;
01171   }
01172 
01173   template<class Num_T> inline
01174   Vec<Num_T>& Vec<Num_T>::operator/=(const Vec<Num_T> &v)
01175   {
01176     it_assert_debug(datasize == v.datasize, "Vec::operator/=(): wrong sizes");
01177     for (int i = 0; i < datasize; ++i) {
01178       data[i] /= v.data[i];
01179     }
01180     return *this;
01181   }
01182 
01183   template<class Num_T> inline
01184   const Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b)
01185   {
01186     Vec<Num_T> out;
01187     elem_div_out(a,b,out);
01188     return out;
01189   }
01190 
01191   template<class Num_T> inline
01192   const Vec<Num_T> elem_div(const Num_T t, const Vec<Num_T> &v)
01193   {
01194     int i;
01195     Vec<Num_T> r(v.datasize);
01196 
01197     for (i=0; i<v.datasize; i++)
01198       r.data[i] = t / v.data[i];
01199 
01200     return r;
01201   }
01202 
01203   template<class Num_T> inline
01204   void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out)
01205   {
01206     it_assert_debug(a.datasize==b.datasize, "Vecelem_div_out: wrong sizes");
01207 
01208     if(out.datasize != a.datasize)
01209       out.set_size(a.size());
01210 
01211     for(int i=0; i<a.datasize; i++)
01212       out.data[i] = a.data[i] / b.data[i];
01213   }
01214 
01215   template<class Num_T> inline
01216   Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b)
01217   {
01218     it_assert_debug(a.datasize==b.datasize, "Vec::elem_div_sum: wrong sizes");
01219 
01220     Num_T acc = 0;
01221 
01222     for(int i=0; i<a.datasize; i++)
01223       acc += a.data[i] / b.data[i];
01224 
01225     return acc;
01226   }
01227 
01228   template<class Num_T>
01229   Vec<Num_T> Vec<Num_T>::get(const Vec<bin> &binlist) const
01230   {
01231     int size = binlist.size();
01232     it_assert_debug(datasize == size, "Vec::get(bvec &): wrong sizes");
01233     Vec<Num_T> temp(size);
01234     int j = 0;
01235     for (int i = 0; i < size; ++i) {
01236       if (binlist(i) == bin(1)) {
01237         temp(j) = data[i];
01238         j++;
01239       }
01240     }
01241     temp.set_size(j, true);
01242     return temp;
01243   }
01244 
01245   template<class Num_T> inline
01246   Vec<Num_T> Vec<Num_T>::right(int nr) const
01247   {
01248     it_assert_debug(nr <= datasize, "Vec::right(): index out of range");
01249     Vec<Num_T> temp(nr);
01250     if (nr > 0) {
01251       copy_vector(nr, &data[datasize-nr], temp.data);
01252     }
01253     return temp;
01254   }
01255 
01256   template<class Num_T> inline
01257   Vec<Num_T> Vec<Num_T>::left(int nr) const
01258   {
01259     it_assert_debug(nr <= datasize, "Vec::left(): index out of range");
01260     Vec<Num_T> temp(nr);
01261     if (nr > 0) {
01262       copy_vector(nr, data, temp.data);
01263     }
01264     return temp;
01265   }
01266 
01267   template<class Num_T> inline
01268   Vec<Num_T> Vec<Num_T>::mid(int start, int nr) const
01269   {
01270     it_assert_debug((start >= 0) && ((start+nr) <= datasize),
01271                     "Vec::mid(): indexing out of range");
01272     Vec<Num_T> temp(nr);
01273     if (nr > 0) {
01274       copy_vector(nr, &data[start], temp.data);
01275     }
01276     return temp;
01277   }
01278 
01279   template<class Num_T>
01280   Vec<Num_T> Vec<Num_T>::split(int pos)
01281   {
01282     it_assert_debug((pos >= 0) && (pos <= datasize), "Vec::split(): index out of range");
01283     Vec<Num_T> temp1(pos);
01284     Vec<Num_T> temp2(datasize-pos);
01285     copy_vector(pos, data, temp1.data);
01286     copy_vector(datasize-pos, &data[pos], temp2.data);
01287     (*this) = temp2;
01288     return temp1;
01289   }
01290 
01291   template<class Num_T>
01292   void Vec<Num_T>::shift_right(const Num_T In, int n)
01293   {
01294     int i=datasize;
01295 
01296     it_assert_debug(n>=0, "Vec::shift_right: index out of range");
01297     while (--i >= n)
01298       data[i] = data[i-n];
01299     while (i >= 0)
01300       data[i--] = In;
01301   }
01302 
01303   template<class Num_T>
01304   void Vec<Num_T>::shift_right(const Vec<Num_T> &In)
01305   {
01306     int i;
01307 
01308     for (i=datasize-1; i>=In.datasize; i--)
01309       data[i]=data[i-In.datasize];
01310     for (i=0; i<In.datasize; i++)
01311       data[i]=In[i];
01312   }
01313 
01314   template<class Num_T>
01315   void Vec<Num_T>::shift_left(const Num_T In, int n)
01316   {
01317     int i;
01318 
01319     it_assert_debug(n>=0, "Vec::shift_left: index out of range");
01320     for (i=0; i<datasize-n; i++)
01321       data[i] = data[i+n];
01322     while (i < datasize)
01323       data[i++] = In;
01324   }
01325 
01326   template<class Num_T>
01327   void Vec<Num_T>::shift_left(const Vec<Num_T> &In)
01328   {
01329     int i;
01330 
01331     for (i=0; i<datasize-In.datasize; i++)
01332       data[i]=data[i+In.datasize];
01333     for (i=datasize-In.datasize; i<datasize; i++)
01334       data[i]=In[i-datasize+In.datasize];
01335   }
01336 
01337   template<class Num_T>
01338   const Vec<Num_T> concat(const Vec<Num_T> &v, const Num_T a)
01339   {
01340     int size = v.size();
01341     Vec<Num_T> temp(size+1);
01342     copy_vector(size, v.data, temp.data);
01343     temp(size) = a;
01344     return temp;
01345   }
01346 
01347   template<class Num_T>
01348   const Vec<Num_T> concat(const Num_T a, const Vec<Num_T> &v)
01349   {
01350     int size = v.size();
01351     Vec<Num_T> temp(size+1);
01352     temp(0) = a;
01353     copy_vector(size, v.data, &temp.data[1]);
01354     return temp;
01355   }
01356 
01357   template<class Num_T>
01358   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2)
01359   {
01360     int size1 = v1.size();
01361     int size2 = v2.size();
01362     Vec<Num_T> temp(size1+size2);
01363     copy_vector(size1, v1.data, temp.data);
01364     copy_vector(size2, v2.data, &temp.data[size1]);
01365     return temp;
01366   }
01367 
01368   template<class Num_T>
01369   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01370                           const Vec<Num_T> &v3)
01371   {
01372     int size1 = v1.size();
01373     int size2 = v2.size();
01374     int size3 = v3.size();
01375     Vec<Num_T> temp(size1+size2+size3);
01376     copy_vector(size1, v1.data, temp.data);
01377     copy_vector(size2, v2.data, &temp.data[size1]);
01378     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01379     return temp;
01380   }
01381 
01382   template<class Num_T>
01383   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01384                           const Vec<Num_T> &v3, const Vec<Num_T> &v4)
01385   {
01386     int size1 = v1.size();
01387     int size2 = v2.size();
01388     int size3 = v3.size();
01389     int size4 = v4.size();
01390     Vec<Num_T> temp(size1+size2+size3+size4);
01391     copy_vector(size1, v1.data, temp.data);
01392     copy_vector(size2, v2.data, &temp.data[size1]);
01393     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01394     copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
01395     return temp;
01396   }
01397 
01398   template<class Num_T>
01399   const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2,
01400                           const Vec<Num_T> &v3, const Vec<Num_T> &v4,
01401                           const Vec<Num_T> &v5)
01402   {
01403     int size1 = v1.size();
01404     int size2 = v2.size();
01405     int size3 = v3.size();
01406     int size4 = v4.size();
01407     int size5 = v5.size();
01408     Vec<Num_T> temp(size1+size2+size3+size4+size5);
01409     copy_vector(size1, v1.data, temp.data);
01410     copy_vector(size2, v2.data, &temp.data[size1]);
01411     copy_vector(size3, v3.data, &temp.data[size1+size2]);
01412     copy_vector(size4, v4.data, &temp.data[size1+size2+size3]);
01413     copy_vector(size5, v5.data, &temp.data[size1+size2+size3+size4]);
01414     return temp;
01415   }
01416 
01417   template<class Num_T> inline
01418   void Vec<Num_T>::set_subvector(int i1, int i2, const Vec<Num_T> &v)
01419   {
01420     if (i1 == -1) i1 = datasize-1;
01421     if (i2 == -1) i2 = datasize-1;
01422 
01423     it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range");
01424     it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary");
01425     it_assert_debug(i2-i1+1 == v.datasize, "Vec::set_subvector(): wrong sizes");
01426 
01427     copy_vector(v.datasize, v.data, data+i1);
01428   }
01429 
01430   template<class Num_T> inline
01431   void Vec<Num_T>:: set_subvector(int i, const Vec<Num_T> &v)
01432   {
01433     it_assert_debug(i>=0, "Vec::set_subvector(): index out of range");
01434     it_assert_debug(i+v.datasize <= datasize, "Vec::set_subvector(): too long input vector");
01435     copy_vector(v.datasize, v.data, data+i);
01436   }
01437 
01438   template<class Num_T>
01439   void Vec<Num_T>::set_subvector(int i1, int i2, const Num_T t)
01440   {
01441     if (i1 == -1) i1 = datasize-1;
01442     if (i2 == -1) i2 = datasize-1;
01443 
01444     it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range");
01445     it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary");
01446 
01447     for (int i=i1;i<=i2;i++)
01448       data[i] = t;
01449   }
01450 
01451   template<class Num_T>
01452   void Vec<Num_T>::replace_mid(int pos, const Vec<Num_T> &v)
01453   {
01454     it_assert_debug((pos>=0) && ((pos+v.length())<=datasize), "Vec::replace_mid: indexing out of range");
01455     copy_vector(v.datasize, v.data, &data[pos]);
01456   }
01457 
01458   template<class Num_T>
01459   void Vec<Num_T>::del(int index)
01460   {
01461     it_assert_debug((index >= 0) && (index < datasize),
01462                     "Vec::del(): index out of range");
01463     Vec<Num_T> temp(*this);
01464     set_size(datasize-1, false);
01465     copy_vector(index, temp.data, data);
01466     copy_vector(datasize-index, &temp.data[index+1], &data[index]);
01467   }
01468 
01469   template<class Num_T>
01470   void  Vec<Num_T>::del(int i1, int i2)
01471   {
01472     it_assert_debug((i1 >= 0) && (i2 < datasize) && (i1 < i2),
01473                     "Vec::del(): index out of range");
01474     Vec<Num_T> temp(*this);
01475     int new_size = datasize-(i2-i1+1);
01476     set_size(new_size, false);
01477     copy_vector(i1, temp.data, data);
01478     copy_vector(datasize-i1, &temp.data[i2+1], &data[i1]);
01479   }
01480 
01481   template<class Num_T>
01482   void Vec<Num_T>::ins(int index, const Num_T in)
01483   {
01484     it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range");
01485     Vec<Num_T> Temp(*this);
01486 
01487     set_size(datasize+1, false);
01488     copy_vector(index, Temp.data, data);
01489     data[index]=in;
01490     copy_vector(Temp.datasize-index, Temp.data+index, data+index+1);
01491   }
01492 
01493   template<class Num_T>
01494   void Vec<Num_T>::ins(int index, const Vec<Num_T> &in)
01495   {
01496     it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range");
01497     Vec<Num_T> Temp(*this);
01498 
01499     set_size(datasize+in.length(), false);
01500     copy_vector(index, Temp.data, data);
01501     copy_vector(in.size(), in.data, &data[index]);
01502     copy_vector(Temp.datasize-index, Temp.data+index, data+index+in.size());
01503   }
01504 
01505   template<class Num_T> inline
01506   Vec<Num_T>& Vec<Num_T>::operator=(const Num_T t)
01507   {
01508     for (int i=0;i<datasize;i++)
01509       data[i] = t;
01510     return *this;
01511   }
01512 
01513   template<class Num_T> inline
01514   Vec<Num_T>& Vec<Num_T>::operator=(const Vec<Num_T> &v)
01515   {
01516     if (this != &v) {
01517       set_size(v.datasize, false);
01518       copy_vector(datasize, v.data, data);
01519     }
01520     return *this;
01521   }
01522 
01523   template<class Num_T> inline
01524   Vec<Num_T>& Vec<Num_T>::operator=(const Mat<Num_T> &m)
01525   {
01526     it_assert_debug( (m.cols() == 1 && datasize == m.rows()) ||
01527                 (m.rows() == 1 && datasize == m.cols()), "Vec::operator=(Mat<Num_T>): wrong size");
01528 
01529     if (m.cols() == 1) {
01530       set_size(m.rows(), false);
01531       copy_vector(m.rows(), m._data(), data);
01532     } else if (m.rows() == 1) {
01533       set_size(m.cols(), false);
01534       copy_vector(m.cols(), m._data(), m.rows(), data, 1);
01535     } else
01536       it_error("Vec::operator=(Mat<Num_T>): wrong size");
01537     return *this;
01538   }
01539 
01540   template<class Num_T> inline
01541   Vec<Num_T>& Vec<Num_T>::operator=(const char *values)
01542   {
01543     set(values);
01544     return *this;
01545   }
01546 
01547   template<>
01548   bvec Vec<std::complex<double> >::operator==(std::complex<double>) const;
01549 
01550   template<class Num_T>
01551   bvec Vec<Num_T>::operator==(const Num_T value) const
01552   {
01553     it_assert(datasize > 0, "Vec::operator==: vector must have size > 0");
01554     Vec<Num_T> invector(*this);
01555     bvec temp(invector.length());
01556 
01557     for (int i=0;i<invector.length();i++)
01558       temp(i)=(invector(i)==value);
01559 
01560     return temp;
01561   }
01562 
01563   template<>
01564   bvec Vec<std::complex<double> >::operator!=(std::complex<double>) const;
01565 
01566   template<class Num_T>
01567   bvec Vec<Num_T>::operator!=(const Num_T value) const
01568   {
01569     it_assert(datasize > 0, "Vec::operator!=: vector must have size > 0");
01570     Vec<Num_T> invector(*this);
01571     bvec temp(invector.length());
01572 
01573     for (int i=0;i<invector.length();i++)
01574       temp(i)=(invector(i)!=value);
01575 
01576     return temp;
01577   }
01578 
01579   template<>
01580   bvec Vec<std::complex<double> >::operator<(std::complex<double>) const;
01581 
01582   template<class Num_T>
01583   bvec Vec<Num_T>::operator<(const Num_T value) const
01584   {
01585     it_assert(datasize > 0, "Vec::operator<: vector must have size > 0");
01586     Vec<Num_T> invector(*this);
01587     bvec temp(invector.length());
01588 
01589     for (int i=0;i<invector.length();i++)
01590       temp(i)=(invector(i)<value);
01591 
01592     return temp;
01593   }
01594 
01595   template<>
01596   bvec Vec<std::complex<double> >::operator<=(std::complex<double>) const;
01597 
01598   template<class Num_T>
01599   bvec Vec<Num_T>::operator<=(const Num_T value) const
01600   {
01601     it_assert(datasize > 0, "Vec::operator<=: vector must have size > 0");
01602     Vec<Num_T> invector(*this);
01603     bvec temp(invector.length());
01604 
01605     for (int i=0;i<invector.length();i++)
01606       temp(i)=(invector(i)<=value);
01607 
01608     return temp;
01609   }
01610 
01611   template<>
01612   bvec Vec<std::complex<double> >::operator>(std::complex<double>) const;
01613 
01614   template<class Num_T>
01615   bvec Vec<Num_T>::operator>(const Num_T value) const
01616   {
01617     it_assert(datasize > 0, "Vec::operator>: vector must have size > 0");
01618     Vec<Num_T> invector(*this);
01619     bvec temp(invector.length());
01620 
01621     for (int i=0;i<invector.length();i++)
01622       temp(i)=(invector(i)>value);
01623 
01624     return temp;
01625   }
01626 
01627   template<>
01628   bvec Vec<std::complex<double> >::operator>=(std::complex<double>) const;
01629 
01630   template<class Num_T>
01631   bvec Vec<Num_T>::operator>=(const Num_T value) const
01632   {
01633     it_assert(datasize > 0, "Vec::operator>=: vector must have size > 0");
01634     Vec<Num_T> invector(*this);
01635     bvec temp(invector.length());
01636 
01637     for (int i=0;i<invector.length();i++)
01638       temp(i)=(invector(i)>=value);
01639 
01640     return temp;
01641   }
01642 
01643   template<class Num_T>
01644   bool Vec<Num_T>::operator==(const Vec<Num_T> &invector) const
01645   {
01646     // OBS ! if wrong size, return false
01647     if (datasize!=invector.datasize) return false;
01648     for (int i=0;i<datasize;i++) {
01649       if (data[i]!=invector.data[i]) return false;
01650     }
01651     return true;
01652   }
01653 
01654   template<class Num_T>
01655   bool Vec<Num_T>::operator!=(const Vec<Num_T> &invector) const
01656   {
01657     if (datasize!=invector.datasize) return true;
01658     for (int i=0;i<datasize;i++) {
01659       if (data[i]!=invector.data[i]) return true;
01660     }
01661     return false;
01662   }
01663 
01665   template<class Num_T>
01666   std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v)
01667   {
01668     int i, sz=v.length();
01669 
01670     os << "[" ;
01671     for (i=0; i<sz; i++) {
01672       os << v(i) ;
01673       if (i < sz-1)
01674         os << " ";
01675     }
01676     os << "]" ;
01677 
01678     return os;
01679   }
01680 
01682   template<class Num_T>
01683   std::istream &operator>>(std::istream &is, Vec<Num_T> &v)
01684   {
01685     std::ostringstream buffer;
01686     bool started = false;
01687     bool finished = false;
01688     bool brackets = false;
01689     char c;
01690 
01691     while (!finished) {
01692       if (is.eof()) {
01693         finished = true;
01694       } else {
01695         c = is.get();
01696 
01697         if (is.eof() || (c == '\n')) {
01698           if (brackets) {
01699             // Right bracket missing
01700             is.setstate(std::ios_base::failbit);
01701             finished = true;
01702           } else if (!((c == '\n') && !started)) {
01703             finished = true;
01704           }
01705         } else if ((c == ' ') || (c == '\t')) {
01706           if (started) {
01707             buffer << ' ';
01708           }
01709         } else if (c == '[') {
01710           if (started) {
01711             // Unexpected left bracket
01712             is.setstate(std::ios_base::failbit);
01713             finished = true;
01714           } else {
01715             started = true;
01716             brackets = true;
01717           }
01718         } else if (c == ']') {
01719           if (!started || !brackets) {
01720             // Unexpected right bracket
01721             is.setstate(std::ios_base::failbit);
01722             finished = true;
01723           } else {
01724             finished = true;
01725           }
01726           while (!is.eof() && (((c = is.peek()) == ' ') || (c == '\t'))) {
01727             is.get();
01728           }
01729           if (!is.eof() && (c == '\n')) {
01730             is.get();
01731           }
01732         } else {
01733           started = true;
01734           buffer << c;
01735         }
01736       }
01737     }
01738 
01739     if (!started) {
01740       v.set_size(0, false);
01741     } else {
01742       v.set(buffer.str());
01743     }
01744 
01745     return is;
01746   }
01747 
01749 
01750   // ----------------------------------------------------------------------
01751   // Instantiations
01752   // ----------------------------------------------------------------------
01753 
01754 #ifdef HAVE_EXTERN_TEMPLATE
01755 
01756   extern template class Vec<double>;
01757   extern template class Vec<int>;
01758   extern template class Vec<short int>;
01759   extern template class Vec<std::complex<double> >;
01760   extern template class Vec<bin>;
01761 
01762   // addition operator
01763 
01764   extern template const vec operator+(const vec &v1, const vec &v2);
01765   extern template const cvec operator+(const cvec &v1, const cvec &v2);
01766   extern template const ivec operator+(const ivec &v1, const ivec &v2);
01767   extern template const svec operator+(const svec &v1, const svec &v2);
01768   extern template const bvec operator+(const bvec &v1, const bvec &v2);
01769 
01770   extern template const vec operator+(const vec &v1, double t);
01771   extern template const cvec operator+(const cvec &v1, std::complex<double> t);
01772   extern template const ivec operator+(const ivec &v1, int t);
01773   extern template const svec operator+(const svec &v1, short t);
01774   extern template const bvec operator+(const bvec &v1, bin t);
01775 
01776   extern template const vec operator+(double t, const vec &v1);
01777   extern template const cvec operator+(std::complex<double> t, const cvec &v1);
01778   extern template const ivec operator+(int t, const ivec &v1);
01779   extern template const svec operator+(short t, const svec &v1);
01780   extern template const bvec operator+(bin t, const bvec &v1);
01781 
01782   // subraction operator
01783 
01784   extern template const vec operator-(const vec &v1, const vec &v2);
01785   extern template const cvec operator-(const cvec &v1, const cvec &v2);
01786   extern template const ivec operator-(const ivec &v1, const ivec &v2);
01787   extern template const svec operator-(const svec &v1, const svec &v2);
01788   extern template const bvec operator-(const bvec &v1, const bvec &v2);
01789 
01790   extern template const vec operator-(const vec &v, double t);
01791   extern template const cvec operator-(const cvec &v, std::complex<double> t);
01792   extern template const ivec operator-(const ivec &v, int t);
01793   extern template const svec operator-(const svec &v, short t);
01794   extern template const bvec operator-(const bvec &v, bin t);
01795 
01796   extern template const vec operator-(double t, const vec &v);
01797   extern template const cvec operator-(std::complex<double> t, const cvec &v);
01798   extern template const ivec operator-(int t, const ivec &v);
01799   extern template const svec operator-(short t, const svec &v);
01800   extern template const bvec operator-(bin t, const bvec &v);
01801 
01802   // unary minus
01803 
01804   extern template const vec operator-(const vec &v);
01805   extern template const cvec operator-(const cvec &v);
01806   extern template const ivec operator-(const ivec &v);
01807   extern template const svec operator-(const svec &v);
01808   extern template const bvec operator-(const bvec &v);
01809 
01810   // multiplication operator
01811 
01812 #if !defined(HAVE_BLAS)
01813   extern template double dot(const vec &v1, const vec &v2);
01814   extern template std::complex<double> dot(const cvec &v1, const cvec &v2);
01815 #endif
01816   extern template int dot(const ivec &v1, const ivec &v2);
01817   extern template short dot(const svec &v1, const svec &v2);
01818   extern template bin dot(const bvec &v1, const bvec &v2);
01819 
01820 #if !defined(HAVE_BLAS)
01821   extern template double operator*(const vec &v1, const vec &v2);
01822   extern template std::complex<double> operator*(const cvec &v1,
01823                                                  const cvec &v2);
01824 #endif
01825   extern template int operator*(const ivec &v1, const ivec &v2);
01826   extern template short operator*(const svec &v1, const svec &v2);
01827   extern template bin operator*(const bvec &v1, const bvec &v2);
01828 
01829 #if !defined(HAVE_BLAS)
01830   extern template const mat outer_product(const vec &v1, const vec &v2,
01831                                           bool hermitian);
01832 #endif
01833   extern template const imat outer_product(const ivec &v1, const ivec &v2,
01834                                            bool hermitian);
01835   extern template const smat outer_product(const svec &v1, const svec &v2,
01836                                            bool hermitian);
01837   extern template const bmat outer_product(const bvec &v1, const bvec &v2,
01838                                            bool hermitian);
01839 
01840   extern template const vec operator*(const vec &v, double t);
01841   extern template const cvec operator*(const cvec &v, std::complex<double> t);
01842   extern template const ivec operator*(const ivec &v, int t);
01843   extern template const svec operator*(const svec &v, short t);
01844   extern template const bvec operator*(const bvec &v, bin t);
01845 
01846   extern template const vec operator*(double t, const vec &v);
01847   extern template const cvec operator*(std::complex<double> t, const cvec &v);
01848   extern template const ivec operator*(int t, const ivec &v);
01849   extern template const svec operator*(short t, const svec &v);
01850   extern template const bvec operator*(bin t, const bvec &v);
01851 
01852   // elementwise multiplication
01853 
01854   extern template const vec elem_mult(const vec &a, const vec &b);
01855   extern template const cvec elem_mult(const cvec &a, const cvec &b);
01856   extern template const ivec elem_mult(const ivec &a, const ivec &b);
01857   extern template const svec elem_mult(const svec &a, const svec &b);
01858   extern template const bvec elem_mult(const bvec &a, const bvec &b);
01859 
01860   extern template void elem_mult_out(const vec &a, const vec &b, vec &out);
01861   extern template void elem_mult_out(const cvec &a, const cvec &b, cvec &out);
01862   extern template void elem_mult_out(const ivec &a, const ivec &b, ivec &out);
01863   extern template void elem_mult_out(const svec &a, const svec &b, svec &out);
01864   extern template void elem_mult_out(const bvec &a, const bvec &b, bvec &out);
01865 
01866   extern template const vec elem_mult(const vec &a, const vec &b,
01867                                       const vec &c);
01868   extern template const cvec elem_mult(const cvec &a, const cvec &b,
01869                                        const cvec &c);
01870   extern template const ivec elem_mult(const ivec &a, const ivec &b,
01871                                        const ivec &c);
01872   extern template const svec elem_mult(const svec &a, const svec &b,
01873                                        const svec &c);
01874   extern template const bvec elem_mult(const bvec &a, const bvec &b,
01875                                        const bvec &c);
01876 
01877   extern template void elem_mult_out(const vec &a, const vec &b, const vec &c,
01878                                      vec &out);
01879   extern template void elem_mult_out(const cvec &a, const cvec &b,
01880                                      const cvec &c, cvec &out);
01881   extern template void elem_mult_out(const ivec &a, const ivec &b,
01882                                      const ivec &c, ivec &out);
01883   extern template void elem_mult_out(const svec &a, const svec &b,
01884                                      const svec &c, svec &out);
01885   extern template void elem_mult_out(const bvec &a, const bvec &b,
01886                                      const bvec &c, bvec &out);
01887 
01888   extern template const vec elem_mult(const vec &a, const vec &b, const vec &c,
01889                                       const vec &d);
01890   extern template const cvec elem_mult(const cvec &a, const cvec &b,
01891                                        const cvec &c, const cvec &d);
01892   extern template const ivec elem_mult(const ivec &a, const ivec &b,
01893                                        const ivec &c, const ivec &d);
01894   extern template const svec elem_mult(const svec &a, const svec &b,
01895                                        const svec &c, const svec &d);
01896   extern template const bvec elem_mult(const bvec &a, const bvec &b,
01897                                        const bvec &c, const bvec &d);
01898 
01899   extern template void elem_mult_out(const vec &a, const vec &b, const vec &c,
01900                                      const vec &d, vec &out);
01901   extern template void elem_mult_out(const cvec &a, const cvec &b,
01902                                      const cvec &c, const cvec &d, cvec &out);
01903   extern template void elem_mult_out(const ivec &a, const ivec &b,
01904                                      const ivec &c, const ivec &d, ivec &out);
01905   extern template void elem_mult_out(const svec &a, const svec &b,
01906                                      const svec &c, const svec &d, svec &out);
01907   extern template void elem_mult_out(const bvec &a, const bvec &b,
01908                                      const bvec &c, const bvec &d, bvec &out);
01909 
01910   // in-place elementwise multiplication
01911 
01912   extern template void elem_mult_inplace(const vec &a, vec &b);
01913   extern template void elem_mult_inplace(const cvec &a, cvec &b);
01914   extern template void elem_mult_inplace(const ivec &a, ivec &b);
01915   extern template void elem_mult_inplace(const svec &a, svec &b);
01916   extern template void elem_mult_inplace(const bvec &a, bvec &b);
01917 
01918   // elementwise multiplication followed by summation
01919 
01920   extern template double elem_mult_sum(const vec &a, const vec &b);
01921   extern template std::complex<double> elem_mult_sum(const cvec &a,
01922                                                      const cvec &b);
01923   extern template int elem_mult_sum(const ivec &a, const ivec &b);
01924   extern template short elem_mult_sum(const svec &a, const svec &b);
01925   extern template bin elem_mult_sum(const bvec &a, const bvec &b);
01926 
01927   // division operator
01928 
01929   extern template const vec operator/(const vec &v, double t);
01930   extern template const cvec operator/(const cvec &v, std::complex<double> t);
01931   extern template const ivec operator/(const ivec &v, int t);
01932   extern template const svec operator/(const svec &v, short t);
01933   extern template const bvec operator/(const bvec &v, bin t);
01934 
01935   extern template const vec operator/(double t, const vec &v);
01936   extern template const cvec operator/(std::complex<double> t, const cvec &v);
01937   extern template const ivec operator/(int t, const ivec &v);
01938   extern template const svec operator/(short t, const svec &v);
01939   extern template const bvec operator/(bin t, const bvec &v);
01940 
01941   // elementwise division operator
01942 
01943   extern template const vec elem_div(const vec &a, const vec &b);
01944   extern template const cvec elem_div(const cvec &a, const cvec &b);
01945   extern template const ivec elem_div(const ivec &a, const ivec &b);
01946   extern template const svec elem_div(const svec &a, const svec &b);
01947   extern template const bvec elem_div(const bvec &a, const bvec &b);
01948 
01949   extern template const vec elem_div(double t, const vec &v);
01950   extern template const cvec elem_div(std::complex<double> t, const cvec &v);
01951   extern template const ivec elem_div(int t, const ivec &v);
01952   extern template const svec elem_div(short t, const svec &v);
01953   extern template const bvec elem_div(bin t, const bvec &v);
01954 
01955   extern template void elem_div_out(const vec &a, const vec &b, vec &out);
01956   extern template void elem_div_out(const cvec &a, const cvec &b, cvec &out);
01957   extern template void elem_div_out(const ivec &a, const ivec &b, ivec &out);
01958   extern template void elem_div_out(const svec &a, const svec &b, svec &out);
01959   extern template void elem_div_out(const bvec &a, const bvec &b, bvec &out);
01960 
01961   // elementwise division followed by summation
01962 
01963   extern template double elem_div_sum(const vec &a, const vec &b);
01964   extern template std::complex<double> elem_div_sum(const cvec &a,
01965                                                     const cvec &b);
01966   extern template int elem_div_sum(const ivec &a, const ivec &b);
01967   extern template short elem_div_sum(const svec &a, const svec &b);
01968   extern template bin elem_div_sum(const bvec &a, const bvec &b);
01969 
01970   // concat operator
01971 
01972   extern template const vec concat(const vec &v, double a);
01973   extern template const cvec concat(const cvec &v, std::complex<double> a);
01974   extern template const ivec concat(const ivec &v, int a);
01975   extern template const svec concat(const svec &v, short a);
01976   extern template const bvec concat(const bvec &v, bin a);
01977 
01978   extern template const vec concat(double a, const vec &v);
01979   extern template const cvec concat(std::complex<double> a, const cvec &v);
01980   extern template const ivec concat(int a, const ivec &v);
01981   extern template const svec concat(short a, const svec &v);
01982   extern template const bvec concat(bin a, const bvec &v);
01983 
01984   extern template const vec concat(const vec &v1, const vec &v2);
01985   extern template const cvec concat(const cvec &v1, const cvec &v2);
01986   extern template const ivec concat(const ivec &v1, const ivec &v2);
01987   extern template const svec concat(const svec &v1, const svec &v2);
01988   extern template const bvec concat(const bvec &v1, const bvec &v2);
01989 
01990   extern template const vec concat(const vec &v1, const vec &v2,
01991                                    const vec &v3);
01992   extern template const cvec concat(const cvec &v1, const cvec &v2,
01993                                     const cvec &v3);
01994   extern template const ivec concat(const ivec &v1, const ivec &v2,
01995                                     const ivec &v3);
01996   extern template const svec concat(const svec &v1, const svec &v2,
01997                                     const svec &v3);
01998   extern template const bvec concat(const bvec &v1, const bvec &v2,
01999                                     const bvec &v3);
02000 
02001   extern template const vec concat(const vec &v1, const vec &v2, const vec &v3,
02002                                    const vec &v4);
02003   extern template const cvec concat(const cvec &v1, const cvec &v2,
02004                                     const cvec &v3, const cvec &v4);
02005   extern template const ivec concat(const ivec &v1, const ivec &v2,
02006                                     const ivec &v3, const ivec &v4);
02007   extern template const svec concat(const svec &v1, const svec &v2,
02008                                     const svec &v3, const svec &v4);
02009   extern template const bvec concat(const bvec &v1, const bvec &v2,
02010                                     const bvec &v3, const bvec &v4);
02011 
02012   extern template const vec concat(const vec &v1, const vec &v2, const vec &v3,
02013                                    const vec &v4, const vec &v5);
02014   extern template const cvec concat(const cvec &v1, const cvec &v2,
02015                                     const cvec &v3, const cvec &v4,
02016                                     const cvec &v5);
02017   extern template const ivec concat(const ivec &v1, const ivec &v2,
02018                                     const ivec &v3, const ivec &v4,
02019                                     const ivec &v5);
02020   extern template const svec concat(const svec &v1, const svec &v2,
02021                                     const svec &v3, const svec &v4,
02022                                     const svec &v5);
02023   extern template const bvec concat(const bvec &v1, const bvec &v2,
02024                                     const bvec &v3, const bvec &v4,
02025                                     const bvec &v5);
02026 
02027   // I/O streams
02028 
02029   extern template std::ostream &operator<<(std::ostream& os, const vec &vect);
02030   extern template std::ostream &operator<<(std::ostream& os, const cvec &vect);
02031   extern template std::ostream &operator<<(std::ostream& os, const svec &vect);
02032   extern template std::ostream &operator<<(std::ostream& os, const ivec &vect);
02033   extern template std::ostream &operator<<(std::ostream& os, const bvec &vect);
02034   extern template std::istream &operator>>(std::istream& is, vec &vect);
02035   extern template std::istream &operator>>(std::istream& is, cvec &vect);
02036   extern template std::istream &operator>>(std::istream& is, svec &vect);
02037   extern template std::istream &operator>>(std::istream& is, ivec &vect);
02038   extern template std::istream &operator>>(std::istream& is, bvec &vect);
02039 
02040 #endif // HAVE_EXTERN_TEMPLATE
02041 
02043 
02044 } // namespace itpp
02045 
02046 #endif // #ifndef VEC_H
SourceForge Logo

Generated on Mon Jan 7 22:28:57 2008 for IT++ by Doxygen 1.5.4