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   // Declaration of class Vec
00047   template <class T> class Vec;
00048   // Declaration of class Sparse_Vec
00049   template <class T> class Sparse_Vec;
00050 
00051   // ----------------------- Sparse_Vec Friends -------------------------------
00052 
00054   template <class T>
00055     Sparse_Vec<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 Sparse_Vec<T> &v2);
00060 
00062   template <class T>
00063     T operator*(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00064 
00066   template <class T>
00067     T operator*(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00068 
00070   template <class T>
00071     Sparse_Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
00072 
00074   template <class T>
00075     Vec<T> elem_mult(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00076 
00078   template <class T>
00079     Sparse_Vec<T> elem_mult_s(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00080 
00082   template <class T>
00083     Vec<T> elem_mult(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00084 
00086   template <class T>
00087     Sparse_Vec<T> elem_mult_s(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00088 
00099   template <class T>
00100     class Sparse_Vec {
00101     public:
00102 
00104     Sparse_Vec();
00105 
00112     Sparse_Vec(int sz, int data_init=200);
00113 
00119     Sparse_Vec(const Sparse_Vec<T> &v);
00120 
00126     Sparse_Vec(const Vec<T> &v);
00127 
00133     Sparse_Vec(const Vec<T> &v, T epsilon);
00134 
00136     ~Sparse_Vec();
00137 
00144     void set_size(int sz, int data_init=-1);
00145 
00147     int size() const { return v_size; }
00148 
00150     inline int nnz()
00151       {
00152         if (check_small_elems_flag) {
00153           remove_small_elements();
00154         }
00155         return used_size;
00156       }
00157 
00159     double density();
00160 
00162     void set_small_element(const T& epsilon);
00163 
00169     void remove_small_elements();
00170 
00176     void resize_data(int new_size);
00177 
00179     void compact();
00180 
00182     void full(Vec<T> &v) const;
00183 
00185     Vec<T> full() const;
00186 
00188     T operator()(int i) const;
00189 
00191     void set(int i, T v);
00192 
00194     void set(const ivec &index_vec, const Vec<T> &v);
00195 
00197     void set_new(int i, T v);
00198 
00200     void set_new(const ivec &index_vec, const Vec<T> &v);
00201 
00203     void add_elem(const int i, const T v);
00204 
00206     void add(const ivec &index_vec, const Vec<T> &v);
00207 
00209     void zeros();
00210 
00212     void zero_elem(const int i);
00213 
00215     void clear();
00216 
00218     void clear_elem(const int i);
00219 
00223     inline void get_nz_data(int p, T& data_out)
00224       {
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       {
00234         if (check_small_elems_flag) {
00235           remove_small_elements();
00236         }
00237         return data[p];
00238       }
00239 
00241     inline int get_nz_index(int p)
00242       {
00243         if (check_small_elems_flag) {
00244           remove_small_elements();
00245         }
00246         return index[p];
00247       }
00248 
00250     inline void get_nz(int p, int &idx, T &dat)
00251       {
00252         if (check_small_elems_flag) {
00253           remove_small_elements();
00254         }
00255         idx=index[p];
00256         dat=data[p];
00257       }
00258 
00260     ivec get_nz_indices();
00261 
00263     Sparse_Vec<T> get_subvector(int i1, int i2) const;
00264 
00266     T sqr() const;
00267 
00269     void operator=(const Sparse_Vec<T> &v);
00271     void operator=(const Vec<T> &v);
00272 
00274     Sparse_Vec<T> operator-() const;
00275 
00277     bool operator==(const Sparse_Vec<T> &v);
00278 
00280     void operator+=(const Sparse_Vec<T> &v);
00281 
00283     void operator+=(const Vec<T> &v);
00284 
00286     void operator-=(const Sparse_Vec<T> &v);
00287 
00289     void operator-=(const Vec<T> &v);
00290 
00292     void operator*=(const T &v);
00293 
00295     void operator/=(const T &v);
00296 
00298     friend Sparse_Vec<T> operator+<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
00300     friend T operator*<>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
00302     friend T operator*<>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00304     friend T operator*<>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00305 
00307     friend Sparse_Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Sparse_Vec<T> &v2);
00308 
00310     friend Vec<T> elem_mult <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00311 
00313     friend Sparse_Vec<T> elem_mult_s <>(const Sparse_Vec<T> &v1, const Vec<T> &v2);
00314 
00316     friend Vec<T> elem_mult <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00317 
00319     friend Sparse_Vec<T> elem_mult_s <>(const Vec<T> &v1, const Sparse_Vec<T> &v2);
00320 
00321     private:
00322     void init();
00323     void alloc();
00324     void free();
00325 
00326     int v_size, used_size, data_size;
00327     T *data;
00328     int *index;
00329     T eps;
00330     bool check_small_elems_flag;
00331   };
00332 
00333 
00338   typedef Sparse_Vec<int> sparse_ivec;
00339 
00344   typedef Sparse_Vec<double> sparse_vec;
00345 
00350   typedef Sparse_Vec<std::complex<double> > sparse_cvec;
00351 
00352   // ----------------------- Implementation starts here --------------------------------
00353 
00354   template <class T>
00355     void Sparse_Vec<T>::init()
00356     {
00357       v_size = 0;
00358       used_size = 0;
00359       data_size = 0;
00360       data = 0;
00361       index = 0;
00362       eps = 0;
00363       check_small_elems_flag = true;
00364     }
00365 
00366   template <class T>
00367     void Sparse_Vec<T>::alloc()
00368     {
00369       if (data_size != 0) {
00370         data = new T[data_size];
00371         index = new int[data_size];
00372       }
00373     }
00374 
00375   template <class T>
00376     void Sparse_Vec<T>::free()
00377     {
00378       delete [] data;
00379       data = 0;
00380       delete [] index;
00381       index = 0;
00382     }
00383 
00384   template <class T>
00385     Sparse_Vec<T>::Sparse_Vec()
00386     {
00387       init();
00388     }
00389 
00390   template <class T>
00391     Sparse_Vec<T>::Sparse_Vec(int sz, int data_init)
00392     {
00393       init();
00394       v_size = sz;
00395       used_size = 0;
00396       data_size = data_init;
00397       alloc();
00398     }
00399 
00400   template <class T>
00401     Sparse_Vec<T>::Sparse_Vec(const Sparse_Vec<T> &v)
00402     {
00403       init();
00404       v_size = v.v_size;
00405       used_size = v.used_size;
00406       data_size = v.data_size;
00407       eps = v.eps;
00408       check_small_elems_flag = v.check_small_elems_flag;
00409       alloc();
00410 
00411       for (int i=0; i<used_size; i++) {
00412         data[i] = v.data[i];
00413         index[i] = v.index[i];
00414       }
00415     }
00416 
00417   template <class T>
00418     Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v)
00419     {
00420       init();
00421       v_size = v.size();
00422       used_size = 0;
00423       data_size = std::min(v.size(), 10000);
00424       alloc();
00425 
00426       for (int i=0; i<v_size; i++) {
00427         if (v(i) != T(0)) {
00428           if (used_size == data_size)
00429             resize_data(data_size*2);
00430           data[used_size] = v(i);
00431           index[used_size] = i;
00432           used_size++;
00433         }
00434       }
00435       compact();
00436     }
00437 
00438   template <class T>
00439     Sparse_Vec<T>::Sparse_Vec(const Vec<T> &v, T epsilon)
00440     {
00441       init();
00442       v_size = v.size();
00443       used_size = 0;
00444       data_size = std::min(v.size(), 10000);
00445       eps = epsilon;
00446       alloc();
00447 
00448     double e = std::abs(epsilon);
00449     for (int i=0; i<v_size; i++) {
00450         if (std::abs(v(i)) > e) {
00451           if (used_size == data_size)
00452             resize_data(data_size*2);
00453           data[used_size] = v(i);
00454           index[used_size] = i;
00455           used_size++;
00456         }
00457       }
00458       compact();
00459     }
00460 
00461   template <class T>
00462     Sparse_Vec<T>::~Sparse_Vec()
00463     {
00464       free();
00465     }
00466 
00467   template <class T>
00468     void Sparse_Vec<T>::set_size(int new_size, int data_init)
00469     {
00470       v_size = new_size;
00471       used_size = 0;
00472       if (data_init != -1){
00473         free();
00474         data_size = data_init;
00475         alloc();
00476       }
00477     }
00478 
00479   template <class T>
00480     double Sparse_Vec<T>::density()
00481     {
00482       if (check_small_elems_flag) {
00483         remove_small_elements();
00484       }
00485       //return static_cast<double>(used_size) / v_size;
00486       return double(used_size) / v_size;
00487     }
00488 
00489   template <class T>
00490     void Sparse_Vec<T>::set_small_element(const T& epsilon)
00491     {
00492       eps=epsilon;
00493       remove_small_elements();
00494     }
00495 
00496   template <class T>
00497     void Sparse_Vec<T>::remove_small_elements()
00498     {
00499       int i;
00500       int nrof_removed_elements = 0;
00501       double e;
00502 
00503       //Remove small elements
00504       e = std::abs(eps);
00505       for (i=0;i<used_size;i++) {
00506         if (std::abs(data[i]) <= e) {
00507           nrof_removed_elements++;
00508         }
00509         else if (nrof_removed_elements > 0) {
00510           data[i-nrof_removed_elements] = data[i];
00511           index[i-nrof_removed_elements] = index[i];
00512         }
00513       }
00514 
00515       //Set new size after small elements have been removed
00516       used_size -= nrof_removed_elements;
00517 
00518       //Set the flag to indicate that all small elements have been removed
00519       check_small_elems_flag = false;
00520     }
00521 
00522 
00523   template <class T>
00524     void Sparse_Vec<T>::resize_data(int new_size)
00525     {
00526       it_assert(new_size>=used_size, "Sparse_Vec<T>::resize_data(int new_size): New size is to small");
00527 
00528       if (new_size != data_size) {
00529         if (new_size == 0)
00530           free();
00531         else {
00532           T *tmp_data = data;
00533           int *tmp_pos = index;
00534           data_size = new_size;
00535           alloc();
00536           for (int p=0; p<used_size; p++) {
00537             data[p] = tmp_data[p];
00538             index[p] = tmp_pos[p];
00539           }
00540           delete [] tmp_data;
00541           delete [] tmp_pos;
00542         }
00543       }
00544     }
00545 
00546   template <class T>
00547     void Sparse_Vec<T>::compact()
00548     {
00549       if (check_small_elems_flag) {
00550         remove_small_elements();
00551       }
00552       resize_data(used_size);
00553     }
00554 
00555   template <class T>
00556     void Sparse_Vec<T>::full(Vec<T> &v) const
00557     {
00558       v.set_size(v_size);
00559 
00560       v = T(0);
00561       for (int p=0; p<used_size; p++)
00562         v(index[p]) = data[p];
00563     }
00564 
00565   template <class T>
00566     Vec<T> Sparse_Vec<T>::full() const
00567     {
00568       Vec<T> r(v_size);
00569       full(r);
00570       return r;
00571     }
00572 
00573   // This is slow. Implement a better search
00574   template <class T>
00575     T Sparse_Vec<T>::operator()(int i) const
00576     {
00577       it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");
00578 
00579       bool found = false;
00580       int p;
00581       for (p=0; p<used_size; p++) {
00582         if (index[p] == i) {
00583           found = true;
00584           break;
00585         }
00586       }
00587       return found ? data[p] : T(0);
00588     }
00589 
00590   template <class T>
00591     void Sparse_Vec<T>::set(int i, T v)
00592     {
00593       it_assert_debug(i >= 0 && i < v_size, "The index of the element is out of range");
00594 
00595       bool found = false;
00596       bool larger_than_eps;
00597       int p;
00598 
00599       for (p=0; p<used_size; p++) {
00600         if (index[p] == i) {
00601           found = true;
00602           break;
00603         }
00604       }
00605 
00606       larger_than_eps = (std::abs(v) > std::abs(eps));
00607 
00608       if (found && larger_than_eps)
00609         data[p] = v;
00610       else if (larger_than_eps) {
00611         if (used_size == data_size)
00612           resize_data(data_size*2+100);
00613         data[used_size] = v;
00614         index[used_size] = i;
00615         used_size++;
00616       }
00617 
00618       //Check if the stored element is smaller than eps. In that case it should be removed.
00619       if (std::abs(v) <= std::abs(eps)) {
00620         remove_small_elements();
00621       }
00622 
00623     }
00624 
00625   template <class T>
00626     void Sparse_Vec<T>::set_new(int i, T v)
00627     {
00628       it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
00629 
00630       //Check that the new element is larger than eps!
00631       if (std::abs(v) > std::abs(eps)) {
00632         if (used_size == data_size)
00633           resize_data(data_size*2+100);
00634         data[used_size] = v;
00635         index[used_size] = i;
00636         used_size++;
00637       }
00638     }
00639 
00640   template <class T>
00641     void Sparse_Vec<T>::add_elem(const int i, const T v)
00642     {
00643       bool found = false;
00644       int p;
00645 
00646       it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
00647 
00648       for (p=0; p<used_size; p++) {
00649         if (index[p] == i) {
00650           found = true;
00651           break;
00652         }
00653       }
00654       if (found)
00655         data[p] += v;
00656       else {
00657         if (used_size == data_size)
00658           resize_data(data_size*2+100);
00659         data[used_size] = v;
00660         index[used_size] = i;
00661         used_size++;
00662       }
00663 
00664       check_small_elems_flag = true;
00665 
00666     }
00667 
00668   template <class T>
00669     void Sparse_Vec<T>::add(const ivec& index_vec, const Vec<T>& v)
00670     {
00671       bool found = false;
00672       int i,p,q;
00673       int nrof_nz=v.size();
00674 
00675       it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
00676 
00677       //Elements are added if they have identical indices
00678       for (q=0; q<nrof_nz; q++){
00679         i=index_vec(q);
00680         for (p=0; p<used_size; p++) {
00681           if (index[p] == i) {
00682             found = true;
00683             break;
00684           }
00685         }
00686         if (found)
00687           data[p] += v(q);
00688         else {
00689           if (used_size == data_size)
00690             resize_data(data_size*2+100);
00691           data[used_size] = v(q);
00692           index[used_size] = i;
00693           used_size++;
00694         }
00695         found = false;
00696       }
00697 
00698       check_small_elems_flag = true;
00699 
00700     }
00701 
00702   template <class T>
00703     void Sparse_Vec<T>::zeros()
00704     {
00705       used_size = 0;
00706       check_small_elems_flag = false;
00707     }
00708 
00709   template <class T>
00710     void Sparse_Vec<T>::zero_elem(const int i)
00711     {
00712       bool found = false;
00713       int p;
00714 
00715       it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
00716 
00717       for (p=0; p<used_size; p++) {
00718         if (index[p] == i) {
00719           found = true;
00720           break;
00721         }
00722       }
00723       if (found) {
00724         data[p] = data[used_size-1];
00725         index[p] = index[used_size-1];
00726         used_size--;
00727       }
00728     }
00729 
00730   template <class T>
00731     void Sparse_Vec<T>::clear()
00732     {
00733       used_size = 0;
00734       check_small_elems_flag = false;
00735     }
00736 
00737   template <class T>
00738     void Sparse_Vec<T>::clear_elem(const int i)
00739     {
00740       bool found = false;
00741       int p;
00742 
00743       it_assert_debug(v_size > i, "The index of the element exceeds the size of the sparse vector");
00744 
00745       for (p=0; p<used_size; p++) {
00746         if (index[p] == i) {
00747           found = true;
00748           break;
00749         }
00750       }
00751       if (found) {
00752         data[p] = data[used_size-1];
00753         index[p] = index[used_size-1];
00754         used_size--;
00755       }
00756     }
00757 
00758   template <class T>
00759     void Sparse_Vec<T>::set(const ivec& index_vec, const Vec<T>& v)
00760     {
00761       it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
00762 
00763       //Clear all old non-zero elements
00764       clear();
00765 
00766       //Add the new non-zero elements
00767       add(index_vec,v);
00768     }
00769 
00770   template <class T>
00771     void Sparse_Vec<T>::set_new(const ivec& index_vec, const Vec<T>& v)
00772     {
00773       int q;
00774       int nrof_nz=v.size();
00775 
00776       it_assert_debug(v_size>max(index_vec),"The indices exceeds the size of the sparse vector");
00777 
00778       //Clear all old non-zero elements
00779       clear();
00780 
00781       for (q=0; q<nrof_nz; q++){
00782         if (std::abs(v[q]) > std::abs(eps)) {
00783           if (used_size == data_size)
00784             resize_data(data_size*2+100);
00785           data[used_size] = v(q);
00786           index[used_size] = index_vec(q);
00787           used_size++;
00788         }
00789       }
00790     }
00791 
00792   template <class T>
00793   ivec Sparse_Vec<T>::get_nz_indices()
00794   {
00795     int n = nnz();
00796     ivec r(n);
00797     for (int i = 0; i < n; i++) {
00798       r(i) = get_nz_index(i);
00799     }
00800     return r;
00801   }
00802 
00803   template <class T>
00804     Sparse_Vec<T> Sparse_Vec<T>::get_subvector(int i1, int i2) const
00805     {
00806       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");
00807 
00808       Sparse_Vec<T> r(i2-i1+1);
00809 
00810       for (int p=0; p<used_size; p++) {
00811         if (index[p] >= i1 && index[p] <= i2) {
00812           if (r.used_size == r.data_size)
00813             r.resize_data(r.data_size*2+100);
00814           r.data[r.used_size] = data[p];
00815           r.index[r.used_size] = index[p]-i1;
00816           r.used_size++;
00817         }
00818       }
00819       r.eps = eps;
00820       r.check_small_elems_flag = check_small_elems_flag;
00821       r.compact();
00822 
00823       return r;
00824     }
00825 
00826   template <class T>
00827     T Sparse_Vec<T>::sqr() const
00828     {
00829       T sum(0);
00830       for (int p=0; p<used_size; p++)
00831         sum += data[p] * data[p];
00832 
00833       return sum;
00834     }
00835 
00836   template <class T>
00837     void Sparse_Vec<T>::operator=(const Sparse_Vec<T> &v)
00838     {
00839       free();
00840       v_size = v.v_size;
00841       used_size = v.used_size;
00842       data_size = v.data_size;
00843       eps = v.eps;
00844       check_small_elems_flag = v.check_small_elems_flag;
00845       alloc();
00846 
00847       for (int i=0; i<used_size; i++) {
00848         data[i] = v.data[i];
00849         index[i] = v.index[i];
00850       }
00851     }
00852 
00853   template <class T>
00854     void Sparse_Vec<T>::operator=(const Vec<T> &v)
00855     {
00856       free();
00857       v_size = v.size();
00858       used_size = 0;
00859       data_size = std::min(v.size(), 10000);
00860       eps = T(0);
00861       check_small_elems_flag = false;
00862       alloc();
00863 
00864       for (int i=0; i<v_size; i++) {
00865         if (v(i) != T(0)) {
00866           if (used_size == data_size)
00867             resize_data(data_size*2);
00868           data[used_size] = v(i);
00869           index[used_size] = i;
00870           used_size++;
00871         }
00872       }
00873       compact();
00874     }
00875 
00876   template <class T>
00877     Sparse_Vec<T> Sparse_Vec<T>::operator-() const
00878     {
00879       Sparse_Vec r(v_size, used_size);
00880 
00881       for (int p=0; p<used_size; p++) {
00882         r.data[p] = -data[p];
00883         r.index[p] = index[p];
00884       }
00885       r.used_size = used_size;
00886 
00887       return r;
00888     }
00889 
00890   template <class T>
00891     bool Sparse_Vec<T>::operator==(const Sparse_Vec<T> &v)
00892     {
00893       int p,q;
00894       bool found=false;
00895 
00896       //Remove small elements before comparing the two sparse_vectors
00897       if (check_small_elems_flag)
00898         remove_small_elements();
00899 
00900       if (v_size!=v.v_size) {
00901         //Return false if vector sizes are unequal
00902         return false;
00903       } else {
00904         for(p=0;p<used_size;p++) {
00905           for(q=0;q<v.used_size;q++) {
00906             if (index[p] == v.index[q]) {
00907               found = true;
00908               break;
00909             }
00910           }
00911           if (found==false)
00912             //Return false if non-zero element not found, or if elements are unequal
00913             return false;
00914           else if (data[p]!=v.data[q])
00915             //Return false if non-zero element not found, or if elements are unequal
00916             return false;
00917           else
00918             //Check next non-zero element
00919             found = false;
00920         }
00921       }
00922 
00923       /*Special handling if sizes do not match.
00924         Required since v may need to do remove_small_elements() for true comparison*/
00925       if (used_size!=v.used_size) {
00926         if (used_size > v.used_size) {
00927           //Return false if number of non-zero elements is less in v
00928           return false;
00929         }
00930         else {
00931           //Ensure that the remaining non-zero elements in v are smaller than v.eps
00932           int nrof_small_elems = 0;
00933           for(q=0;q<v.used_size;q++) {
00934             if (std::abs(v.data[q]) <= std::abs(v.eps))
00935               nrof_small_elems++;
00936           }
00937           if (v.used_size-nrof_small_elems != used_size)
00938             //Return false if the number of "true" non-zero elements are unequal
00939             return false;
00940         }
00941       }
00942 
00943       //All elements checks => return true
00944       return true;
00945     }
00946 
00947   template <class T>
00948     void Sparse_Vec<T>::operator+=(const Sparse_Vec<T> &v)
00949     {
00950       int i,p;
00951       T tmp_data;
00952       int nrof_nz_v=v.used_size;
00953 
00954       it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
00955 
00956       for (p=0; p<nrof_nz_v; p++) {
00957         i = v.index[p];
00958         tmp_data = v.data[p];
00959         //get_nz(p,i,tmp_data);
00960         add_elem(i,tmp_data);
00961       }
00962 
00963       check_small_elems_flag = true;
00964     }
00965 
00966   template <class T>
00967     void Sparse_Vec<T>::operator+=(const Vec<T> &v)
00968     {
00969       int i;
00970 
00971       it_assert_debug(v_size == v.size(), "Attempted addition of unequal sized sparse vectors");
00972 
00973       for (i=0; i<v.size(); i++)
00974         if (v(i)!=T(0))
00975           add_elem(i,v(i));
00976 
00977       check_small_elems_flag = true;
00978     }
00979 
00980 
00981   template <class T>
00982     void Sparse_Vec<T>::operator-=(const Sparse_Vec<T> &v)
00983     {
00984       int i,p;
00985       T tmp_data;
00986       int nrof_nz_v=v.used_size;
00987 
00988       it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
00989 
00990       for (p=0; p<nrof_nz_v; p++) {
00991         i = v.index[p];
00992         tmp_data = v.data[p];
00993         //v.get_nz(p,i,tmp_data);
00994         add_elem(i,-tmp_data);
00995       }
00996 
00997       check_small_elems_flag = true;
00998     }
00999 
01000   template <class T>
01001     void Sparse_Vec<T>::operator-=(const Vec<T> &v)
01002     {
01003       int i;
01004 
01005       it_assert_debug(v_size == v.size(), "Attempted subtraction of unequal sized sparse vectors");
01006 
01007       for (i=0; i<v.size(); i++)
01008         if (v(i)!=T(0))
01009           add_elem(i,-v(i));
01010 
01011       check_small_elems_flag = true;
01012     }
01013 
01014   template <class T>
01015     void Sparse_Vec<T>::operator*=(const T &v)
01016     {
01017       int p;
01018 
01019       for (p=0; p<used_size; p++) {
01020         data[p]*=v;}
01021 
01022       check_small_elems_flag = true;
01023     }
01024 
01025   template <class T>
01026     void Sparse_Vec<T>::operator/=(const T &v)
01027     {
01028       int p;
01029       for (p=0; p<used_size; p++) {
01030         data[p]/=v;}
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 
SourceForge Logo

Generated on Mon Jan 7 22:28:56 2008 for IT++ by Doxygen 1.5.4