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