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
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 
01609 template<class Num_T>
01610 bvec Vec<Num_T>::operator==(Num_T t) const
01611 {
01612   it_assert_debug(datasize > 0, "Vec<>::operator==(): Wrong size");
01613   bvec temp(datasize);
01614   for (int i = 0; i < datasize; i++)
01615     temp(i) = (data[i] == t);
01616   return temp;
01617 }
01618 
01619 template<class Num_T>
01620 bvec Vec<Num_T>::operator!=(Num_T t) const
01621 {
01622   it_assert_debug(datasize > 0, "Vec<>::operator!=(): Wrong size");
01623   bvec temp(datasize);
01624   for (int i = 0; i < datasize; i++)
01625     temp(i) = (data[i] != t);
01626   return temp;
01627 }
01628 
01630 template<>
01631 bvec Vec<std::complex<double> >::operator<(std::complex<double>) const;
01633 
01634 template<class Num_T>
01635 bvec Vec<Num_T>::operator<(Num_T t) const
01636 {
01637   it_assert_debug(datasize > 0, "Vec<>::operator<(): Wrong size");
01638   bvec temp(datasize);
01639   for (int i = 0; i < datasize; i++)
01640     temp(i) = (data[i] < t);
01641   return temp;
01642 }
01643 
01645 template<>
01646 bvec Vec<std::complex<double> >::operator<=(std::complex<double>) const;
01648 
01649 template<class Num_T>
01650 bvec Vec<Num_T>::operator<=(Num_T t) const
01651 {
01652   it_assert_debug(datasize > 0, "Vec<>::operator<=(): Wrong size");
01653   bvec temp(datasize);
01654   for (int i = 0; i < datasize; i++)
01655     temp(i) = (data[i] <= t);
01656   return temp;
01657 }
01658 
01660 template<>
01661 bvec Vec<std::complex<double> >::operator>(std::complex<double>) const;
01663 
01664 template<class Num_T>
01665 bvec Vec<Num_T>::operator>(Num_T t) const
01666 {
01667   it_assert_debug(datasize > 0, "Vec<>::operator>(): Wrong size");
01668   bvec temp(datasize);
01669   for (int i = 0; i < datasize; i++)
01670     temp(i) = (data[i] > t);
01671   return temp;
01672 }
01673 
01675 template<>
01676 bvec Vec<std::complex<double> >::operator>=(std::complex<double>) const;
01678 
01679 template<class Num_T>
01680 bvec Vec<Num_T>::operator>=(Num_T t) const
01681 {
01682   it_assert_debug(datasize > 0, "Vec<>::operator>=(): Wrong size");
01683   bvec temp(datasize);
01684   for (int i = 0; i < datasize; i++)
01685     temp(i) = (data[i] >= t);
01686   return temp;
01687 }
01688 
01689 template<class Num_T>
01690 bool Vec<Num_T>::operator==(const Vec<Num_T> &invector) const
01691 {
01692   // OBS ! if wrong size, return false
01693   if (datasize != invector.datasize) return false;
01694   for (int i = 0;i < datasize;i++) {
01695     if (data[i] != invector.data[i]) return false;
01696   }
01697   return true;
01698 }
01699 
01700 template<class Num_T>
01701 bool Vec<Num_T>::operator!=(const Vec<Num_T> &invector) const
01702 {
01703   if (datasize != invector.datasize) return true;
01704   for (int i = 0;i < datasize;i++) {
01705     if (data[i] != invector.data[i]) return true;
01706   }
01707   return false;
01708 }
01709 
01711 template<class Num_T>
01712 std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v)
01713 {
01714   int i, sz = v.length();
01715 
01716   os << "[" ;
01717   for (i = 0; i < sz; i++) {
01718     os << v(i) ;
01719     if (i < sz - 1)
01720       os << " ";
01721   }
01722   os << "]" ;
01723 
01724   return os;
01725 }
01726 
01728 template<class Num_T>
01729 std::istream &operator>>(std::istream &is, Vec<Num_T> &v)
01730 {
01731   std::ostringstream buffer;
01732   bool started = false;
01733   bool finished = false;
01734   bool brackets = false;
01735   char c;
01736 
01737   while (!finished) {
01738     if (is.eof()) {
01739       finished = true;
01740     }
01741     else {
01742       is.get(c);
01743 
01744       if (is.eof() || (c == '\n')) {
01745         if (brackets) {
01746           // Right bracket missing
01747           is.setstate(std::ios_base::failbit);
01748           finished = true;
01749         }
01750         else if (!((c == '\n') && !started)) {
01751           finished = true;
01752         }
01753       }
01754       else if ((c == ' ') || (c == '\t')) {
01755         if (started) {
01756           buffer << ' ';
01757         }
01758       }
01759       else if (c == '[') {
01760         if (started) {
01761           // Unexpected left bracket
01762           is.setstate(std::ios_base::failbit);
01763           finished = true;
01764         }
01765         else {
01766           started = true;
01767           brackets = true;
01768         }
01769       }
01770       else if (c == ']') {
01771         if (!started || !brackets) {
01772           // Unexpected right bracket
01773           is.setstate(std::ios_base::failbit);
01774           finished = true;
01775         }
01776         else {
01777           finished = true;
01778         }
01779         while (!is.eof() && (((c = static_cast<char>(is.peek())) == ' ')
01780                              || (c == '\t'))) {
01781           is.get();
01782         }
01783         if (!is.eof() && (c == '\n')) {
01784           is.get();
01785         }
01786       }
01787       else {
01788         started = true;
01789         buffer << c;
01790       }
01791     }
01792   }
01793 
01794   if (!started) {
01795     v.set_size(0, false);
01796   }
01797   else {
01798     v.set(buffer.str());
01799   }
01800 
01801   return is;
01802 }
01803 
01805 
01806 // ----------------------------------------------------------------------
01807 // Instantiations
01808 // ----------------------------------------------------------------------
01809 
01810 #ifdef HAVE_EXTERN_TEMPLATE
01811 
01812 extern template class Vec<double>;
01813 extern template class Vec<int>;
01814 extern template class Vec<short int>;
01815 extern template class Vec<std::complex<double> >;
01816 extern template class Vec<bin>;
01817 
01818 // addition operator
01819 
01820 extern template vec operator+(const vec &v1, const vec &v2);
01821 extern template cvec operator+(const cvec &v1, const cvec &v2);
01822 extern template ivec operator+(const ivec &v1, const ivec &v2);
01823 extern template svec operator+(const svec &v1, const svec &v2);
01824 extern template bvec operator+(const bvec &v1, const bvec &v2);
01825 
01826 extern template vec operator+(const vec &v1, double t);
01827 extern template cvec operator+(const cvec &v1, std::complex<double> t);
01828 extern template ivec operator+(const ivec &v1, int t);
01829 extern template svec operator+(const svec &v1, short t);
01830 extern template bvec operator+(const bvec &v1, bin t);
01831 
01832 extern template vec operator+(double t, const vec &v1);
01833 extern template cvec operator+(std::complex<double> t, const cvec &v1);
01834 extern template ivec operator+(int t, const ivec &v1);
01835 extern template svec operator+(short t, const svec &v1);
01836 extern template bvec operator+(bin t, const bvec &v1);
01837 
01838 // subtraction operator
01839 
01840 extern template vec operator-(const vec &v1, const vec &v2);
01841 extern template cvec operator-(const cvec &v1, const cvec &v2);
01842 extern template ivec operator-(const ivec &v1, const ivec &v2);
01843 extern template svec operator-(const svec &v1, const svec &v2);
01844 extern template bvec operator-(const bvec &v1, const bvec &v2);
01845 
01846 extern template vec operator-(const vec &v, double t);
01847 extern template cvec operator-(const cvec &v, std::complex<double> t);
01848 extern template ivec operator-(const ivec &v, int t);
01849 extern template svec operator-(const svec &v, short t);
01850 extern template bvec operator-(const bvec &v, bin t);
01851 
01852 extern template vec operator-(double t, const vec &v);
01853 extern template cvec operator-(std::complex<double> t, const cvec &v);
01854 extern template ivec operator-(int t, const ivec &v);
01855 extern template svec operator-(short t, const svec &v);
01856 extern template bvec operator-(bin t, const bvec &v);
01857 
01858 // unary minus
01859 
01860 extern template vec operator-(const vec &v);
01861 extern template cvec operator-(const cvec &v);
01862 extern template ivec operator-(const ivec &v);
01863 extern template svec operator-(const svec &v);
01864 extern template bvec operator-(const bvec &v);
01865 
01866 // multiplication operator
01867 
01868 #if !defined(HAVE_BLAS)
01869 extern template double dot(const vec &v1, const vec &v2);
01870 #if !(defined(HAVE_ZDOTUSUB) || defined(HAVE_ZDOTU_VOID))
01871 extern template std::complex<double> dot(const cvec &v1, const cvec &v2);
01872 #endif // !(HAVE_ZDOTUSUB || HAVE_ZDOTU_VOID)
01873 #endif // HAVE_BLAS
01874 extern template int dot(const ivec &v1, const ivec &v2);
01875 extern template short dot(const svec &v1, const svec &v2);
01876 extern template bin dot(const bvec &v1, const bvec &v2);
01877 
01878 #if !defined(HAVE_BLAS)
01879 extern template double operator*(const vec &v1, const vec &v2);
01880 extern template std::complex<double> operator*(const cvec &v1,
01881     const cvec &v2);
01882 #endif
01883 extern template int operator*(const ivec &v1, const ivec &v2);
01884 extern template short operator*(const svec &v1, const svec &v2);
01885 extern template bin operator*(const bvec &v1, const bvec &v2);
01886 
01887 #if !defined(HAVE_BLAS)
01888 extern template mat outer_product(const vec &v1, const vec &v2,
01889                                     bool hermitian);
01890 #endif
01891 extern template imat outer_product(const ivec &v1, const ivec &v2,
01892                                      bool hermitian);
01893 extern template smat outer_product(const svec &v1, const svec &v2,
01894                                      bool hermitian);
01895 extern template bmat outer_product(const bvec &v1, const bvec &v2,
01896                                      bool hermitian);
01897 
01898 extern template vec operator*(const vec &v, double t);
01899 extern template cvec operator*(const cvec &v, std::complex<double> t);
01900 extern template ivec operator*(const ivec &v, int t);
01901 extern template svec operator*(const svec &v, short t);
01902 extern template bvec operator*(const bvec &v, bin t);
01903 
01904 extern template vec operator*(double t, const vec &v);
01905 extern template cvec operator*(std::complex<double> t, const cvec &v);
01906 extern template ivec operator*(int t, const ivec &v);
01907 extern template svec operator*(short t, const svec &v);
01908 extern template bvec operator*(bin t, const bvec &v);
01909 
01910 // elementwise multiplication
01911 
01912 extern template vec elem_mult(const vec &a, const vec &b);
01913 extern template cvec elem_mult(const cvec &a, const cvec &b);
01914 extern template ivec elem_mult(const ivec &a, const ivec &b);
01915 extern template svec elem_mult(const svec &a, const svec &b);
01916 extern template bvec elem_mult(const bvec &a, const bvec &b);
01917 
01918 extern template void elem_mult_out(const vec &a, const vec &b, vec &out);
01919 extern template void elem_mult_out(const cvec &a, const cvec &b, cvec &out);
01920 extern template void elem_mult_out(const ivec &a, const ivec &b, ivec &out);
01921 extern template void elem_mult_out(const svec &a, const svec &b, svec &out);
01922 extern template void elem_mult_out(const bvec &a, const bvec &b, bvec &out);
01923 
01924 extern template vec elem_mult(const vec &a, const vec &b, const vec &c);
01925 extern template cvec elem_mult(const cvec &a, const cvec &b, const cvec &c);
01926 extern template ivec elem_mult(const ivec &a, const ivec &b, const ivec &c);
01927 extern template svec elem_mult(const svec &a, const svec &b, const svec &c);
01928 extern template bvec elem_mult(const bvec &a, const bvec &b, const bvec &c);
01929 
01930 extern template void elem_mult_out(const vec &a, const vec &b,
01931                                      const vec &c, vec &out);
01932 extern template void elem_mult_out(const cvec &a, const cvec &b,
01933                                      const cvec &c, cvec &out);
01934 extern template void elem_mult_out(const ivec &a, const ivec &b,
01935                                      const ivec &c, ivec &out);
01936 extern template void elem_mult_out(const svec &a, const svec &b,
01937                                      const svec &c, svec &out);
01938 extern template void elem_mult_out(const bvec &a, const bvec &b,
01939                                      const bvec &c, bvec &out);
01940 
01941 extern template vec elem_mult(const vec &a, const vec &b,
01942                                 const vec &c, const vec &d);
01943 extern template cvec elem_mult(const cvec &a, const cvec &b,
01944                                  const cvec &c, const cvec &d);
01945 extern template ivec elem_mult(const ivec &a, const ivec &b,
01946                                  const ivec &c, const ivec &d);
01947 extern template svec elem_mult(const svec &a, const svec &b,
01948                                  const svec &c, const svec &d);
01949 extern template bvec elem_mult(const bvec &a, const bvec &b,
01950                                  const bvec &c, const bvec &d);
01951 
01952 extern template void elem_mult_out(const vec &a, const vec &b, const vec &c,
01953                                      const vec &d, vec &out);
01954 extern template void elem_mult_out(const cvec &a, const cvec &b,
01955                                      const cvec &c, const cvec &d, cvec &out);
01956 extern template void elem_mult_out(const ivec &a, const ivec &b,
01957                                      const ivec &c, const ivec &d, ivec &out);
01958 extern template void elem_mult_out(const svec &a, const svec &b,
01959                                      const svec &c, const svec &d, svec &out);
01960 extern template void elem_mult_out(const bvec &a, const bvec &b,
01961                                      const bvec &c, const bvec &d, bvec &out);
01962 
01963 // in-place elementwise multiplication
01964 
01965 extern template void elem_mult_inplace(const vec &a, vec &b);
01966 extern template void elem_mult_inplace(const cvec &a, cvec &b);
01967 extern template void elem_mult_inplace(const ivec &a, ivec &b);
01968 extern template void elem_mult_inplace(const svec &a, svec &b);
01969 extern template void elem_mult_inplace(const bvec &a, bvec &b);
01970 
01971 // elementwise multiplication followed by summation
01972 
01973 extern template double elem_mult_sum(const vec &a, const vec &b);
01974 extern template std::complex<double> elem_mult_sum(const cvec &a,
01975     const cvec &b);
01976 extern template int elem_mult_sum(const ivec &a, const ivec &b);
01977 extern template short elem_mult_sum(const svec &a, const svec &b);
01978 extern template bin elem_mult_sum(const bvec &a, const bvec &b);
01979 
01980 // division operator
01981 
01982 extern template vec operator/(const vec &v, double t);
01983 extern template cvec operator/(const cvec &v, std::complex<double> t);
01984 extern template ivec operator/(const ivec &v, int t);
01985 extern template svec operator/(const svec &v, short t);
01986 extern template bvec operator/(const bvec &v, bin t);
01987 
01988 extern template vec operator/(double t, const vec &v);
01989 extern template cvec operator/(std::complex<double> t, const cvec &v);
01990 extern template ivec operator/(int t, const ivec &v);
01991 extern template svec operator/(short t, const svec &v);
01992 extern template bvec operator/(bin t, const bvec &v);
01993 
01994 // elementwise division operator
01995 
01996 extern template vec elem_div(const vec &a, const vec &b);
01997 extern template cvec elem_div(const cvec &a, const cvec &b);
01998 extern template ivec elem_div(const ivec &a, const ivec &b);
01999 extern template svec elem_div(const svec &a, const svec &b);
02000 extern template bvec elem_div(const bvec &a, const bvec &b);
02001 
02002 extern template vec elem_div(double t, const vec &v);
02003 extern template cvec elem_div(std::complex<double> t, const cvec &v);
02004 extern template ivec elem_div(int t, const ivec &v);
02005 extern template svec elem_div(short t, const svec &v);
02006 extern template bvec elem_div(bin t, const bvec &v);
02007 
02008 extern template void elem_div_out(const vec &a, const vec &b, vec &out);
02009 extern template void elem_div_out(const cvec &a, const cvec &b, cvec &out);
02010 extern template void elem_div_out(const ivec &a, const ivec &b, ivec &out);
02011 extern template void elem_div_out(const svec &a, const svec &b, svec &out);
02012 extern template void elem_div_out(const bvec &a, const bvec &b, bvec &out);
02013 
02014 // elementwise division followed by summation
02015 
02016 extern template double elem_div_sum(const vec &a, const vec &b);
02017 extern template std::complex<double> elem_div_sum(const cvec &a,
02018     const cvec &b);
02019 extern template int elem_div_sum(const ivec &a, const ivec &b);
02020 extern template short elem_div_sum(const svec &a, const svec &b);
02021 extern template bin elem_div_sum(const bvec &a, const bvec &b);
02022 
02023 // concat operator
02024 
02025 extern template vec concat(const vec &v, double a);
02026 extern template cvec concat(const cvec &v, std::complex<double> a);
02027 extern template ivec concat(const ivec &v, int a);
02028 extern template svec concat(const svec &v, short a);
02029 extern template bvec concat(const bvec &v, bin a);
02030 
02031 extern template vec concat(double a, const vec &v);
02032 extern template cvec concat(std::complex<double> a, const cvec &v);
02033 extern template ivec concat(int a, const ivec &v);
02034 extern template svec concat(short a, const svec &v);
02035 extern template bvec concat(bin a, const bvec &v);
02036 
02037 extern template vec concat(const vec &v1, const vec &v2);
02038 extern template cvec concat(const cvec &v1, const cvec &v2);
02039 extern template ivec concat(const ivec &v1, const ivec &v2);
02040 extern template svec concat(const svec &v1, const svec &v2);
02041 extern template bvec concat(const bvec &v1, const bvec &v2);
02042 
02043 extern template vec concat(const vec &v1, const vec &v2, const vec &v3);
02044 extern template cvec concat(const cvec &v1, const cvec &v2, const cvec &v3);
02045 extern template ivec concat(const ivec &v1, const ivec &v2, const ivec &v3);
02046 extern template svec concat(const svec &v1, const svec &v2, const svec &v3);
02047 extern template bvec concat(const bvec &v1, const bvec &v2, const bvec &v3);
02048 
02049 extern template vec concat(const vec &v1, const vec &v2,
02050                              const vec &v3, const vec &v4);
02051 extern template cvec concat(const cvec &v1, const cvec &v2,
02052                               const cvec &v3, const cvec &v4);
02053 extern template ivec concat(const ivec &v1, const ivec &v2,
02054                               const ivec &v3, const ivec &v4);
02055 extern template svec concat(const svec &v1, const svec &v2,
02056                               const svec &v3, const svec &v4);
02057 extern template bvec concat(const bvec &v1, const bvec &v2,
02058                               const bvec &v3, const bvec &v4);
02059 
02060 extern template vec concat(const vec &v1, const vec &v2, const vec &v3,
02061                              const vec &v4, const vec &v5);
02062 extern template cvec concat(const cvec &v1, const cvec &v2, const cvec &v3,
02063                               const cvec &v4, const cvec &v5);
02064 extern template ivec concat(const ivec &v1, const ivec &v2, const ivec &v3,
02065                               const ivec &v4, const ivec &v5);
02066 extern template svec concat(const svec &v1, const svec &v2, const svec &v3,
02067                               const svec &v4, const svec &v5);
02068 extern template bvec concat(const bvec &v1, const bvec &v2, const bvec &v3,
02069                               const bvec &v4, const bvec &v5);
02070 
02071 // I/O streams
02072 
02073 extern template std::ostream &operator<<(std::ostream& os, const vec &vect);
02074 extern template std::ostream &operator<<(std::ostream& os, const cvec &vect);
02075 extern template std::ostream &operator<<(std::ostream& os, const svec &vect);
02076 extern template std::ostream &operator<<(std::ostream& os, const ivec &vect);
02077 extern template std::ostream &operator<<(std::ostream& os, const bvec &vect);
02078 extern template std::istream &operator>>(std::istream& is, vec &vect);
02079 extern template std::istream &operator>>(std::istream& is, cvec &vect);
02080 extern template std::istream &operator>>(std::istream& is, svec &vect);
02081 extern template std::istream &operator>>(std::istream& is, ivec &vect);
02082 extern template std::istream &operator>>(std::istream& is, bvec &vect);
02083 
02084 #endif // HAVE_EXTERN_TEMPLATE
02085 
02087 
02088 } // namespace itpp
02089 
02090 #endif // #ifndef VEC_H
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Wed Jan 20 23:03:04 2010 for IT++ by Doxygen 1.6.2