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

Generated on Wed Mar 2 2011 22:05:00 for IT++ by Doxygen 1.7.3