IT++ Logo Newcom Logo

svec.h

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

Generated on Sat Aug 25 23:40:25 2007 for IT++ by Doxygen 1.5.2