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