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

Generated on Thu Apr 19 14:23:57 2007 for IT++ by Doxygen 1.4.6