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

Generated on Sat Apr 19 10:41:12 2008 for IT++ by Doxygen 1.5.5