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 00047 // Declaration of class Vec 00048 template <class T> class Vec; 00049 // Declaration of class Sparse_Vec 00050 template <class T> class Sparse_Vec; 00051 00052 // ----------------------- Sparse_Vec Friends ------------------------------- 00053 00055 template <class T> 00056 Sparse_Vec<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 Sparse_Vec<T> &v2); 00061 00063 template <class T> 00064 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00065 00067 template <class T> 00068 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00069 00071 template <class T> 00072 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00073 00075 template <class T> 00076 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00077 00079 template <class T> 00080 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00081 00083 template <class T> 00084 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00085 00087 template <class T> 00088 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00089 00100 template <class T> 00101 class Sparse_Vec 00102 { 00103 public: 00104 00106 Sparse_Vec(); 00107 00114 Sparse_Vec(int sz, int data_init = 200); 00115 00121 Sparse_Vec(const Sparse_Vec<T> &v); 00122 00128 Sparse_Vec(const Vec<T> &v); 00129 00135 Sparse_Vec(const Vec<T> &v, T epsilon); 00136 00138 ~Sparse_Vec(); 00139 00146 void set_size(int sz, int data_init = -1); 00147 00149 int size() const { return v_size; } 00150 00152 inline int nnz() { 00153 if (check_small_elems_flag) { 00154 remove_small_elements(); 00155 } 00156 return used_size; 00157 } 00158 00160 double density(); 00161 00163 void set_small_element(const T& epsilon); 00164 00170 void remove_small_elements(); 00171 00177 void resize_data(int new_size); 00178 00180 void compact(); 00181 00183 void full(Vec<T> &v) const; 00184 00186 Vec<T> full() const; 00187 00189 T operator()(int i) const; 00190 00192 void set(int i, T v); 00193 00195 void set(const ivec &index_vec, const Vec<T> &v); 00196 00198 void set_new(int i, T v); 00199 00201 void set_new(const ivec &index_vec, const Vec<T> &v); 00202 00204 void add_elem(const int i, const T v); 00205 00207 void add(const ivec &index_vec, const Vec<T> &v); 00208 00210 void zeros(); 00211 00213 void zero_elem(const int i); 00214 00216 void clear(); 00217 00219 void clear_elem(const int i); 00220 00224 inline void get_nz_data(int p, T& data_out) { 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 if (check_small_elems_flag) { 00234 remove_small_elements(); 00235 } 00236 return data[p]; 00237 } 00238 00240 inline int get_nz_index(int p) { 00241 if (check_small_elems_flag) { 00242 remove_small_elements(); 00243 } 00244 return index[p]; 00245 } 00246 00248 inline void get_nz(int p, int &idx, T &dat) { 00249 if (check_small_elems_flag) { 00250 remove_small_elements(); 00251 } 00252 idx = index[p]; 00253 dat = data[p]; 00254 } 00255 00257 ivec get_nz_indices(); 00258 00260 Sparse_Vec<T> get_subvector(int i1, int i2) const; 00261 00263 T sqr() const; 00264 00266 void operator=(const Sparse_Vec<T> &v); 00268 void operator=(const Vec<T> &v); 00269 00271 Sparse_Vec<T> operator-() const; 00272 00274 bool operator==(const Sparse_Vec<T> &v); 00275 00277 void operator+=(const Sparse_Vec<T> &v); 00278 00280 void operator+=(const Vec<T> &v); 00281 00283 void operator-=(const Sparse_Vec<T> &v); 00284 00286 void operator-=(const Vec<T> &v); 00287 00289 void operator*=(const T &v); 00290 00292 void operator/=(const T &v); 00293 00295 friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00297 friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00299 friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00301 friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00302 00304 friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00305 00307 friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00308 00310 friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00311 00313 friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00314 00316 friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00317 00318 private: 00319 void init(); 00320 void alloc(); 00321 void free(); 00322 00323 int v_size, used_size, data_size; 00324 T *data; 00325 int *index; 00326 T eps; 00327 bool check_small_elems_flag; 00328 }; 00329 00330 00335 typedef Sparse_Vec<int> sparse_ivec; 00336 00341 typedef Sparse_Vec<double> sparse_vec; 00342 00347 typedef Sparse_Vec<std::complex<double> > sparse_cvec; 00348 00349 // ----------------------- Implementation starts here -------------------------------- 00350 00351 template <class T> 00352 void Sparse_Vec<T>::init() 00353 { 00354 v_size = 0; 00355 used_size = 0; 00356 data_size = 0; 00357 data = 0; 00358 index = 0; 00359 eps = 0; 00360 check_small_elems_flag = true; 00361 } 00362 00363 template <class T> 00364 void Sparse_Vec<T>::alloc() 00365 { 00366 if (data_size != 0) { 00367 data = new T[data_size]; 00368 index = new int[data_size]; 00369 } 00370 } 00371 00372 template <class T> 00373 void Sparse_Vec<T>::free() 00374 { 00375 delete [] data; 00376 data = 0; 00377 delete [] index; 00378 index = 0; 00379 } 00380 00381 template <class T> 00382 Sparse_Vec<T>::Sparse_Vec() 00383 { 00384 init(); 00385 } 00386 00387 template <class T> 00388 Sparse_Vec<T>::Sparse_Vec(int sz, int data_init) 00389 { 00390 init(); 00391 v_size = sz; 00392 used_size = 0; 00393 data_size = data_init; 00394 alloc(); 00395 } 00396 00397 template <class T> 00398 Sparse_Vec<T>::Sparse_Vec(const Sparse_Vec<T> &v) 00399 { 00400 init(); 00401 v_size = v.v_size; 00402 used_size = v.used_size; 00403 data_size = v.data_size; 00404 eps = v.eps; 00405 check_small_elems_flag = v.check_small_elems_flag; 00406 alloc(); 00407 00408 for (int i = 0; i < used_size; i++) { 00409 data[i] = v.data[i]; 00410 index[i] = v.index[i]; 00411 } 00412 } 00413 00414 template <class T> 00415 Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v) 00416 { 00417 init(); 00418 v_size = v.size(); 00419 used_size = 0; 00420 data_size = std::min(v.size(), 10000); 00421 alloc(); 00422 00423 for (int i = 0; i < v_size; i++) { 00424 if (v(i) != T(0)) { 00425 if (used_size == data_size) 00426 resize_data(data_size*2); 00427 data[used_size] = v(i); 00428 index[used_size] = i; 00429 used_size++; 00430 } 00431 } 00432 compact(); 00433 } 00434 00435 template <class T> 00436 Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon) 00437 { 00438 init(); 00439 v_size = v.size(); 00440 used_size = 0; 00441 data_size = std::min(v.size(), 10000); 00442 eps = epsilon; 00443 alloc(); 00444 00445 double e = std::abs(epsilon); 00446 for (int i = 0; i < v_size; i++) { 00447 if (std::abs(v(i)) > e) { 00448 if (used_size == data_size) 00449 resize_data(data_size*2); 00450 data[used_size] = v(i); 00451 index[used_size] = i; 00452 used_size++; 00453 } 00454 } 00455 compact(); 00456 } 00457 00458 template <class T> 00459 Sparse_Vec<T>::~Sparse_Vec() 00460 { 00461 free(); 00462 } 00463 00464 template <class T> 00465 void Sparse_Vec<T>::set_size(int new_size, int data_init) 00466 { 00467 v_size = new_size; 00468 used_size = 0; 00469 if (data_init != -1) { 00470 free(); 00471 data_size = data_init; 00472 alloc(); 00473 } 00474 } 00475 00476 template <class T> 00477 double Sparse_Vec<T>::density() 00478 { 00479 if (check_small_elems_flag) { 00480 remove_small_elements(); 00481 } 00482 //return static_cast<double>(used_size) / v_size; 00483 return double(used_size) / v_size; 00484 } 00485 00486 template <class T> 00487 void Sparse_Vec<T>::set_small_element(const T& epsilon) 00488 { 00489 eps = epsilon; 00490 remove_small_elements(); 00491 } 00492 00493 template <class T> 00494 void Sparse_Vec<T>::remove_small_elements() 00495 { 00496 int i; 00497 int nrof_removed_elements = 0; 00498 double e; 00499 00500 //Remove small elements 00501 e = std::abs(eps); 00502 for (i = 0;i < used_size;i++) { 00503 if (std::abs(data[i]) <= e) { 00504 nrof_removed_elements++; 00505 } 00506 else if (nrof_removed_elements > 0) { 00507 data[i-nrof_removed_elements] = data[i]; 00508 index[i-nrof_removed_elements] = index[i]; 00509 } 00510 } 00511 00512 //Set new size after small elements have been removed 00513 used_size -= nrof_removed_elements; 00514 00515 //Set the flag to indicate that all small elements have been removed 00516 check_small_elems_flag = false; 00517 } 00518 00519 00520 template <class T> 00521 void Sparse_Vec<T>::resize_data(int new_size) 00522 { 00523 it_assert(new_size >= used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small"); 00524 00525 if (new_size != data_size) { 00526 if (new_size == 0) 00527 free(); 00528 else { 00529 T *tmp_data = data; 00530 int *tmp_pos = index; 00531 data_size = new_size; 00532 alloc(); 00533 for (int p = 0; p < used_size; p++) { 00534 data[p] = tmp_data[p]; 00535 index[p] = tmp_pos[p]; 00536 } 00537 delete [] tmp_data; 00538 delete [] tmp_pos; 00539 } 00540 } 00541 } 00542 00543 template <class T> 00544 void Sparse_Vec<T>::compact() 00545 { 00546 if (check_small_elems_flag) { 00547 remove_small_elements(); 00548 } 00549 resize_data(used_size); 00550 } 00551 00552 template <class T> 00553 void Sparse_Vec<T>::full(Vec<T> &v) const 00554 { 00555 v.set_size(v_size); 00556 00557 v = T(0); 00558 for (int p = 0; p < used_size; p++) 00559 v(index[p]) = data[p]; 00560 } 00561 00562 template <class T> 00563 Vec<T> Sparse_Vec<T>::full() const 00564 { 00565 Vec<T> r(v_size); 00566 full(r); 00567 return r; 00568 } 00569 00570 // This is slow. Implement a better search 00571 template <class T> 00572 T Sparse_Vec<T>::operator()(int i) const 00573 { 00574 it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range"); 00575 00576 bool found = false; 00577 int p; 00578 for (p = 0; p < used_size; p++) { 00579 if (index[p] == i) { 00580 found = true; 00581 break; 00582 } 00583 } 00584 return found ? data[p] : T(0); 00585 } 00586 00587 template <class T> 00588 void Sparse_Vec<T>::set(int i, T v) 00589 { 00590 it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range"); 00591 00592 bool found = false; 00593 bool larger_than_eps; 00594 int p; 00595 00596 for (p = 0; p < used_size; p++) { 00597 if (index[p] == i) { 00598 found = true; 00599 break; 00600 } 00601 } 00602 00603 larger_than_eps = (std::abs(v) > std::abs(eps)); 00604 00605 if (found && larger_than_eps) 00606 data[p] = v; 00607 else if (larger_than_eps) { 00608 if (used_size == data_size) 00609 resize_data(data_size*2 + 100); 00610 data[used_size] = v; 00611 index[used_size] = i; 00612 used_size++; 00613 } 00614 00615 //Check if the stored element is smaller than eps. In that case it should be removed. 00616 if (std::abs(v) <= std::abs(eps)) { 00617 remove_small_elements(); 00618 } 00619 00620 } 00621 00622 template <class T> 00623 void Sparse_Vec<T>::set_new(int i, T v) 00624 { 00625 it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00626 00627 //Check that the new element is larger than eps! 00628 if (std::abs(v) > std::abs(eps)) { 00629 if (used_size == data_size) 00630 resize_data(data_size*2 + 100); 00631 data[used_size] = v; 00632 index[used_size] = i; 00633 used_size++; 00634 } 00635 } 00636 00637 template <class T> 00638 void Sparse_Vec<T>::add_elem(const int i, const T v) 00639 { 00640 bool found = false; 00641 int p; 00642 00643 it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00644 00645 for (p = 0; p < used_size; p++) { 00646 if (index[p] == i) { 00647 found = true; 00648 break; 00649 } 00650 } 00651 if (found) 00652 data[p] += v; 00653 else { 00654 if (used_size == data_size) 00655 resize_data(data_size*2 + 100); 00656 data[used_size] = v; 00657 index[used_size] = i; 00658 used_size++; 00659 } 00660 00661 check_small_elems_flag = true; 00662 00663 } 00664 00665 template <class T> 00666 void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v) 00667 { 00668 bool found = false; 00669 int i, p, q; 00670 int nrof_nz = v.size(); 00671 00672 it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector"); 00673 00674 //Elements are added if they have identical indices 00675 for (q = 0; q < nrof_nz; q++) { 00676 i = index_vec(q); 00677 for (p = 0; p < used_size; p++) { 00678 if (index[p] == i) { 00679 found = true; 00680 break; 00681 } 00682 } 00683 if (found) 00684 data[p] += v(q); 00685 else { 00686 if (used_size == data_size) 00687 resize_data(data_size*2 + 100); 00688 data[used_size] = v(q); 00689 index[used_size] = i; 00690 used_size++; 00691 } 00692 found = false; 00693 } 00694 00695 check_small_elems_flag = true; 00696 00697 } 00698 00699 template <class T> 00700 void Sparse_Vec<T>::zeros() 00701 { 00702 used_size = 0; 00703 check_small_elems_flag = false; 00704 } 00705 00706 template <class T> 00707 void Sparse_Vec<T>::zero_elem(const int i) 00708 { 00709 bool found = false; 00710 int p; 00711 00712 it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00713 00714 for (p = 0; p < used_size; p++) { 00715 if (index[p] == i) { 00716 found = true; 00717 break; 00718 } 00719 } 00720 if (found) { 00721 data[p] = data[used_size-1]; 00722 index[p] = index[used_size-1]; 00723 used_size--; 00724 } 00725 } 00726 00727 template <class T> 00728 void Sparse_Vec<T>::clear() 00729 { 00730 used_size = 0; 00731 check_small_elems_flag = false; 00732 } 00733 00734 template <class T> 00735 void Sparse_Vec<T>::clear_elem(const int i) 00736 { 00737 bool found = false; 00738 int p; 00739 00740 it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector"); 00741 00742 for (p = 0; p < used_size; p++) { 00743 if (index[p] == i) { 00744 found = true; 00745 break; 00746 } 00747 } 00748 if (found) { 00749 data[p] = data[used_size-1]; 00750 index[p] = index[used_size-1]; 00751 used_size--; 00752 } 00753 } 00754 00755 template <class T> 00756 void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v) 00757 { 00758 it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector"); 00759 00760 //Clear all old non-zero elements 00761 clear(); 00762 00763 //Add the new non-zero elements 00764 add(index_vec, v); 00765 } 00766 00767 template <class T> 00768 void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v) 00769 { 00770 int q; 00771 int nrof_nz = v.size(); 00772 00773 it_assert_debug(v_size > max(index_vec), "The indices exceeds the size of the sparse vector"); 00774 00775 //Clear all old non-zero elements 00776 clear(); 00777 00778 for (q = 0; q < nrof_nz; q++) { 00779 if (std::abs(v[q]) > std::abs(eps)) { 00780 if (used_size == data_size) 00781 resize_data(data_size*2 + 100); 00782 data[used_size] = v(q); 00783 index[used_size] = index_vec(q); 00784 used_size++; 00785 } 00786 } 00787 } 00788 00789 template <class T> 00790 ivec Sparse_Vec<T>::get_nz_indices() 00791 { 00792 int n = nnz(); 00793 ivec r(n); 00794 for (int i = 0; i < n; i++) { 00795 r(i) = get_nz_index(i); 00796 } 00797 return r; 00798 } 00799 00800 template <class T> 00801 Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const 00802 { 00803 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"); 00804 00805 Sparse_Vec<T> r(i2 - i1 + 1); 00806 00807 for (int p = 0; p < used_size; p++) { 00808 if (index[p] >= i1 && index[p] <= i2) { 00809 if (r.used_size == r.data_size) 00810 r.resize_data(r.data_size*2 + 100); 00811 r.data[r.used_size] = data[p]; 00812 r.index[r.used_size] = index[p] - i1; 00813 r.used_size++; 00814 } 00815 } 00816 r.eps = eps; 00817 r.check_small_elems_flag = check_small_elems_flag; 00818 r.compact(); 00819 00820 return r; 00821 } 00822 00823 template <class T> 00824 T Sparse_Vec<T>::sqr() const 00825 { 00826 T sum(0); 00827 for (int p = 0; p < used_size; p++) 00828 sum += data[p] * data[p]; 00829 00830 return sum; 00831 } 00832 00833 template <class T> 00834 void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v) 00835 { 00836 free(); 00837 v_size = v.v_size; 00838 used_size = v.used_size; 00839 data_size = v.data_size; 00840 eps = v.eps; 00841 check_small_elems_flag = v.check_small_elems_flag; 00842 alloc(); 00843 00844 for (int i = 0; i < used_size; i++) { 00845 data[i] = v.data[i]; 00846 index[i] = v.index[i]; 00847 } 00848 } 00849 00850 template <class T> 00851 void Sparse_Vec<T>::operator=(const Vec<T> &v) 00852 { 00853 free(); 00854 v_size = v.size(); 00855 used_size = 0; 00856 data_size = std::min(v.size(), 10000); 00857 eps = T(0); 00858 check_small_elems_flag = false; 00859 alloc(); 00860 00861 for (int i = 0; i < v_size; i++) { 00862 if (v(i) != T(0)) { 00863 if (used_size == data_size) 00864 resize_data(data_size*2); 00865 data[used_size] = v(i); 00866 index[used_size] = i; 00867 used_size++; 00868 } 00869 } 00870 compact(); 00871 } 00872 00873 template <class T> 00874 Sparse_Vec<T> Sparse_Vec<T>::operator-() const 00875 { 00876 Sparse_Vec r(v_size, used_size); 00877 00878 for (int p = 0; p < used_size; p++) { 00879 r.data[p] = -data[p]; 00880 r.index[p] = index[p]; 00881 } 00882 r.used_size = used_size; 00883 00884 return r; 00885 } 00886 00887 template <class T> 00888 bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v) 00889 { 00890 int p, q; 00891 bool found = false; 00892 00893 //Remove small elements before comparing the two sparse_vectors 00894 if (check_small_elems_flag) 00895 remove_small_elements(); 00896 00897 if (v_size != v.v_size) { 00898 //Return false if vector sizes are unequal 00899 return false; 00900 } 00901 else { 00902 for (p = 0;p < used_size;p++) { 00903 for (q = 0;q < v.used_size;q++) { 00904 if (index[p] == v.index[q]) { 00905 found = true; 00906 break; 00907 } 00908 } 00909 if (found == false) 00910 //Return false if non-zero element not found, or if elements are unequal 00911 return false; 00912 else if (data[p] != v.data[q]) 00913 //Return false if non-zero element not found, or if elements are unequal 00914 return false; 00915 else 00916 //Check next non-zero element 00917 found = false; 00918 } 00919 } 00920 00921 /*Special handling if sizes do not match. 00922 Required since v may need to do remove_small_elements() for true comparison*/ 00923 if (used_size != v.used_size) { 00924 if (used_size > v.used_size) { 00925 //Return false if number of non-zero elements is less in v 00926 return false; 00927 } 00928 else { 00929 //Ensure that the remaining non-zero elements in v are smaller than v.eps 00930 int nrof_small_elems = 0; 00931 for (q = 0;q < v.used_size;q++) { 00932 if (std::abs(v.data[q]) <= std::abs(v.eps)) 00933 nrof_small_elems++; 00934 } 00935 if (v.used_size - nrof_small_elems != used_size) 00936 //Return false if the number of "true" non-zero elements are unequal 00937 return false; 00938 } 00939 } 00940 00941 //All elements checks => return true 00942 return true; 00943 } 00944 00945 template <class T> 00946 void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v) 00947 { 00948 int i, p; 00949 T tmp_data; 00950 int nrof_nz_v = v.used_size; 00951 00952 it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00953 00954 for (p = 0; p < nrof_nz_v; p++) { 00955 i = v.index[p]; 00956 tmp_data = v.data[p]; 00957 //get_nz(p,i,tmp_data); 00958 add_elem(i, tmp_data); 00959 } 00960 00961 check_small_elems_flag = true; 00962 } 00963 00964 template <class T> 00965 void Sparse_Vec<T>::operator+=(const Vec<T> &v) 00966 { 00967 int i; 00968 00969 it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00970 00971 for (i = 0; i < v.size(); i++) 00972 if (v(i) != T(0)) 00973 add_elem(i, v(i)); 00974 00975 check_small_elems_flag = true; 00976 } 00977 00978 00979 template <class T> 00980 void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v) 00981 { 00982 int i, p; 00983 T tmp_data; 00984 int nrof_nz_v = v.used_size; 00985 00986 it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 00987 00988 for (p = 0; p < nrof_nz_v; p++) { 00989 i = v.index[p]; 00990 tmp_data = v.data[p]; 00991 //v.get_nz(p,i,tmp_data); 00992 add_elem(i, -tmp_data); 00993 } 00994 00995 check_small_elems_flag = true; 00996 } 00997 00998 template <class T> 00999 void Sparse_Vec<T>::operator-=(const Vec<T> &v) 01000 { 01001 int i; 01002 01003 it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 01004 01005 for (i = 0; i < v.size(); i++) 01006 if (v(i) != T(0)) 01007 add_elem(i, -v(i)); 01008 01009 check_small_elems_flag = true; 01010 } 01011 01012 template <class T> 01013 void Sparse_Vec<T>::operator*=(const T &v) 01014 { 01015 int p; 01016 01017 for (p = 0; p < used_size; p++) { 01018 data[p] *= v; 01019 } 01020 01021 check_small_elems_flag = true; 01022 } 01023 01024 template <class T> 01025 void Sparse_Vec<T>::operator/=(const T &v) 01026 { 01027 int p; 01028 for (p = 0; p < used_size; p++) { 01029 data[p] /= v; 01030 } 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 Wed Feb 9 2011 13:47:23 for IT++ by Doxygen 1.7.3