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 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
Generated on Sun Jul 26 08:54:55 2009 for IT++ by Doxygen 1.5.9