IT++ Logo

svec.h

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

Generated on Tue Nov 23 08:45:10 2010 for IT++ by Doxygen 1.6.1