00001 00033 #ifndef SVEC_H 00034 #define SVEC_H 00035 00036 #include <itpp/base/vec.h> 00037 #include <itpp/base/stat.h> 00038 00039 00040 namespace itpp { 00041 00042 // Declaration of class Vec 00043 template <class T> class Vec; 00044 // Declaration of class Sparse_Vec 00045 template <class T> class Sparse_Vec; 00046 00047 // ----------------------- Sparse_Vec Friends ------------------------------- 00048 00050 template <class T> 00051 Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00052 00054 template <class T> 00055 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 Vec<T> &v2); 00060 00062 template <class T> 00063 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00064 00066 template <class T> 00067 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00068 00070 template <class T> 00071 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00072 00074 template <class T> 00075 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00076 00078 template <class T> 00079 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00080 00082 template <class T> 00083 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00084 00093 template <class T> 00094 class Sparse_Vec { 00095 public: 00096 00098 Sparse_Vec(); 00099 00106 Sparse_Vec(int sz, int data_init=200); 00107 00113 Sparse_Vec(const Sparse_Vec<T> &v); 00114 00120 Sparse_Vec(const Vec<T> &v); 00121 00127 Sparse_Vec(const Vec<T> &v, T epsilon); 00128 00130 ~Sparse_Vec(); 00131 00138 void set_size(int sz, int data_init=-1); 00139 00141 int size() const { return v_size; } 00142 00144 inline int nnz() 00145 { 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 { 00219 if (check_small_elems_flag) { 00220 remove_small_elements(); 00221 } 00222 data_out = data[p]; 00223 } 00224 00226 inline T get_nz_data(int p) 00227 { 00228 if (check_small_elems_flag) { 00229 remove_small_elements(); 00230 } 00231 return data[p]; 00232 } 00233 00235 inline int get_nz_index(int p) 00236 { 00237 if (check_small_elems_flag) { 00238 remove_small_elements(); 00239 } 00240 return index[p]; 00241 } 00242 00244 inline void get_nz(int p, int &idx, T &dat) 00245 { 00246 if (check_small_elems_flag) { 00247 remove_small_elements(); 00248 } 00249 idx=index[p]; 00250 dat=data[p]; 00251 } 00252 00254 Sparse_Vec<T> get_subvector(int i1, int i2) const; 00255 00257 T sqr() const; 00258 00260 void operator=(const Sparse_Vec<T> &v); 00262 void operator=(const Vec<T> &v); 00263 00265 Sparse_Vec<T> operator-() const; 00266 00268 bool operator==(const Sparse_Vec<T> &v); 00269 00271 void operator+=(const Sparse_Vec<T> &v); 00272 00274 void operator+=(const 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 T &v); 00284 00286 void operator/=(const T &v); 00287 00289 friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00291 friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00293 friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00295 friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00296 00298 friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2); 00299 00301 friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00302 00304 friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2); 00305 00307 friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00308 00310 friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2); 00311 00312 private: 00313 void init(); 00314 void alloc(); 00315 void free(); 00316 00317 int v_size, used_size, data_size; 00318 T *data; 00319 int *index; 00320 T eps; 00321 bool check_small_elems_flag; 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_assert0(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_assert0(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_assert0(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_assert0(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_assert0(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_assert0(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_assert0(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_assert0(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_assert0(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 Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const 00784 { 00785 it_assert0(v_size > i1 && v_size > i2 && i1<=i2 && i1>=0, "The index of the element exceeds the size of the sparse vector"); 00786 00787 Sparse_Vec<T> r(i2-i1+1); 00788 00789 for (int p=0; p<used_size; p++) { 00790 if (index[p] >= i1 && index[p] <= i2) { 00791 if (r.used_size == r.data_size) 00792 r.resize_data(r.data_size*2+100); 00793 r.data[r.used_size] = data[p]; 00794 r.index[r.used_size] = index[p]-i1; 00795 r.used_size++; 00796 } 00797 } 00798 r.eps = eps; 00799 r.check_small_elems_flag = check_small_elems_flag; 00800 r.compact(); 00801 00802 return r; 00803 } 00804 00805 template <class T> 00806 T Sparse_Vec<T>::sqr() const 00807 { 00808 T sum(0); 00809 for (int p=0; p<used_size; p++) 00810 sum += data[p] * data[p]; 00811 00812 return sum; 00813 } 00814 00815 template <class T> 00816 void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v) 00817 { 00818 free(); 00819 v_size = v.v_size; 00820 used_size = v.used_size; 00821 data_size = v.data_size; 00822 eps = v.eps; 00823 check_small_elems_flag = v.check_small_elems_flag; 00824 alloc(); 00825 00826 for (int i=0; i<used_size; i++) { 00827 data[i] = v.data[i]; 00828 index[i] = v.index[i]; 00829 } 00830 } 00831 00832 template <class T> 00833 void Sparse_Vec<T>::operator=(const Vec<T> &v) 00834 { 00835 free(); 00836 v_size = v.size(); 00837 used_size = 0; 00838 data_size = std::min(v.size(), 10000); 00839 eps = T(0); 00840 check_small_elems_flag = false; 00841 alloc(); 00842 00843 for (int i=0; i<v_size; i++) { 00844 if (v(i) != T(0)) { 00845 if (used_size == data_size) 00846 resize_data(data_size*2); 00847 data[used_size] = v(i); 00848 index[used_size] = i; 00849 used_size++; 00850 } 00851 } 00852 compact(); 00853 } 00854 00855 template <class T> 00856 Sparse_Vec<T> Sparse_Vec<T>::operator-() const 00857 { 00858 Sparse_Vec r(v_size, used_size); 00859 00860 for (int p=0; p<used_size; p++) { 00861 r.data[p] = -data[p]; 00862 r.index[p] = index[p]; 00863 } 00864 r.used_size = used_size; 00865 00866 return r; 00867 } 00868 00869 template <class T> 00870 bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v) 00871 { 00872 int p,q; 00873 bool found=false; 00874 00875 //Remove small elements before comparing the two sparse_vectors 00876 if (check_small_elems_flag) 00877 remove_small_elements(); 00878 00879 if (v_size!=v.v_size) { 00880 //Return false if vector sizes are unequal 00881 return false; 00882 } else { 00883 for(p=0;p<used_size;p++) { 00884 for(q=0;q<v.used_size;q++) { 00885 if (index[p] == v.index[q]) { 00886 found = true; 00887 break; 00888 } 00889 } 00890 if (found==false) 00891 //Return false if non-zero element not found, or if elements are unequal 00892 return false; 00893 else if (data[p]!=v.data[q]) 00894 //Return false if non-zero element not found, or if elements are unequal 00895 return false; 00896 else 00897 //Check next non-zero element 00898 found = false; 00899 } 00900 } 00901 00902 /*Special handling if sizes do not match. 00903 Required since v may need to do remove_small_elements() for true comparison*/ 00904 if (used_size!=v.used_size) { 00905 if (used_size > v.used_size) { 00906 //Return false if number of non-zero elements is less in v 00907 return false; 00908 } 00909 else { 00910 //Ensure that the remaining non-zero elements in v are smaller than v.eps 00911 int nrof_small_elems = 0; 00912 for(q=0;q<v.used_size;q++) { 00913 if (std::abs(v.data[q]) <= std::abs(v.eps)) 00914 nrof_small_elems++; 00915 } 00916 if (v.used_size-nrof_small_elems != used_size) 00917 //Return false if the number of "true" non-zero elements are unequal 00918 return false; 00919 } 00920 } 00921 00922 //All elements checks => return true 00923 return true; 00924 } 00925 00926 template <class T> 00927 void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v) 00928 { 00929 int i,p; 00930 T tmp_data; 00931 int nrof_nz_v=v.used_size; 00932 00933 it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00934 00935 for (p=0; p<nrof_nz_v; p++) { 00936 i = v.index[p]; 00937 tmp_data = v.data[p]; 00938 //get_nz(p,i,tmp_data); 00939 add_elem(i,tmp_data); 00940 } 00941 00942 check_small_elems_flag = true; 00943 } 00944 00945 template <class T> 00946 void Sparse_Vec<T>::operator+=(const Vec<T> &v) 00947 { 00948 int i; 00949 00950 it_assert0(v_size == v.size(), "Attempted addition of unequal sized sparse vectors"); 00951 00952 for (i=0; i<v.size(); i++) 00953 if (v(i)!=T(0)) 00954 add_elem(i,v(i)); 00955 00956 check_small_elems_flag = true; 00957 } 00958 00959 00960 template <class T> 00961 void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v) 00962 { 00963 int i,p; 00964 T tmp_data; 00965 int nrof_nz_v=v.used_size; 00966 00967 it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 00968 00969 for (p=0; p<nrof_nz_v; p++) { 00970 i = v.index[p]; 00971 tmp_data = v.data[p]; 00972 //v.get_nz(p,i,tmp_data); 00973 add_elem(i,-tmp_data); 00974 } 00975 00976 check_small_elems_flag = true; 00977 } 00978 00979 template <class T> 00980 void Sparse_Vec<T>::operator-=(const Vec<T> &v) 00981 { 00982 int i; 00983 00984 it_assert0(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors"); 00985 00986 for (i=0; i<v.size(); i++) 00987 if (v(i)!=T(0)) 00988 add_elem(i,-v(i)); 00989 00990 check_small_elems_flag = true; 00991 } 00992 00993 template <class T> 00994 void Sparse_Vec<T>::operator*=(const T &v) 00995 { 00996 int p; 00997 00998 for (p=0; p<used_size; p++) { 00999 data[p]*=v;} 01000 01001 check_small_elems_flag = true; 01002 } 01003 01004 template <class T> 01005 void Sparse_Vec<T>::operator/=(const T &v) 01006 { 01007 int p; 01008 for (p=0; p<used_size; p++) { 01009 data[p]/=v;} 01010 01011 if (std::abs(eps) > 0) { 01012 check_small_elems_flag = true; 01013 } 01014 } 01015 01016 template <class T> 01017 T operator*(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01018 { 01019 it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> * Sparse_Vec<T>"); 01020 01021 T sum(0); 01022 Vec<T> v1f(v1.v_size); 01023 v1.full(v1f); 01024 for (int p=0; p<v2.used_size; p++) { 01025 if (v1f[v2.index[p]] != T(0)) 01026 sum += v1f[v2.index[p]] * v2.data[p]; 01027 } 01028 01029 return sum; 01030 } 01031 01032 template <class T> 01033 T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01034 { 01035 it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted"); 01036 01037 T sum(0); 01038 for (int p1=0; p1<v1.used_size; p1++) 01039 sum += v1.data[p1] * v2[v1.index[p1]]; 01040 01041 return sum; 01042 } 01043 01044 template <class T> 01045 T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01046 { 01047 it_assert0(v1.size() == v2.size(), "Multiplication of unequal sized vectors attempted"); 01048 01049 T sum(0); 01050 for (int p2=0; p2<v2.used_size; p2++) 01051 sum += v1[v2.index[p2]] * v2.data[p2]; 01052 01053 return sum; 01054 } 01055 01056 template <class T> 01057 Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01058 { 01059 it_assert0(v1.v_size == v2.v_size, "elem_mult(Sparse_Vec<T>, Sparse_Vec<T>)"); 01060 01061 Sparse_Vec<T> r(v1.v_size); 01062 ivec pos(v1.v_size); 01063 pos = -1; 01064 for (int p1=0; p1<v1.used_size; p1++) 01065 pos[v1.index[p1]] = p1; 01066 for (int p2=0; p2<v2.used_size; p2++) { 01067 if (pos[v2.index[p2]] != -1) { 01068 if (r.used_size == r.data_size) 01069 r.resize_data(r.used_size*2+100); 01070 r.data[r.used_size] = v1.data[pos[v2.index[p2]]] * v2.data[p2]; 01071 r.index[r.used_size] = v2.index[p2]; 01072 r.used_size++; 01073 } 01074 } 01075 r.compact(); 01076 01077 return r; 01078 } 01079 01080 template <class T> 01081 Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01082 { 01083 it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)"); 01084 01085 Vec<T> r(v1.v_size); 01086 r = T(0); 01087 for (int p1=0; p1<v1.used_size; p1++) 01088 r[v1.index[p1]] = v1.data[p1] * v2[v1.index[p1]]; 01089 01090 return r; 01091 } 01092 01093 template <class T> 01094 Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2) 01095 { 01096 it_assert0(v1.v_size == v2.size(), "elem_mult(Sparse_Vec<T>, Vec<T>)"); 01097 01098 Sparse_Vec<T> r(v1.v_size); 01099 for (int p1=0; p1<v1.used_size; p1++) { 01100 if (v2[v1.index[p1]] != T(0)) { 01101 if (r.used_size == r.data_size) 01102 r.resize_data(r.used_size*2+100); 01103 r.data[r.used_size] = v1.data[p1] * v2[v1.index[p1]]; 01104 r.index[r.used_size] = v1.index[p1]; 01105 r.used_size++; 01106 } 01107 } 01108 r.compact(); 01109 01110 return r; 01111 } 01112 01113 template <class T> 01114 Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01115 { 01116 it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)"); 01117 01118 Vec<T> r(v2.v_size); 01119 r = T(0); 01120 for (int p2=0; p2<v2.used_size; p2++) 01121 r[v2.index[p2]] = v1[v2.index[p2]] * v2.data[p2]; 01122 01123 return r; 01124 } 01125 01126 template <class T> 01127 Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2) 01128 { 01129 it_assert0(v1.size() == v2.v_size, "elem_mult(Vec<T>, Sparse_Vec<T>)"); 01130 01131 Sparse_Vec<T> r(v2.v_size); 01132 for (int p2=0; p2<v2.used_size; p2++) { 01133 if (v1[v2.index[p2]] != T(0)) { 01134 if (r.used_size == r.data_size) 01135 r.resize_data(r.used_size*2+100); 01136 r.data[r.used_size] = v1[v2.index[p2]] * v2.data[p2]; 01137 r.index[r.used_size] = v2.index[p2]; 01138 r.used_size++; 01139 } 01140 } 01141 r.compact(); 01142 01143 return r; 01144 } 01145 01146 template <class T> 01147 Sparse_Vec<T> operator+(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2) 01148 { 01149 it_assert0(v1.v_size == v2.v_size, "Sparse_Vec<T> + Sparse_Vec<T>"); 01150 01151 Sparse_Vec<T> r(v1); 01152 ivec pos(v1.v_size); 01153 pos = -1; 01154 for (int p1=0; p1<v1.used_size; p1++) 01155 pos[v1.index[p1]] = p1; 01156 for (int p2=0; p2<v2.used_size; p2++) { 01157 if (pos[v2.index[p2]] == -1) {// A new entry 01158 if (r.used_size == r.data_size) 01159 r.resize_data(r.used_size*2+100); 01160 r.data[r.used_size] = v2.data[p2]; 01161 r.index[r.used_size] = v2.index[p2]; 01162 r.used_size++; 01163 } 01164 else 01165 r.data[pos[v2.index[p2]]] += v2.data[p2]; 01166 } 01167 r.compact(); 01168 01169 return r; 01170 } 01171 01172 template <class T> 01173 inline Sparse_Vec<T> sparse(const Vec<T> &v) 01174 { 01175 Sparse_Vec<T> s(v); 01176 return s; 01177 } 01178 01179 template <class T> 01180 inline Sparse_Vec<T> sparse(const Vec<T> &v, T epsilon) 01181 { 01182 Sparse_Vec<T> s(v, epsilon); 01183 return s; 01184 } 01185 01186 template <class T> 01187 inline Vec<T> full(const Sparse_Vec<T> &s) 01188 { 01189 Vec<T> v; 01190 s.full(v); 01191 return v; 01192 } 01193 01194 /*-------------------------------------------------------------- 01195 * Explicit initializations 01196 *--------------------------------------------------------------*/ 01197 #ifndef _MSC_VER 01198 01200 extern template class Sparse_Vec<int>; 01202 extern template class Sparse_Vec<double>; 01204 extern template class Sparse_Vec<std::complex<double> >; 01205 01206 01208 extern template sparse_ivec operator+(const sparse_ivec &, const sparse_ivec &); 01210 extern template sparse_vec operator+(const sparse_vec &, const sparse_vec &); 01212 extern template sparse_cvec operator+(const sparse_cvec &, const sparse_cvec &); 01213 01215 extern template int operator*(const sparse_ivec &, const sparse_ivec &); 01217 extern template double operator*(const sparse_vec &, const sparse_vec &); 01219 extern template std::complex<double> operator*(const sparse_cvec &, const sparse_cvec &); 01220 01222 extern template int operator*(const sparse_ivec &, const ivec &); 01224 extern template double operator*(const sparse_vec &, const vec &); 01226 extern template std::complex<double> operator*(const sparse_cvec &, const cvec &); 01227 01229 extern template int operator*(const ivec &, const sparse_ivec &); 01231 extern template double operator*(const vec &, const sparse_vec &); 01233 extern template std::complex<double> operator*(const cvec &, const sparse_cvec &); 01234 01236 extern template sparse_ivec elem_mult(const sparse_ivec &, const sparse_ivec &); 01238 extern template sparse_vec elem_mult(const sparse_vec &, const sparse_vec &); 01240 extern template sparse_cvec elem_mult(const sparse_cvec &, const sparse_cvec &); 01241 01243 extern template ivec elem_mult(const sparse_ivec &, const ivec &); 01245 extern template vec elem_mult(const sparse_vec &, const vec &); 01247 extern template cvec elem_mult(const sparse_cvec &, const cvec &); 01248 01250 extern template sparse_ivec elem_mult_s(const sparse_ivec &, const ivec &); 01252 extern template sparse_vec elem_mult_s(const sparse_vec &, const vec &); 01254 extern template sparse_cvec elem_mult_s(const sparse_cvec &, const cvec &); 01255 01257 extern template ivec elem_mult(const ivec &, const sparse_ivec &); 01259 extern template vec elem_mult(const vec &, const sparse_vec &); 01261 extern template cvec elem_mult(const cvec &, const sparse_cvec &); 01262 01264 extern template sparse_ivec elem_mult_s(const ivec &, const sparse_ivec &); 01266 extern template sparse_vec elem_mult_s(const vec &, const sparse_vec &); 01268 extern template sparse_cvec elem_mult_s(const cvec &, const sparse_cvec &); 01269 #endif 01270 01271 } // namespace itpp 01272 01273 #endif // #ifndef SVEC_H 01274
Generated on Thu Apr 19 14:20:36 2007 for IT++ by Doxygen 1.4.6