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