00001 00033 #ifndef SVEC_H 00034 #define SVEC_H 00035 00036 #include <itpp/base/vec.h> 00037 #include <itpp/base/stat.h> 00038 #include <cstdlib> 00039 00040 00041 namespace itpp { 00042 00043 // Declaration of class Vec 00044 template <class T> class Vec; 00045 // Declaration of class Sparse_Vec 00046 template <class T> class Sparse_Vec; 00047 00048 // ----------------------- Sparse_Vec Friends ------------------------------- 00049 00051 template <class T> 00052 Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00053 00055 template <class T> 00056 T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00057 00059 template <class T> 00060 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00061 00063 template <class T> 00064 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00065 00067 template <class T> 00068 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00069 00071 template <class T> 00072 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00073 00075 template <class T> 00076 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00077 00079 template <class T> 00080 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00081 00083 template <class T> 00084 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00085 00094 template <class T> 00095 class Sparse_Vec { 00096 public: 00097 00099 Sparse_Vec(); 00100 00107 Sparse_Vec(int sz, int data_init=200); 00108 00114 Sparse_Vec(const Sparse_Vec<T> &v); 00115 00121 Sparse_Vec(const Vec<T> &v); 00122 00128 Sparse_Vec(const Vec<T> &v, T epsilon); 00129 00131 ~Sparse_Vec(); 00132 00139 void set_size(int sz, int data_init=-1); 00140 00142 int size() const { return v_size; } 00143 00145 inline int nnz() 00146 { 00147 if (check_small_elems_flag) { 00148 remove_small_elements(); 00149 } 00150 return used_size; 00151 } 00152 00154 double density(); 00155 00157 void set_small_element(const T& epsilon); 00158 00164 void remove_small_elements(); 00165 00171 void resize_data(int new_size); 00172 00174 void compact(); 00175 00177 void full(Vec<T> &v) const; 00178 00180 Vec<T> full() const; 00181 00183 T operator()(int i) const; 00184 00186 void set(int i, T v); 00187 00189 void set(const ivec &index_vec, const Vec<T> &v); 00190 00192 void set_new(int i, T v); 00193 00195 void set_new(const ivec &index_vec, const Vec<T> &v); 00196 00198 void add_elem(const int i, const T v); 00199 00201 void add(const ivec &index_vec, const Vec<T> &v); 00202 00204 void zeros(); 00205 00207 void zero_elem(const int i); 00208 00210 void clear(); 00211 00213 void clear_elem(const int i); 00214 00218 inline void get_nz_data(int p, T& data_out) 00219 { 00220 if (check_small_elems_flag) { 00221 remove_small_elements(); 00222 } 00223 data_out = data[p]; 00224 } 00225 00227 inline T get_nz_data(int p) 00228 { 00229 if (check_small_elems_flag) { 00230 remove_small_elements(); 00231 } 00232 return data[p]; 00233 } 00234 00236 inline int get_nz_index(int p) 00237 { 00238 if (check_small_elems_flag) { 00239 remove_small_elements(); 00240 } 00241 return index[p]; 00242 } 00243 00245 inline void get_nz(int p, int &idx, T &dat) 00246 { 00247 if (check_small_elems_flag) { 00248 remove_small_elements(); 00249 } 00250 idx=index[p]; 00251 dat=data[p]; 00252 } 00253 00255 Sparse_Vec<T> get_subvector(int i1, int i2) const; 00256 00258 T sqr() const; 00259 00261 void operator=(const Sparse_Vec<T> &v); 00263 void operator=(const Vec<T> &v); 00264 00266 Sparse_Vec<T> operator-() const; 00267 00269 bool operator==(const Sparse_Vec<T> &v); 00270 00272 void operator+=(const Sparse_Vec<T> &v); 00273 00275 void operator+=(const Vec<T> &v); 00276 00278 void operator-=(const Sparse_Vec<T> &v); 00279 00281 void operator-=(const Vec<T> &v); 00282 00284 void operator*=(const T &v); 00285 00287 void operator/=(const T &v); 00288 00290 friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00292 friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00294 friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00296 friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00297 00299 friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00300 00302 friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00303 00305 friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00306 00308 friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00309 00311 friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00312 00313 private: 00314 void init(); 00315 void alloc(); 00316 void free(); 00317 00318 int v_size, used_size, data_size; 00319 T *data; 00320 int *index; 00321 T eps; 00322 bool check_small_elems_flag; 00323 }; 00324 00329 typedef Sparse_Vec<int> sparse_ivec; 00330 00335 typedef Sparse_Vec<double> sparse_vec; 00336 00341 typedef Sparse_Vec<std::complex<double> > sparse_cvec; 00342 00343 // ----------------------- Implementation starts here -------------------------------- 00344 00345 template <class T> 00346 void Sparse_Vec<T>::init() 00347 { 00348 v_size = 0; 00349 used_size = 0; 00350 data_size = 0; 00351 data = 0; 00352 index = 0; 00353 eps = 0; 00354 check_small_elems_flag = true; 00355 } 00356 00357 template <class T> 00358 void Sparse_Vec<T>::alloc() 00359 { 00360 if (data_size != 0) { 00361 data = new T[data_size]; 00362 index = new int[data_size]; 00363 } 00364 } 00365 00366 template <class T> 00367 void Sparse_Vec<T>::free() 00368 { 00369 delete [] data; 00370 data = 0; 00371 delete [] index; 00372 index = 0; 00373 } 00374 00375 template <class T> 00376 Sparse_Vec<T>::Sparse_Vec() 00377 { 00378 init(); 00379 } 00380 00381 template <class T> 00382 Sparse_Vec<T>::Sparse_Vec(int sz, int data_init) 00383 { 00384 init(); 00385 v_size = sz; 00386 used_size = 0; 00387 data_size = data_init; 00388 alloc(); 00389 } 00390 00391 template <class T> 00392 Sparse_Vec<T>::Sparse_Vec(const Sparse_Vec<T> &v) 00393 { 00394 init(); 00395 v_size = v.v_size; 00396 used_size = v.used_size; 00397 data_size = v.data_size; 00398 eps = v.eps; 00399 check_small_elems_flag = v.check_small_elems_flag; 00400 alloc(); 00401 00402 for (int i=0; i<used_size; i++) { 00403 data[i] = v.data[i]; 00404 index[i] = v.index[i]; 00405 } 00406 } 00407 00408 template <class T> 00409 Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v) 00410 { 00411 init(); 00412 v_size = v.size(); 00413 used_size = 0; 00414 data_size = std::min(v.size(), 10000); 00415 alloc(); 00416 00417 for (int i=0; i<v_size; i++) { 00418 if (v(i) != T(0)) { 00419 if (used_size == data_size) 00420 resize_data(data_size*2); 00421 data[used_size] = v(i); 00422 index[used_size] = i; 00423 used_size++; 00424 } 00425 } 00426 compact(); 00427 } 00428 00429 template <class T> 00430 Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon) 00431 { 00432 init(); 00433 v_size = v.size(); 00434 used_size = 0; 00435 data_size = std::min(v.size(), 10000); 00436 eps = epsilon; 00437 alloc(); 00438 00439 double e = std::abs(epsilon); 00440 for (int i=0; i<v_size; i++) { 00441 if (std::abs(v(i)) > e) { 00442 if (used_size == data_size) 00443 resize_data(data_size*2); 00444 data[used_size] = v(i); 00445 index[used_size] = i; 00446 used_size++; 00447 } 00448 } 00449 compact(); 00450 } 00451 00452 template <class T> 00453 Sparse_Vec<T>::~Sparse_Vec() 00454 { 00455 free(); 00456 } 00457 00458 template <class T> 00459 void Sparse_Vec<T>::set_size(int new_size, int data_init) 00460 { 00461 v_size = new_size; 00462 used_size = 0; 00463 if (data_init != -1){ 00464 free(); 00465 data_size = data_init; 00466 alloc(); 00467 } 00468 } 00469 00470 template <class T> 00471 double Sparse_Vec<T>::density() 00472 { 00473 if (check_small_elems_flag) { 00474 remove_small_elements(); 00475 } 00476 //return static_cast<double>(used_size) / v_size; 00477 return double(used_size) / v_size; 00478 } 00479 00480 template <class T> 00481 void Sparse_Vec<T>::set_small_element(const T& epsilon) 00482 { 00483 eps=epsilon; 00484 remove_small_elements(); 00485 } 00486 00487 template <class T> 00488 void Sparse_Vec<T>::remove_small_elements() 00489 { 00490 int i; 00491 int nrof_removed_elements = 0; 00492 double e; 00493 00494 //Remove small elements 00495 e = std::abs(eps); 00496 for (i=0;i<used_size;i++) { 00497 if (std::abs(data[i]) <= e) { 00498 nrof_removed_elements++; 00499 } 00500 else if (nrof_removed_elements > 0) { 00501 data[i-nrof_removed_elements] = data[i]; 00502 index[i-nrof_removed_elements] = index[i]; 00503 } 00504 } 00505 00506 //Set new size after small elements have been removed 00507 used_size -= nrof_removed_elements; 00508 00509 //Set the flag to indicate that all small elements have been removed 00510 check_small_elems_flag = false; 00511 } 00512 00513 00514 template <class T> 00515 void Sparse_Vec<T>::resize_data(int new_size) 00516 { 00517 it_assert(new_size>=used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small"); 00518 00519 if (new_size != data_size) { 00520 if (new_size == 0) 00521 free(); 00522 else { 00523 T *tmp_data = data; 00524 int *tmp_pos = index; 00525 data_size = new_size; 00526 alloc(); 00527 for (int p=0; p<used_size; p++) { 00528 data[p] = tmp_data[p]; 00529 index[p] = tmp_pos[p]; 00530 } 00531 delete [] tmp_data; 00532 delete [] tmp_pos; 00533 } 00534 } 00535 } 00536 00537 template <class T> 00538 void Sparse_Vec<T>::compact() 00539 { 00540 if (check_small_elems_flag) { 00541 remove_small_elements(); 00542 } 00543 resize_data(used_size); 00544 } 00545 00546 template <class T> 00547 void Sparse_Vec<T>::full(Vec<T> &v) const 00548 { 00549 v.set_size(v_size); 00550 00551 v = T(0); 00552 for (int p=0; p<used_size; p++) 00553 v(index[p]) = data[p]; 00554 } 00555 00556 template <class T> 00557 Vec<T> Sparse_Vec<T>::full() const 00558 { 00559 Vec<T> r(v_size); 00560 full(r); 00561 return r; 00562 } 00563 00564 // This is slow. Implement a better search 00565 template <class T> 00566 T Sparse_Vec<T>::operator()(int i) const 00567 { 00568 it_assert0(i >= 0 && i < v_size, "The index of the element is out of range"); 00569 00570 bool found = false; 00571 int p; 00572 for (p=0; p<used_size; p++) { 00573 if (index[p] == i) { 00574 found = true; 00575 break; 00576 } 00577 } 00578 return found ? data[p] : T(0); 00579 } 00580 00581 template <class T> 00582 void Sparse_Vec<T>::set(int i, T v) 00583 { 00584 it_assert0(i >= 0 && i < v_size, "The index of the element is out of range"); 00585 00586 bool found = false; 00587 bool larger_than_eps; 00588 int p; 00589 00590 for (p=0; p<used_size; p++) { 00591 if (index[p] == i) { 00592 found = true; 00593 break; 00594 } 00595 } 00596 00597 larger_than_eps = (std::abs(v) > std::abs(eps)); 00598 00599 if (found && larger_than_eps) 00600 data[p] = v; 00601 else if (larger_than_eps) { 00602 if (used_size == data_size) 00603 resize_data(data_size*2+100); 00604 data[used_size] = v; 00605 index[used_size] = i; 00606 used_size++; 00607 } 00608 00609 //Check if the stored element is smaller than eps. In that case it should be removed. 00610 if (std::abs(v) <= std::abs(eps)) { 00611 remove_small_elements(); 00612 } 00613 00614 } 00615 00616 template <class T> 00617 void Sparse_Vec<T>::set_new(int i, T v) 00618 { 00619 it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00620 00621 //Check that the new element is larger than eps! 00622 if (std::abs(v) > std::abs(eps)) { 00623 if (used_size == data_size) 00624 resize_data(data_size*2+100); 00625 data[used_size] = v; 00626 index[used_size] = i; 00627 used_size++; 00628 } 00629 } 00630 00631 template <class T> 00632 void Sparse_Vec<T>::add_elem(const int i, const T v) 00633 { 00634 bool found = false; 00635 int p; 00636 00637 it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00638 00639 for (p=0; p<used_size; p++) { 00640 if (index[p] == i) { 00641 found = true; 00642 break; 00643 } 00644 } 00645 if (found) 00646 data[p] += v; 00647 else { 00648 if (used_size == data_size) 00649 resize_data(data_size*2+100); 00650 data[used_size] = v; 00651 index[used_size] = i; 00652 used_size++; 00653 } 00654 00655 check_small_elems_flag = true; 00656 00657 } 00658 00659 template <class T> 00660 void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v) 00661 { 00662 bool found = false; 00663 int i,p,q; 00664 int nrof_nz=v.size(); 00665 00666 it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector"); 00667 00668 //Elements are added if they have identical indices 00669 for (q=0; q<nrof_nz; q++){ 00670 i=index_vec(q); 00671 for (p=0; p<used_size; p++) { 00672 if (index[p] == i) { 00673 found = true; 00674 break; 00675 } 00676 } 00677 if (found) 00678 data[p] += v(q); 00679 else { 00680 if (used_size == data_size) 00681 resize_data(data_size*2+100); 00682 data[used_size] = v(q); 00683 index[used_size] = i; 00684 used_size++; 00685 } 00686 found = false; 00687 } 00688 00689 check_small_elems_flag = true; 00690 00691 } 00692 00693 template <class T> 00694 void Sparse_Vec<T>::zeros() 00695 { 00696 used_size = 0; 00697 check_small_elems_flag = false; 00698 } 00699 00700 template <class T> 00701 void Sparse_Vec<T>::zero_elem(const int i) 00702 { 00703 bool found = false; 00704 int p; 00705 00706 it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00707 00708 for (p=0; p<used_size; p++) { 00709 if (index[p] == i) { 00710 found = true; 00711 break; 00712 } 00713 } 00714 if (found) { 00715 data[p] = data[used_size-1]; 00716 index[p] = index[used_size-1]; 00717 used_size--; 00718 } 00719 } 00720 00721 template <class T> 00722 void Sparse_Vec<T>::clear() 00723 { 00724 used_size = 0; 00725 check_small_elems_flag = false; 00726 } 00727 00728 template <class T> 00729 void Sparse_Vec<T>::clear_elem(const int i) 00730 { 00731 bool found = false; 00732 int p; 00733 00734 it_assert0(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00735 00736 for (p=0; p<used_size; p++) { 00737 if (index[p] == i) { 00738 found = true; 00739 break; 00740 } 00741 } 00742 if (found) { 00743 data[p] = data[used_size-1]; 00744 index[p] = index[used_size-1]; 00745 used_size--; 00746 } 00747 } 00748 00749 template <class T> 00750 void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v) 00751 { 00752 it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector"); 00753 00754 //Clear all old non-zero elements 00755 clear(); 00756 00757 //Add the new non-zero elements 00758 add(index_vec,v); 00759 } 00760 00761 template <class T> 00762 void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v) 00763 { 00764 int q; 00765 int nrof_nz=v.size(); 00766 00767 it_assert0(v_size>max(index_vec),"The indices exceeds the size of the sparse vector"); 00768 00769 //Clear all old non-zero elements 00770 clear(); 00771 00772 for (q=0; q<nrof_nz; q++){ 00773 if (std::abs(v[q]) > std::abs(eps)) { 00774 if (used_size == data_size) 00775 resize_data(data_size*2+100); 00776 data[used_size] = v(q); 00777 index[used_size] = index_vec(q); 00778 used_size++; 00779 } 00780 } 00781 } 00782 00783 template <class T> 00784 Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const 00785 { 00786 it_assert0(v_size > i1 && v_size > i2 && i1<=i2 && i1>=0, "The index of the element exceeds the size of the sparse vector"); 00787 00788 Sparse_Vec<T> r(i2-i1+1); 00789 00790 for (int p=0; p<used_size; p++) { 00791 if (index[p] >= i1 && index[p] <= i2) { 00792 if (r.used_size == r.data_size) 00793 r.resize_data(r.data_size*2+100); 00794 r.data[r.used_size] = data[p]; 00795 r.index[r.used_size] = index[p]-i1; 00796 r.used_size++; 00797 } 00798 } 00799 r.eps = eps; 00800 r.check_small_elems_flag = check_small_elems_flag; 00801 r.compact(); 00802 00803 return r; 00804 } 00805 00806 template <class T> 00807 T Sparse_Vec<T>::sqr() const 00808 { 00809 T sum(0); 00810 for (int p=0; p<used_size; p++) 00811 sum += data[p] * data[p]; 00812 00813 return sum; 00814 } 00815 00816 template <class T> 00817 void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v) 00818 { 00819 free(); 00820 v_size = v.v_size; 00821 used_size = v.used_size; 00822 data_size = v.data_size; 00823 eps = v.eps; 00824 check_small_elems_flag = v.check_small_elems_flag; 00825 alloc(); 00826 00827 for (int i=0; i<used_size; i++) { 00828 data[i] = v.data[i]; 00829 index[i] = v.index[i]; 00830 } 00831 } 00832 00833 template <class T> 00834 void Sparse_Vec<T>::operator=(const Vec<T> &v) 00835 { 00836 free(); 00837 v_size = v.size(); 00838 used_size = 0; 00839 data_size = std::min(v.size(), 10000); 00840 eps = T(0); 00841 check_small_elems_flag = false; 00842 alloc(); 00843 00844 for (int i=0; i<v_size; i++) { 00845 if (v(i) != T(0)) { 00846 if (used_size == data_size) 00847 resize_data(data_size*2); 00848 data[used_size] = v(i); 00849 index[used_size] = i; 00850 used_size++; 00851 } 00852 } 00853 compact(); 00854 } 00855 00856 template <class T> 00857 Sparse_Vec<T> Sparse_Vec<T>::operator-() const 00858 { 00859 Sparse_Vec r(v_size, used_size); 00860 00861 for (int p=0; p<used_size; p++) { 00862 r.data[p] = -data[p]; 00863 r.index[p] = index[p]; 00864 } 00865 r.used_size = used_size; 00866 00867 return r; 00868 } 00869 00870 template <class T> 00871 bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v) 00872 { 00873 int p,q; 00874 bool found=false; 00875 00876 //Remove small elements before comparing the two sparse_vectors 00877 if (check_small_elems_flag) 00878 remove_small_elements(); 00879 00880 if (v_size!=v.v_size) { 00881 //Return false if vector sizes are unequal 00882 return false; 00883 } else { 00884 for(p=0;p<used_size;p++) { 00885 for(q=0;q<v.used_size;q++) { 00886 if (index[p] == v.index[q]) { 00887 found = true; 00888 break; 00889 } 00890 } 00891 if (found==false) 00892 //Return false if non-zero element not found, or if elements are unequal 00893 return false; 00894 else if (data[p]!=v.data[q]) 00895 //Return false if non-zero element not found, or if elements are unequal 00896 return false; 00897 else 00898 //Check next non-zero element 00899 found = false; 00900 } 00901 } 00902 00903 /*Special handling if sizes do not match. 00904 Required since v may need to do remove_small_elements() for true comparison*/ 00905 if (used_size!=v.used_size) { 00906 if (used_size > v.used_size) { 00907 //Return false if number of non-zero elements is less in v 00908 return false; 00909 } 00910 else { 00911 //Ensure that the remaining non-zero elements in v are smaller than v.eps 00912 int nrof_small_elems = 0; 00913 for(q=0;q<v.used_size;q++) { 00914 if (std::abs(v.data[q]) <= std::abs(v.eps)) 00915 nrof_small_elems++; 00916 } 00917 if (v.used_size-nrof_small_elems != used_size) 00918 //Return false if the number of "true" non-zero elements are unequal 00919 return false; 00920 } 00921 } 00922 00923 //All elements checks => return true 00924 return true; 00925 } 00926 00927 template <class T> 00928 void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v) 00929 { 00930 int i,p; 00931 T tmp_data; 00932 int nrof_nz_v=v.used_size; 00933 00934 it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00935 00936 for (p=0; p<nrof_nz_v; p++) { 00937 i = v.index[p]; 00938 tmp_data = v.data[p]; 00939 //get_nz(p,i,tmp_data); 00940 add_elem(i,tmp_data); 00941 } 00942 00943 check_small_elems_flag = true; 00944 } 00945 00946 template <class T> 00947 void Sparse_Vec<T>::operator+=(const Vec<T> &v) 00948 { 00949 int i; 00950 00951 it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00952 00953 for (i=0; i<v.size(); i++) 00954 if (v(i)!=T(0)) 00955 add_elem(i,v(i)); 00956 00957 check_small_elems_flag = true; 00958 } 00959 00960 00961 template <class T> 00962 void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v) 00963 { 00964 int i,p; 00965 T tmp_data; 00966 int nrof_nz_v=v.used_size; 00967 00968 it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 00969 00970 for (p=0; p<nrof_nz_v; p++) { 00971 i = v.index[p]; 00972 tmp_data = v.data[p]; 00973 //v.get_nz(p,i,tmp_data); 00974 add_elem(i,-tmp_data); 00975 } 00976 00977 check_small_elems_flag = true; 00978 } 00979 00980 template <class T> 00981 void Sparse_Vec<T>::operator-=(const Vec<T> &v) 00982 { 00983 int i; 00984 00985 it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 00986 00987 for (i=0; i<v.size(); i++) 00988 if (v(i)!=T(0)) 00989 add_elem(i,-v(i)); 00990 00991 check_small_elems_flag = true; 00992 } 00993 00994 template <class T> 00995 void Sparse_Vec<T>::operator*=(const T &v) 00996 { 00997 int p; 00998 00999 for (p=0; p<used_size; p++) { 01000 data[p]*=v;} 01001 01002 check_small_elems_flag = true; 01003 } 01004 01005 template <class T> 01006 void Sparse_Vec<T>::operator/=(const T &v) 01007 { 01008 int p; 01009 for (p=0; p<used_size; p++) { 01010 data[p]/=v;} 01011 01012 if (std::abs(eps) > 0) { 01013 check_small_elems_flag = true; 01014 } 01015 } 01016 01017 template <class T> 01018 T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01019 { 01020 it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> * Sparse_Vec<T>"); 01021 01022 T sum(0); 01023 Vec<T> v1f(v1.v_size); 01024 v1.full(v1f); 01025 for (int p=0; p<v2.used_size; p++) { 01026 if (v1f[v2.index[p]] != T(0)) 01027 sum += v1f[v2.index[p]] * v2.data[p]; 01028 } 01029 01030 return sum; 01031 } 01032 01033 template <class T> 01034 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01035 { 01036 it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted"); 01037 01038 T sum(0); 01039 for (int p1=0; p1<v1.used_size; p1++) 01040 sum += v1.data[p1] * v2[v1.index[p1]]; 01041 01042 return sum; 01043 } 01044 01045 template <class T> 01046 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01047 { 01048 it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted"); 01049 01050 T sum(0); 01051 for (int p2=0; p2<v2.used_size; p2++) 01052 sum += v1[v2.index[p2]] * v2.data[p2]; 01053 01054 return sum; 01055 } 01056 01057 template <class T> 01058 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01059 { 01060 it_assert0(v1.v_size == v2.v_size, "elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)"); 01061 01062 Sparse_Vec<T> r(v1.v_size); 01063 ivec pos(v1.v_size); 01064 pos = -1; 01065 for (int p1=0; p1<v1.used_size; p1++) 01066 pos[v1.index[p1]] = p1; 01067 for (int p2=0; p2<v2.used_size; p2++) { 01068 if (pos[v2.index[p2]] != -1) { 01069 if (r.used_size == r.data_size) 01070 r.resize_data(r.used_size*2+100); 01071 r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2]; 01072 r.index[r.used_size] = v2.index[p2]; 01073 r.used_size++; 01074 } 01075 } 01076 r.compact(); 01077 01078 return r; 01079 } 01080 01081 template <class T> 01082 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01083 { 01084 it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)"); 01085 01086 Vec<T> r(v1.v_size); 01087 r = T(0); 01088 for (int p1=0; p1<v1.used_size; p1++) 01089 r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]]; 01090 01091 return r; 01092 } 01093 01094 template <class T> 01095 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01096 { 01097 it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)"); 01098 01099 Sparse_Vec<T> r(v1.v_size); 01100 for (int p1=0; p1<v1.used_size; p1++) { 01101 if (v2[v1.index[p1]] != T(0)) { 01102 if (r.used_size == r.data_size) 01103 r.resize_data(r.used_size*2+100); 01104 r.data[r.used_size] = v1.data[p1] * v2[v1.index[p1]]; 01105 r.index[r.used_size] = v1.index[p1]; 01106 r.used_size++; 01107 } 01108 } 01109 r.compact(); 01110 01111 return r; 01112 } 01113 01114 template <class T> 01115 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01116 { 01117 it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)"); 01118 01119 Vec<T> r(v2.v_size); 01120 r = T(0); 01121 for (int p2=0; p2<v2.used_size; p2++) 01122 r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2]; 01123 01124 return r; 01125 } 01126 01127 template <class T> 01128 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01129 { 01130 it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)"); 01131 01132 Sparse_Vec<T> r(v2.v_size); 01133 for (int p2=0; p2<v2.used_size; p2++) { 01134 if (v1[v2.index[p2]] != T(0)) { 01135 if (r.used_size == r.data_size) 01136 r.resize_data(r.used_size*2+100); 01137 r.data[r.used_size] = v1[v2.index[p2]] * v2.data[p2]; 01138 r.index[r.used_size] = v2.index[p2]; 01139 r.used_size++; 01140 } 01141 } 01142 r.compact(); 01143 01144 return r; 01145 } 01146 01147 template <class T> 01148 Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01149 { 01150 it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> + Sparse_Vec<T>"); 01151 01152 Sparse_Vec<T> r(v1); 01153 ivec pos(v1.v_size); 01154 pos = -1; 01155 for (int p1=0; p1<v1.used_size; p1++) 01156 pos[v1.index[p1]] = p1; 01157 for (int p2=0; p2<v2.used_size; p2++) { 01158 if (pos[v2.index[p2]] == -1) {// A new entry 01159 if (r.used_size == r.data_size) 01160 r.resize_data(r.used_size*2+100); 01161 r.data[r.used_size] = v2.data[p2]; 01162 r.index[r.used_size] = v2.index[p2]; 01163 r.used_size++; 01164 } 01165 else 01166 r.data[pos[v2.index[p2]]] += v2.data[p2]; 01167 } 01168 r.compact(); 01169 01170 return r; 01171 } 01172 01173 template <class T> 01174 inline Sparse_Vec<T> sparse(const Vec<T> &v) 01175 { 01176 Sparse_Vec<T> s(v); 01177 return s; 01178 } 01179 01180 template <class T> 01181 inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon) 01182 { 01183 Sparse_Vec<T> s(v, epsilon); 01184 return s; 01185 } 01186 01187 template <class T> 01188 inline Vec<T> full(const Sparse_Vec<T> &s) 01189 { 01190 Vec<T> v; 01191 s.full(v); 01192 return v; 01193 } 01194 01195 /*-------------------------------------------------------------- 01196 * Explicit initializations 01197 *--------------------------------------------------------------*/ 01198 #ifndef _MSC_VER 01199 01201 extern template class Sparse_Vec<int>; 01203 extern template class Sparse_Vec<double>; 01205 extern template class Sparse_Vec<std::complex<double> >; 01206 01207 01209 extern template sparse_ivec operator+(const sparse_ivec &, const sparse_ivec &); 01211 extern template sparse_vec operator+(const sparse_vec &, const sparse_vec &); 01213 extern template sparse_cvec operator+(const sparse_cvec &, const sparse_cvec &); 01214 01216 extern template int operator*(const sparse_ivec &, const sparse_ivec &); 01218 extern template double operator*(const sparse_vec &, const sparse_vec &); 01220 extern template std::complex<double> operator*(const sparse_cvec &, const sparse_cvec &); 01221 01223 extern template int operator*(const sparse_ivec &, const ivec &); 01225 extern template double operator*(const sparse_vec &, const vec &); 01227 extern template std::complex<double> operator*(const sparse_cvec &, const cvec &); 01228 01230 extern template int operator*(const ivec &, const sparse_ivec &); 01232 extern template double operator*(const vec &, const sparse_vec &); 01234 extern template std::complex<double> operator*(const cvec &, const sparse_cvec &); 01235 01237 extern template sparse_ivec elem_mult(const sparse_ivec &, const sparse_ivec &); 01239 extern template sparse_vec elem_mult(const sparse_vec &, const sparse_vec &); 01241 extern template sparse_cvec elem_mult(const sparse_cvec &, const sparse_cvec &); 01242 01244 extern template ivec elem_mult(const sparse_ivec &, const ivec &); 01246 extern template vec elem_mult(const sparse_vec &, const vec &); 01248 extern template cvec elem_mult(const sparse_cvec &, const cvec &); 01249 01251 extern template sparse_ivec elem_mult_s(const sparse_ivec &, const ivec &); 01253 extern template sparse_vec elem_mult_s(const sparse_vec &, const vec &); 01255 extern template sparse_cvec elem_mult_s(const sparse_cvec &, const cvec &); 01256 01258 extern template ivec elem_mult(const ivec &, const sparse_ivec &); 01260 extern template vec elem_mult(const vec &, const sparse_vec &); 01262 extern template cvec elem_mult(const cvec &, const sparse_cvec &); 01263 01265 extern template sparse_ivec elem_mult_s(const ivec &, const sparse_ivec &); 01267 extern template sparse_vec elem_mult_s(const vec &, const sparse_vec &); 01269 extern template sparse_cvec elem_mult_s(const cvec &, const sparse_cvec &); 01270 #endif 01271 01272 } // namespace itpp 01273 01274 #endif // #ifndef SVEC_H 01275
Generated on Sat Aug 25 23:40:25 2007 for IT++ by Doxygen 1.5.2