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 #if defined(HAVE_ZDOTU_VOID) 00929 blas::zdotu_(&output, &v1.datasize, v1.data, &incr, v2.data, &incr); 00930 #elif defined(HAVE_ZDOTU_RETURN) 00931 output = blas::zdotu_(&v1.datasize, v1.data, &incr, v2.data, &incr); 00932 #else 00933 blas::zdotusub_(&output, &v1.datasize, v1.data, &incr, v2.data, &incr); 00934 #endif // HAVE_ZDOTU_VOID 00935 return output; 00936 } 00937 #endif // HAVE_BLAS 00938 00939 template<class Num_T> inline 00940 Num_T dot(const Vec<Num_T> &v1, const Vec<Num_T> &v2) 00941 { 00942 int i; 00943 Num_T r=Num_T(0); 00944 00945 it_assert_debug(v1.datasize==v2.datasize, "Vec::dot: wrong sizes"); 00946 for (i=0; i<v1.datasize; i++) 00947 r += v1.data[i] * v2.data[i]; 00948 00949 return r; 00950 } 00951 00952 template<class Num_T> inline 00953 Num_T operator*(const Vec<Num_T> &v1, const Vec<Num_T> &v2) 00954 { 00955 return dot(v1, v2); 00956 } 00957 00958 00959 #if defined(HAVE_BLAS) 00960 template<> inline 00961 const mat outer_product(const vec &v1, const vec &v2, bool hermitian) 00962 { 00963 it_assert_debug((v1.datasize > 0) && (v2.datasize > 0), 00964 "Vec::outer_product():: Input vector of zero size"); 00965 00966 mat out(v1.datasize, v2.datasize); 00967 out.zeros(); 00968 double alpha = 1.0; 00969 int incr = 1; 00970 blas::dger_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr, 00971 v2.data, &incr, out._data(), &v1.datasize); 00972 return out; 00973 } 00974 00975 template<> inline 00976 const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian) 00977 { 00978 it_assert_debug((v1.datasize > 0) && (v2.datasize > 0), 00979 "Vec::outer_product():: Input vector of zero size"); 00980 00981 cmat out(v1.datasize, v2.datasize); 00982 out.zeros(); 00983 std::complex<double> alpha(1.0); 00984 int incr = 1; 00985 if (hermitian) { 00986 blas::zgerc_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr, 00987 v2.data, &incr, out._data(), &v1.datasize); 00988 } 00989 else { 00990 blas::zgeru_(&v1.datasize, &v2.datasize, &alpha, v1.data, &incr, 00991 v2.data, &incr, out._data(), &v1.datasize); 00992 } 00993 return out; 00994 } 00995 #else 00997 template<> inline 00998 const cmat outer_product(const cvec &v1, const cvec &v2, bool hermitian) 00999 { 01000 it_assert_debug((v1.datasize > 0) && (v2.datasize > 0), 01001 "Vec::outer_product():: Input vector of zero size"); 01002 01003 cmat out(v1.datasize, v2.datasize); 01004 if (hermitian) { 01005 for (int i = 0; i < v1.datasize; ++i) { 01006 for (int j = 0; j < v2.datasize; ++j) { 01007 out(i, j) = v1.data[i] * conj(v2.data[j]); 01008 } 01009 } 01010 } 01011 else { 01012 for (int i = 0; i < v1.datasize; ++i) { 01013 for (int j = 0; j < v2.datasize; ++j) { 01014 out(i, j) = v1.data[i] * v2.data[j]; 01015 } 01016 } 01017 } 01018 return out; 01019 } 01020 #endif // HAVE_BLAS 01021 01022 template<class Num_T> inline 01023 const Mat<Num_T> outer_product(const Vec<Num_T> &v1, const Vec<Num_T> &v2, 01024 bool hermitian) 01025 { 01026 int i, j; 01027 01028 it_assert_debug((v1.datasize > 0) && (v2.datasize > 0), 01029 "Vec::outer_product:: Input vector of zero size"); 01030 01031 Mat<Num_T> r(v1.datasize, v2.datasize); 01032 for (i=0; i<v1.datasize; i++) { 01033 for (j=0; j<v2.datasize; j++) { 01034 r(i,j) = v1.data[i] * v2.data[j]; 01035 } 01036 } 01037 01038 return r; 01039 } 01040 01041 template<class Num_T> inline 01042 const Vec<Num_T> operator*(const Vec<Num_T> &v, const Num_T t) 01043 { 01044 int i; 01045 Vec<Num_T> r(v.datasize); 01046 for (i=0; i<v.datasize; i++) 01047 r.data[i] = v.data[i] * t; 01048 01049 return r; 01050 } 01051 01052 template<class Num_T> inline 01053 const Vec<Num_T> operator*(const Num_T t, const Vec<Num_T> &v) 01054 { 01055 return operator*(v, t); 01056 } 01057 01058 template<class Num_T> inline 01059 const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b) 01060 { 01061 Vec<Num_T> out; 01062 elem_mult_out(a,b,out); 01063 return out; 01064 } 01065 01066 template<class Num_T> inline 01067 const Vec<Num_T> elem_mult(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c) 01068 { 01069 Vec<Num_T> out; 01070 elem_mult_out(a,b,c,out); 01071 return out; 01072 } 01073 01074 template<class Num_T> inline 01075 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) 01076 { 01077 Vec<Num_T> out; 01078 elem_mult_out(a,b,c,d,out); 01079 return out; 01080 } 01081 01082 template<class Num_T> inline 01083 void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out) 01084 { 01085 it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_out: wrong sizes"); 01086 01087 if(out.datasize != a.datasize) 01088 out.set_size(a.size()); 01089 01090 for(int i=0; i<a.datasize; i++) 01091 out.data[i] = a.data[i] * b.data[i]; 01092 } 01093 01094 template<class Num_T> inline 01095 void elem_mult_out(const Vec<Num_T> &a, const Vec<Num_T> &b, const Vec<Num_T> &c, Vec<Num_T> &out) 01096 { 01097 it_assert_debug(a.datasize==b.datasize==c.datasize, "Vec::elem_mult_out: wrong sizes"); 01098 01099 if(out.datasize != a.datasize) 01100 out.set_size(a.size()); 01101 01102 for(int i=0; i<a.datasize; i++) 01103 out.data[i] = a.data[i] * b.data[i] * c.data[i]; 01104 } 01105 01106 template<class Num_T> inline 01107 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) 01108 { 01109 it_assert_debug(a.datasize==b.datasize==c.datasize==d.datasize, "Vec::elem_mult_out: wrong sizes"); 01110 01111 if(out.datasize != a.datasize) 01112 out.set_size(a.size()); 01113 01114 for(int i=0; i<a.datasize; i++) 01115 out.data[i] = a.data[i] * b.data[i] * c.data[i] * d.data[i]; 01116 } 01117 01118 template<class Num_T> inline 01119 void elem_mult_inplace(const Vec<Num_T> &a, Vec<Num_T> &b) 01120 { 01121 it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_inplace: wrong sizes"); 01122 01123 for(int i=0; i<a.datasize; i++) 01124 b.data[i] *= a.data[i]; 01125 } 01126 01127 template<class Num_T> inline 01128 Num_T elem_mult_sum(const Vec<Num_T> &a, const Vec<Num_T> &b) 01129 { 01130 it_assert_debug(a.datasize==b.datasize, "Vec::elem_mult_sum: wrong sizes"); 01131 01132 Num_T acc = 0; 01133 01134 for(int i=0; i<a.datasize; i++) 01135 acc += a.data[i] * b.data[i]; 01136 01137 return acc; 01138 } 01139 01140 template<class Num_T> inline 01141 const Vec<Num_T> operator/(const Vec<Num_T> &v, const Num_T t) 01142 { 01143 int i; 01144 Vec<Num_T> r(v.datasize); 01145 01146 for (i=0; i<v.datasize; i++) 01147 r.data[i] = v.data[i] / t; 01148 01149 return r; 01150 } 01151 01152 template<class Num_T> inline 01153 const Vec<Num_T> operator/(const Num_T t, const Vec<Num_T> &v) 01154 { 01155 int i; 01156 Vec<Num_T> r(v.datasize); 01157 01158 for (i=0; i<v.datasize; i++) 01159 r.data[i] = t / v.data[i]; 01160 01161 return r; 01162 } 01163 01164 template<class Num_T> inline 01165 Vec<Num_T>& Vec<Num_T>::operator/=(const Num_T t) 01166 { 01167 for (int i = 0; i < datasize; ++i) { 01168 data[i] /= t; 01169 } 01170 return *this; 01171 } 01172 01173 template<class Num_T> inline 01174 Vec<Num_T>& Vec<Num_T>::operator/=(const Vec<Num_T> &v) 01175 { 01176 it_assert_debug(datasize == v.datasize, "Vec::operator/=(): wrong sizes"); 01177 for (int i = 0; i < datasize; ++i) { 01178 data[i] /= v.data[i]; 01179 } 01180 return *this; 01181 } 01182 01183 template<class Num_T> inline 01184 const Vec<Num_T> elem_div(const Vec<Num_T> &a, const Vec<Num_T> &b) 01185 { 01186 Vec<Num_T> out; 01187 elem_div_out(a,b,out); 01188 return out; 01189 } 01190 01191 template<class Num_T> inline 01192 const Vec<Num_T> elem_div(const Num_T t, const Vec<Num_T> &v) 01193 { 01194 int i; 01195 Vec<Num_T> r(v.datasize); 01196 01197 for (i=0; i<v.datasize; i++) 01198 r.data[i] = t / v.data[i]; 01199 01200 return r; 01201 } 01202 01203 template<class Num_T> inline 01204 void elem_div_out(const Vec<Num_T> &a, const Vec<Num_T> &b, Vec<Num_T> &out) 01205 { 01206 it_assert_debug(a.datasize==b.datasize, "Vecelem_div_out: wrong sizes"); 01207 01208 if(out.datasize != a.datasize) 01209 out.set_size(a.size()); 01210 01211 for(int i=0; i<a.datasize; i++) 01212 out.data[i] = a.data[i] / b.data[i]; 01213 } 01214 01215 template<class Num_T> inline 01216 Num_T elem_div_sum(const Vec<Num_T> &a, const Vec<Num_T> &b) 01217 { 01218 it_assert_debug(a.datasize==b.datasize, "Vec::elem_div_sum: wrong sizes"); 01219 01220 Num_T acc = 0; 01221 01222 for(int i=0; i<a.datasize; i++) 01223 acc += a.data[i] / b.data[i]; 01224 01225 return acc; 01226 } 01227 01228 template<class Num_T> 01229 Vec<Num_T> Vec<Num_T>::get(const Vec<bin> &binlist) const 01230 { 01231 int size = binlist.size(); 01232 it_assert_debug(datasize == size, "Vec::get(bvec &): wrong sizes"); 01233 Vec<Num_T> temp(size); 01234 int j = 0; 01235 for (int i = 0; i < size; ++i) { 01236 if (binlist(i) == bin(1)) { 01237 temp(j) = data[i]; 01238 j++; 01239 } 01240 } 01241 temp.set_size(j, true); 01242 return temp; 01243 } 01244 01245 template<class Num_T> inline 01246 Vec<Num_T> Vec<Num_T>::right(int nr) const 01247 { 01248 it_assert_debug(nr <= datasize, "Vec::right(): index out of range"); 01249 Vec<Num_T> temp(nr); 01250 if (nr > 0) { 01251 copy_vector(nr, &data[datasize-nr], temp.data); 01252 } 01253 return temp; 01254 } 01255 01256 template<class Num_T> inline 01257 Vec<Num_T> Vec<Num_T>::left(int nr) const 01258 { 01259 it_assert_debug(nr <= datasize, "Vec::left(): index out of range"); 01260 Vec<Num_T> temp(nr); 01261 if (nr > 0) { 01262 copy_vector(nr, data, temp.data); 01263 } 01264 return temp; 01265 } 01266 01267 template<class Num_T> inline 01268 Vec<Num_T> Vec<Num_T>::mid(int start, int nr) const 01269 { 01270 it_assert_debug((start >= 0) && ((start+nr) <= datasize), 01271 "Vec::mid(): indexing out of range"); 01272 Vec<Num_T> temp(nr); 01273 if (nr > 0) { 01274 copy_vector(nr, &data[start], temp.data); 01275 } 01276 return temp; 01277 } 01278 01279 template<class Num_T> 01280 Vec<Num_T> Vec<Num_T>::split(int pos) 01281 { 01282 it_assert_debug((pos >= 0) && (pos <= datasize), "Vec::split(): index out of range"); 01283 Vec<Num_T> temp1(pos); 01284 Vec<Num_T> temp2(datasize-pos); 01285 copy_vector(pos, data, temp1.data); 01286 copy_vector(datasize-pos, &data[pos], temp2.data); 01287 (*this) = temp2; 01288 return temp1; 01289 } 01290 01291 template<class Num_T> 01292 void Vec<Num_T>::shift_right(const Num_T In, int n) 01293 { 01294 int i=datasize; 01295 01296 it_assert_debug(n>=0, "Vec::shift_right: index out of range"); 01297 while (--i >= n) 01298 data[i] = data[i-n]; 01299 while (i >= 0) 01300 data[i--] = In; 01301 } 01302 01303 template<class Num_T> 01304 void Vec<Num_T>::shift_right(const Vec<Num_T> &In) 01305 { 01306 int i; 01307 01308 for (i=datasize-1; i>=In.datasize; i--) 01309 data[i]=data[i-In.datasize]; 01310 for (i=0; i<In.datasize; i++) 01311 data[i]=In[i]; 01312 } 01313 01314 template<class Num_T> 01315 void Vec<Num_T>::shift_left(const Num_T In, int n) 01316 { 01317 int i; 01318 01319 it_assert_debug(n>=0, "Vec::shift_left: index out of range"); 01320 for (i=0; i<datasize-n; i++) 01321 data[i] = data[i+n]; 01322 while (i < datasize) 01323 data[i++] = In; 01324 } 01325 01326 template<class Num_T> 01327 void Vec<Num_T>::shift_left(const Vec<Num_T> &In) 01328 { 01329 int i; 01330 01331 for (i=0; i<datasize-In.datasize; i++) 01332 data[i]=data[i+In.datasize]; 01333 for (i=datasize-In.datasize; i<datasize; i++) 01334 data[i]=In[i-datasize+In.datasize]; 01335 } 01336 01337 template<class Num_T> 01338 const Vec<Num_T> concat(const Vec<Num_T> &v, const Num_T a) 01339 { 01340 int size = v.size(); 01341 Vec<Num_T> temp(size+1); 01342 copy_vector(size, v.data, temp.data); 01343 temp(size) = a; 01344 return temp; 01345 } 01346 01347 template<class Num_T> 01348 const Vec<Num_T> concat(const Num_T a, const Vec<Num_T> &v) 01349 { 01350 int size = v.size(); 01351 Vec<Num_T> temp(size+1); 01352 temp(0) = a; 01353 copy_vector(size, v.data, &temp.data[1]); 01354 return temp; 01355 } 01356 01357 template<class Num_T> 01358 const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2) 01359 { 01360 int size1 = v1.size(); 01361 int size2 = v2.size(); 01362 Vec<Num_T> temp(size1+size2); 01363 copy_vector(size1, v1.data, temp.data); 01364 copy_vector(size2, v2.data, &temp.data[size1]); 01365 return temp; 01366 } 01367 01368 template<class Num_T> 01369 const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, 01370 const Vec<Num_T> &v3) 01371 { 01372 int size1 = v1.size(); 01373 int size2 = v2.size(); 01374 int size3 = v3.size(); 01375 Vec<Num_T> temp(size1+size2+size3); 01376 copy_vector(size1, v1.data, temp.data); 01377 copy_vector(size2, v2.data, &temp.data[size1]); 01378 copy_vector(size3, v3.data, &temp.data[size1+size2]); 01379 return temp; 01380 } 01381 01382 template<class Num_T> 01383 const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, 01384 const Vec<Num_T> &v3, const Vec<Num_T> &v4) 01385 { 01386 int size1 = v1.size(); 01387 int size2 = v2.size(); 01388 int size3 = v3.size(); 01389 int size4 = v4.size(); 01390 Vec<Num_T> temp(size1+size2+size3+size4); 01391 copy_vector(size1, v1.data, temp.data); 01392 copy_vector(size2, v2.data, &temp.data[size1]); 01393 copy_vector(size3, v3.data, &temp.data[size1+size2]); 01394 copy_vector(size4, v4.data, &temp.data[size1+size2+size3]); 01395 return temp; 01396 } 01397 01398 template<class Num_T> 01399 const Vec<Num_T> concat(const Vec<Num_T> &v1, const Vec<Num_T> &v2, 01400 const Vec<Num_T> &v3, const Vec<Num_T> &v4, 01401 const Vec<Num_T> &v5) 01402 { 01403 int size1 = v1.size(); 01404 int size2 = v2.size(); 01405 int size3 = v3.size(); 01406 int size4 = v4.size(); 01407 int size5 = v5.size(); 01408 Vec<Num_T> temp(size1+size2+size3+size4+size5); 01409 copy_vector(size1, v1.data, temp.data); 01410 copy_vector(size2, v2.data, &temp.data[size1]); 01411 copy_vector(size3, v3.data, &temp.data[size1+size2]); 01412 copy_vector(size4, v4.data, &temp.data[size1+size2+size3]); 01413 copy_vector(size5, v5.data, &temp.data[size1+size2+size3+size4]); 01414 return temp; 01415 } 01416 01417 template<class Num_T> inline 01418 void Vec<Num_T>::set_subvector(int i1, int i2, const Vec<Num_T> &v) 01419 { 01420 if (i1 == -1) i1 = datasize-1; 01421 if (i2 == -1) i2 = datasize-1; 01422 01423 it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range"); 01424 it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary"); 01425 it_assert_debug(i2-i1+1 == v.datasize, "Vec::set_subvector(): wrong sizes"); 01426 01427 copy_vector(v.datasize, v.data, data+i1); 01428 } 01429 01430 template<class Num_T> inline 01431 void Vec<Num_T>:: set_subvector(int i, const Vec<Num_T> &v) 01432 { 01433 it_assert_debug(i>=0, "Vec::set_subvector(): index out of range"); 01434 it_assert_debug(i+v.datasize <= datasize, "Vec::set_subvector(): too long input vector"); 01435 copy_vector(v.datasize, v.data, data+i); 01436 } 01437 01438 template<class Num_T> 01439 void Vec<Num_T>::set_subvector(int i1, int i2, const Num_T t) 01440 { 01441 if (i1 == -1) i1 = datasize-1; 01442 if (i2 == -1) i2 = datasize-1; 01443 01444 it_assert_debug(i1>=0 && i2>=0 && i1<datasize && i2<datasize, "Vec::set_subvector(): indicies out of range"); 01445 it_assert_debug(i2>=i1, "Vec::set_subvector(): i2 >= i1 necessary"); 01446 01447 for (int i=i1;i<=i2;i++) 01448 data[i] = t; 01449 } 01450 01451 template<class Num_T> 01452 void Vec<Num_T>::replace_mid(int pos, const Vec<Num_T> &v) 01453 { 01454 it_assert_debug((pos>=0) && ((pos+v.length())<=datasize), "Vec::replace_mid: indexing out of range"); 01455 copy_vector(v.datasize, v.data, &data[pos]); 01456 } 01457 01458 template<class Num_T> 01459 void Vec<Num_T>::del(int index) 01460 { 01461 it_assert_debug((index >= 0) && (index < datasize), 01462 "Vec::del(): index out of range"); 01463 Vec<Num_T> temp(*this); 01464 set_size(datasize-1, false); 01465 copy_vector(index, temp.data, data); 01466 copy_vector(datasize-index, &temp.data[index+1], &data[index]); 01467 } 01468 01469 template<class Num_T> 01470 void Vec<Num_T>::del(int i1, int i2) 01471 { 01472 it_assert_debug((i1 >= 0) && (i2 < datasize) && (i1 < i2), 01473 "Vec::del(): index out of range"); 01474 Vec<Num_T> temp(*this); 01475 int new_size = datasize-(i2-i1+1); 01476 set_size(new_size, false); 01477 copy_vector(i1, temp.data, data); 01478 copy_vector(datasize-i1, &temp.data[i2+1], &data[i1]); 01479 } 01480 01481 template<class Num_T> 01482 void Vec<Num_T>::ins(int index, const Num_T in) 01483 { 01484 it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range"); 01485 Vec<Num_T> Temp(*this); 01486 01487 set_size(datasize+1, false); 01488 copy_vector(index, Temp.data, data); 01489 data[index]=in; 01490 copy_vector(Temp.datasize-index, Temp.data+index, data+index+1); 01491 } 01492 01493 template<class Num_T> 01494 void Vec<Num_T>::ins(int index, const Vec<Num_T> &in) 01495 { 01496 it_assert_debug((index>=0) && (index<=datasize), "Vec::ins: index out of range"); 01497 Vec<Num_T> Temp(*this); 01498 01499 set_size(datasize+in.length(), false); 01500 copy_vector(index, Temp.data, data); 01501 copy_vector(in.size(), in.data, &data[index]); 01502 copy_vector(Temp.datasize-index, Temp.data+index, data+index+in.size()); 01503 } 01504 01505 template<class Num_T> inline 01506 Vec<Num_T>& Vec<Num_T>::operator=(const Num_T t) 01507 { 01508 for (int i=0;i<datasize;i++) 01509 data[i] = t; 01510 return *this; 01511 } 01512 01513 template<class Num_T> inline 01514 Vec<Num_T>& Vec<Num_T>::operator=(const Vec<Num_T> &v) 01515 { 01516 if (this != &v) { 01517 set_size(v.datasize, false); 01518 copy_vector(datasize, v.data, data); 01519 } 01520 return *this; 01521 } 01522 01523 template<class Num_T> inline 01524 Vec<Num_T>& Vec<Num_T>::operator=(const Mat<Num_T> &m) 01525 { 01526 it_assert_debug( (m.cols() == 1 && datasize == m.rows()) || 01527 (m.rows() == 1 && datasize == m.cols()), "Vec::operator=(Mat<Num_T>): wrong size"); 01528 01529 if (m.cols() == 1) { 01530 set_size(m.rows(), false); 01531 copy_vector(m.rows(), m._data(), data); 01532 } else if (m.rows() == 1) { 01533 set_size(m.cols(), false); 01534 copy_vector(m.cols(), m._data(), m.rows(), data, 1); 01535 } else 01536 it_error("Vec::operator=(Mat<Num_T>): wrong size"); 01537 return *this; 01538 } 01539 01540 template<class Num_T> inline 01541 Vec<Num_T>& Vec<Num_T>::operator=(const char *values) 01542 { 01543 set(values); 01544 return *this; 01545 } 01546 01547 template<> 01548 bvec Vec<std::complex<double> >::operator==(std::complex<double>) const; 01549 01550 template<class Num_T> 01551 bvec Vec<Num_T>::operator==(const Num_T value) const 01552 { 01553 it_assert(datasize > 0, "Vec::operator==: vector must have size > 0"); 01554 Vec<Num_T> invector(*this); 01555 bvec temp(invector.length()); 01556 01557 for (int i=0;i<invector.length();i++) 01558 temp(i)=(invector(i)==value); 01559 01560 return temp; 01561 } 01562 01563 template<> 01564 bvec Vec<std::complex<double> >::operator!=(std::complex<double>) const; 01565 01566 template<class Num_T> 01567 bvec Vec<Num_T>::operator!=(const Num_T value) const 01568 { 01569 it_assert(datasize > 0, "Vec::operator!=: vector must have size > 0"); 01570 Vec<Num_T> invector(*this); 01571 bvec temp(invector.length()); 01572 01573 for (int i=0;i<invector.length();i++) 01574 temp(i)=(invector(i)!=value); 01575 01576 return temp; 01577 } 01578 01579 template<> 01580 bvec Vec<std::complex<double> >::operator<(std::complex<double>) const; 01581 01582 template<class Num_T> 01583 bvec Vec<Num_T>::operator<(const Num_T value) const 01584 { 01585 it_assert(datasize > 0, "Vec::operator<: vector must have size > 0"); 01586 Vec<Num_T> invector(*this); 01587 bvec temp(invector.length()); 01588 01589 for (int i=0;i<invector.length();i++) 01590 temp(i)=(invector(i)<value); 01591 01592 return temp; 01593 } 01594 01595 template<> 01596 bvec Vec<std::complex<double> >::operator<=(std::complex<double>) const; 01597 01598 template<class Num_T> 01599 bvec Vec<Num_T>::operator<=(const Num_T value) const 01600 { 01601 it_assert(datasize > 0, "Vec::operator<=: vector must have size > 0"); 01602 Vec<Num_T> invector(*this); 01603 bvec temp(invector.length()); 01604 01605 for (int i=0;i<invector.length();i++) 01606 temp(i)=(invector(i)<=value); 01607 01608 return temp; 01609 } 01610 01611 template<> 01612 bvec Vec<std::complex<double> >::operator>(std::complex<double>) const; 01613 01614 template<class Num_T> 01615 bvec Vec<Num_T>::operator>(const Num_T value) const 01616 { 01617 it_assert(datasize > 0, "Vec::operator>: vector must have size > 0"); 01618 Vec<Num_T> invector(*this); 01619 bvec temp(invector.length()); 01620 01621 for (int i=0;i<invector.length();i++) 01622 temp(i)=(invector(i)>value); 01623 01624 return temp; 01625 } 01626 01627 template<> 01628 bvec Vec<std::complex<double> >::operator>=(std::complex<double>) const; 01629 01630 template<class Num_T> 01631 bvec Vec<Num_T>::operator>=(const Num_T value) const 01632 { 01633 it_assert(datasize > 0, "Vec::operator>=: vector must have size > 0"); 01634 Vec<Num_T> invector(*this); 01635 bvec temp(invector.length()); 01636 01637 for (int i=0;i<invector.length();i++) 01638 temp(i)=(invector(i)>=value); 01639 01640 return temp; 01641 } 01642 01643 template<class Num_T> 01644 bool Vec<Num_T>::operator==(const Vec<Num_T> &invector) const 01645 { 01646 // OBS ! if wrong size, return false 01647 if (datasize!=invector.datasize) return false; 01648 for (int i=0;i<datasize;i++) { 01649 if (data[i]!=invector.data[i]) return false; 01650 } 01651 return true; 01652 } 01653 01654 template<class Num_T> 01655 bool Vec<Num_T>::operator!=(const Vec<Num_T> &invector) const 01656 { 01657 if (datasize!=invector.datasize) return true; 01658 for (int i=0;i<datasize;i++) { 01659 if (data[i]!=invector.data[i]) return true; 01660 } 01661 return false; 01662 } 01663 01665 template<class Num_T> 01666 std::ostream &operator<<(std::ostream &os, const Vec<Num_T> &v) 01667 { 01668 int i, sz=v.length(); 01669 01670 os << "[" ; 01671 for (i=0; i<sz; i++) { 01672 os << v(i) ; 01673 if (i < sz-1) 01674 os << " "; 01675 } 01676 os << "]" ; 01677 01678 return os; 01679 } 01680 01682 template<class Num_T> 01683 std::istream &operator>>(std::istream &is, Vec<Num_T> &v) 01684 { 01685 std::ostringstream buffer; 01686 bool started = false; 01687 bool finished = false; 01688 bool brackets = false; 01689 char c; 01690 01691 while (!finished) { 01692 if (is.eof()) { 01693 finished = true; 01694 } else { 01695 c = is.get(); 01696 01697 if (is.eof() || (c == '\n')) { 01698 if (brackets) { 01699 // Right bracket missing 01700 is.setstate(std::ios_base::failbit); 01701 finished = true; 01702 } else if (!((c == '\n') && !started)) { 01703 finished = true; 01704 } 01705 } else if ((c == ' ') || (c == '\t')) { 01706 if (started) { 01707 buffer << ' '; 01708 } 01709 } else if (c == '[') { 01710 if (started) { 01711 // Unexpected left bracket 01712 is.setstate(std::ios_base::failbit); 01713 finished = true; 01714 } else { 01715 started = true; 01716 brackets = true; 01717 } 01718 } else if (c == ']') { 01719 if (!started || !brackets) { 01720 // Unexpected right bracket 01721 is.setstate(std::ios_base::failbit); 01722 finished = true; 01723 } else { 01724 finished = true; 01725 } 01726 while (!is.eof() && (((c = is.peek()) == ' ') || (c == '\t'))) { 01727 is.get(); 01728 } 01729 if (!is.eof() && (c == '\n')) { 01730 is.get(); 01731 } 01732 } else { 01733 started = true; 01734 buffer << c; 01735 } 01736 } 01737 } 01738 01739 if (!started) { 01740 v.set_size(0, false); 01741 } else { 01742 v.set(buffer.str()); 01743 } 01744 01745 return is; 01746 } 01747 01749 01750 // ---------------------------------------------------------------------- 01751 // Instantiations 01752 // ---------------------------------------------------------------------- 01753 01754 #ifdef HAVE_EXTERN_TEMPLATE 01755 01756 extern template class Vec<double>; 01757 extern template class Vec<int>; 01758 extern template class Vec<short int>; 01759 extern template class Vec<std::complex<double> >; 01760 extern template class Vec<bin>; 01761 01762 // addition operator 01763 01764 extern template const vec operator+(const vec &v1, const vec &v2); 01765 extern template const cvec operator+(const cvec &v1, const cvec &v2); 01766 extern template const ivec operator+(const ivec &v1, const ivec &v2); 01767 extern template const svec operator+(const svec &v1, const svec &v2); 01768 extern template const bvec operator+(const bvec &v1, const bvec &v2); 01769 01770 extern template const vec operator+(const vec &v1, double t); 01771 extern template const cvec operator+(const cvec &v1, std::complex<double> t); 01772 extern template const ivec operator+(const ivec &v1, int t); 01773 extern template const svec operator+(const svec &v1, short t); 01774 extern template const bvec operator+(const bvec &v1, bin t); 01775 01776 extern template const vec operator+(double t, const vec &v1); 01777 extern template const cvec operator+(std::complex<double> t, const cvec &v1); 01778 extern template const ivec operator+(int t, const ivec &v1); 01779 extern template const svec operator+(short t, const svec &v1); 01780 extern template const bvec operator+(bin t, const bvec &v1); 01781 01782 // subraction operator 01783 01784 extern template const vec operator-(const vec &v1, const vec &v2); 01785 extern template const cvec operator-(const cvec &v1, const cvec &v2); 01786 extern template const ivec operator-(const ivec &v1, const ivec &v2); 01787 extern template const svec operator-(const svec &v1, const svec &v2); 01788 extern template const bvec operator-(const bvec &v1, const bvec &v2); 01789 01790 extern template const vec operator-(const vec &v, double t); 01791 extern template const cvec operator-(const cvec &v, std::complex<double> t); 01792 extern template const ivec operator-(const ivec &v, int t); 01793 extern template const svec operator-(const svec &v, short t); 01794 extern template const bvec operator-(const bvec &v, bin t); 01795 01796 extern template const vec operator-(double t, const vec &v); 01797 extern template const cvec operator-(std::complex<double> t, const cvec &v); 01798 extern template const ivec operator-(int t, const ivec &v); 01799 extern template const svec operator-(short t, const svec &v); 01800 extern template const bvec operator-(bin t, const bvec &v); 01801 01802 // unary minus 01803 01804 extern template const vec operator-(const vec &v); 01805 extern template const cvec operator-(const cvec &v); 01806 extern template const ivec operator-(const ivec &v); 01807 extern template const svec operator-(const svec &v); 01808 extern template const bvec operator-(const bvec &v); 01809 01810 // multiplication operator 01811 01812 #if !defined(HAVE_BLAS) 01813 extern template double dot(const vec &v1, const vec &v2); 01814 extern template std::complex<double> dot(const cvec &v1, const cvec &v2); 01815 #endif 01816 extern template int dot(const ivec &v1, const ivec &v2); 01817 extern template short dot(const svec &v1, const svec &v2); 01818 extern template bin dot(const bvec &v1, const bvec &v2); 01819 01820 #if !defined(HAVE_BLAS) 01821 extern template double operator*(const vec &v1, const vec &v2); 01822 extern template std::complex<double> operator*(const cvec &v1, 01823 const cvec &v2); 01824 #endif 01825 extern template int operator*(const ivec &v1, const ivec &v2); 01826 extern template short operator*(const svec &v1, const svec &v2); 01827 extern template bin operator*(const bvec &v1, const bvec &v2); 01828 01829 #if !defined(HAVE_BLAS) 01830 extern template const mat outer_product(const vec &v1, const vec &v2, 01831 bool hermitian); 01832 #endif 01833 extern template const imat outer_product(const ivec &v1, const ivec &v2, 01834 bool hermitian); 01835 extern template const smat outer_product(const svec &v1, const svec &v2, 01836 bool hermitian); 01837 extern template const bmat outer_product(const bvec &v1, const bvec &v2, 01838 bool hermitian); 01839 01840 extern template const vec operator*(const vec &v, double t); 01841 extern template const cvec operator*(const cvec &v, std::complex<double> t); 01842 extern template const ivec operator*(const ivec &v, int t); 01843 extern template const svec operator*(const svec &v, short t); 01844 extern template const bvec operator*(const bvec &v, bin t); 01845 01846 extern template const vec operator*(double t, const vec &v); 01847 extern template const cvec operator*(std::complex<double> t, const cvec &v); 01848 extern template const ivec operator*(int t, const ivec &v); 01849 extern template const svec operator*(short t, const svec &v); 01850 extern template const bvec operator*(bin t, const bvec &v); 01851 01852 // elementwise multiplication 01853 01854 extern template const vec elem_mult(const vec &a, const vec &b); 01855 extern template const cvec elem_mult(const cvec &a, const cvec &b); 01856 extern template const ivec elem_mult(const ivec &a, const ivec &b); 01857 extern template const svec elem_mult(const svec &a, const svec &b); 01858 extern template const bvec elem_mult(const bvec &a, const bvec &b); 01859 01860 extern template void elem_mult_out(const vec &a, const vec &b, vec &out); 01861 extern template void elem_mult_out(const cvec &a, const cvec &b, cvec &out); 01862 extern template void elem_mult_out(const ivec &a, const ivec &b, ivec &out); 01863 extern template void elem_mult_out(const svec &a, const svec &b, svec &out); 01864 extern template void elem_mult_out(const bvec &a, const bvec &b, bvec &out); 01865 01866 extern template const vec elem_mult(const vec &a, const vec &b, 01867 const vec &c); 01868 extern template const cvec elem_mult(const cvec &a, const cvec &b, 01869 const cvec &c); 01870 extern template const ivec elem_mult(const ivec &a, const ivec &b, 01871 const ivec &c); 01872 extern template const svec elem_mult(const svec &a, const svec &b, 01873 const svec &c); 01874 extern template const bvec elem_mult(const bvec &a, const bvec &b, 01875 const bvec &c); 01876 01877 extern template void elem_mult_out(const vec &a, const vec &b, const vec &c, 01878 vec &out); 01879 extern template void elem_mult_out(const cvec &a, const cvec &b, 01880 const cvec &c, cvec &out); 01881 extern template void elem_mult_out(const ivec &a, const ivec &b, 01882 const ivec &c, ivec &out); 01883 extern template void elem_mult_out(const svec &a, const svec &b, 01884 const svec &c, svec &out); 01885 extern template void elem_mult_out(const bvec &a, const bvec &b, 01886 const bvec &c, bvec &out); 01887 01888 extern template const vec elem_mult(const vec &a, const vec &b, const vec &c, 01889 const vec &d); 01890 extern template const cvec elem_mult(const cvec &a, const cvec &b, 01891 const cvec &c, const cvec &d); 01892 extern template const ivec elem_mult(const ivec &a, const ivec &b, 01893 const ivec &c, const ivec &d); 01894 extern template const svec elem_mult(const svec &a, const svec &b, 01895 const svec &c, const svec &d); 01896 extern template const bvec elem_mult(const bvec &a, const bvec &b, 01897 const bvec &c, const bvec &d); 01898 01899 extern template void elem_mult_out(const vec &a, const vec &b, const vec &c, 01900 const vec &d, vec &out); 01901 extern template void elem_mult_out(const cvec &a, const cvec &b, 01902 const cvec &c, const cvec &d, cvec &out); 01903 extern template void elem_mult_out(const ivec &a, const ivec &b, 01904 const ivec &c, const ivec &d, ivec &out); 01905 extern template void elem_mult_out(const svec &a, const svec &b, 01906 const svec &c, const svec &d, svec &out); 01907 extern template void elem_mult_out(const bvec &a, const bvec &b, 01908 const bvec &c, const bvec &d, bvec &out); 01909 01910 // in-place elementwise multiplication 01911 01912 extern template void elem_mult_inplace(const vec &a, vec &b); 01913 extern template void elem_mult_inplace(const cvec &a, cvec &b); 01914 extern template void elem_mult_inplace(const ivec &a, ivec &b); 01915 extern template void elem_mult_inplace(const svec &a, svec &b); 01916 extern template void elem_mult_inplace(const bvec &a, bvec &b); 01917 01918 // elementwise multiplication followed by summation 01919 01920 extern template double elem_mult_sum(const vec &a, const vec &b); 01921 extern template std::complex<double> elem_mult_sum(const cvec &a, 01922 const cvec &b); 01923 extern template int elem_mult_sum(const ivec &a, const ivec &b); 01924 extern template short elem_mult_sum(const svec &a, const svec &b); 01925 extern template bin elem_mult_sum(const bvec &a, const bvec &b); 01926 01927 // division operator 01928 01929 extern template const vec operator/(const vec &v, double t); 01930 extern template const cvec operator/(const cvec &v, std::complex<double> t); 01931 extern template const ivec operator/(const ivec &v, int t); 01932 extern template const svec operator/(const svec &v, short t); 01933 extern template const bvec operator/(const bvec &v, bin t); 01934 01935 extern template const vec operator/(double t, const vec &v); 01936 extern template const cvec operator/(std::complex<double> t, const cvec &v); 01937 extern template const ivec operator/(int t, const ivec &v); 01938 extern template const svec operator/(short t, const svec &v); 01939 extern template const bvec operator/(bin t, const bvec &v); 01940 01941 // elementwise division operator 01942 01943 extern template const vec elem_div(const vec &a, const vec &b); 01944 extern template const cvec elem_div(const cvec &a, const cvec &b); 01945 extern template const ivec elem_div(const ivec &a, const ivec &b); 01946 extern template const svec elem_div(const svec &a, const svec &b); 01947 extern template const bvec elem_div(const bvec &a, const bvec &b); 01948 01949 extern template const vec elem_div(double t, const vec &v); 01950 extern template const cvec elem_div(std::complex<double> t, const cvec &v); 01951 extern template const ivec elem_div(int t, const ivec &v); 01952 extern template const svec elem_div(short t, const svec &v); 01953 extern template const bvec elem_div(bin t, const bvec &v); 01954 01955 extern template void elem_div_out(const vec &a, const vec &b, vec &out); 01956 extern template void elem_div_out(const cvec &a, const cvec &b, cvec &out); 01957 extern template void elem_div_out(const ivec &a, const ivec &b, ivec &out); 01958 extern template void elem_div_out(const svec &a, const svec &b, svec &out); 01959 extern template void elem_div_out(const bvec &a, const bvec &b, bvec &out); 01960 01961 // elementwise division followed by summation 01962 01963 extern template double elem_div_sum(const vec &a, const vec &b); 01964 extern template std::complex<double> elem_div_sum(const cvec &a, 01965 const cvec &b); 01966 extern template int elem_div_sum(const ivec &a, const ivec &b); 01967 extern template short elem_div_sum(const svec &a, const svec &b); 01968 extern template bin elem_div_sum(const bvec &a, const bvec &b); 01969 01970 // concat operator 01971 01972 extern template const vec concat(const vec &v, double a); 01973 extern template const cvec concat(const cvec &v, std::complex<double> a); 01974 extern template const ivec concat(const ivec &v, int a); 01975 extern template const svec concat(const svec &v, short a); 01976 extern template const bvec concat(const bvec &v, bin a); 01977 01978 extern template const vec concat(double a, const vec &v); 01979 extern template const cvec concat(std::complex<double> a, const cvec &v); 01980 extern template const ivec concat(int a, const ivec &v); 01981 extern template const svec concat(short a, const svec &v); 01982 extern template const bvec concat(bin a, const bvec &v); 01983 01984 extern template const vec concat(const vec &v1, const vec &v2); 01985 extern template const cvec concat(const cvec &v1, const cvec &v2); 01986 extern template const ivec concat(const ivec &v1, const ivec &v2); 01987 extern template const svec concat(const svec &v1, const svec &v2); 01988 extern template const bvec concat(const bvec &v1, const bvec &v2); 01989 01990 extern template const vec concat(const vec &v1, const vec &v2, 01991 const vec &v3); 01992 extern template const cvec concat(const cvec &v1, const cvec &v2, 01993 const cvec &v3); 01994 extern template const ivec concat(const ivec &v1, const ivec &v2, 01995 const ivec &v3); 01996 extern template const svec concat(const svec &v1, const svec &v2, 01997 const svec &v3); 01998 extern template const bvec concat(const bvec &v1, const bvec &v2, 01999 const bvec &v3); 02000 02001 extern template const vec concat(const vec &v1, const vec &v2, const vec &v3, 02002 const vec &v4); 02003 extern template const cvec concat(const cvec &v1, const cvec &v2, 02004 const cvec &v3, const cvec &v4); 02005 extern template const ivec concat(const ivec &v1, const ivec &v2, 02006 const ivec &v3, const ivec &v4); 02007 extern template const svec concat(const svec &v1, const svec &v2, 02008 const svec &v3, const svec &v4); 02009 extern template const bvec concat(const bvec &v1, const bvec &v2, 02010 const bvec &v3, const bvec &v4); 02011 02012 extern template const vec concat(const vec &v1, const vec &v2, const vec &v3, 02013 const vec &v4, const vec &v5); 02014 extern template const cvec concat(const cvec &v1, const cvec &v2, 02015 const cvec &v3, const cvec &v4, 02016 const cvec &v5); 02017 extern template const ivec concat(const ivec &v1, const ivec &v2, 02018 const ivec &v3, const ivec &v4, 02019 const ivec &v5); 02020 extern template const svec concat(const svec &v1, const svec &v2, 02021 const svec &v3, const svec &v4, 02022 const svec &v5); 02023 extern template const bvec concat(const bvec &v1, const bvec &v2, 02024 const bvec &v3, const bvec &v4, 02025 const bvec &v5); 02026 02027 // I/O streams 02028 02029 extern template std::ostream &operator<<(std::ostream& os, const vec &vect); 02030 extern template std::ostream &operator<<(std::ostream& os, const cvec &vect); 02031 extern template std::ostream &operator<<(std::ostream& os, const svec &vect); 02032 extern template std::ostream &operator<<(std::ostream& os, const ivec &vect); 02033 extern template std::ostream &operator<<(std::ostream& os, const bvec &vect); 02034 extern template std::istream &operator>>(std::istream& is, vec &vect); 02035 extern template std::istream &operator>>(std::istream& is, cvec &vect); 02036 extern template std::istream &operator>>(std::istream& is, svec &vect); 02037 extern template std::istream &operator>>(std::istream& is, ivec &vect); 02038 extern template std::istream &operator>>(std::istream& is, bvec &vect); 02039 02040 #endif // HAVE_EXTERN_TEMPLATE 02041 02043 02044 } // namespace itpp 02045 02046 #endif // #ifndef VEC_H
Generated on Mon Jan 7 22:28:57 2008 for IT++ by Doxygen 1.5.4