IT++ Logo
svec.h
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
SourceForge Logo

Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4