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