IT++ Logo Newcom Logo

smat.h

Go to the documentation of this file.
00001 
00033 #ifndef SMAT_H
00034 #define SMAT_H
00035 
00036 #include <itpp/base/svec.h>
00037 
00038 
00039 namespace itpp {
00040 
00041   // Declaration of class Vec
00042   template <class T> class Vec;
00043   // Declaration of class Mat
00044   template <class T> class Mat;
00045   // Declaration of class Sparse_Vec
00046   template <class T> class Sparse_Vec;
00047   // Declaration of class Sparse_Mat
00048   template <class T> class Sparse_Mat;
00049 
00050   // ------------------------ Sparse_Mat Friends -------------------------------------
00051 
00053   template <class T>
00054     Sparse_Mat<T> operator+(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00055 
00057   template <class T>
00058     Sparse_Mat<T> operator*(const T &c, const Sparse_Mat<T> &m);
00059 
00061   template <class T>
00062     Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00063 
00065   template <class T>
00066     Sparse_Vec<T> operator*(const Sparse_Mat<T> &m, const Sparse_Vec<T> &v);
00067 
00069   template <class T>
00070     Vec<T> operator*(const Sparse_Mat<T> &m, const Vec<T> &v);
00071 
00073   template <class T>
00074     Vec<T> operator*(const Vec<T> &v, const Sparse_Mat<T> &m);
00075 
00077   template <class T>
00078     Mat<T> trans_mult(const Sparse_Mat<T> &m);
00079 
00081   template <class T>
00082     Sparse_Mat<T> trans_mult_s(const Sparse_Mat<T> &m);
00083 
00085   template <class T>
00086     Sparse_Mat<T> trans_mult(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00087 
00089   template <class T>
00090     Vec<T> trans_mult(const Sparse_Mat<T> &m, const Vec<T> &v);
00091 
00093   template <class T>
00094     Sparse_Mat<T> mult_trans(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00095 
00103   template <class T>
00104     class Sparse_Mat {
00105     public:
00106 
00108     Sparse_Mat();
00109 
00120     Sparse_Mat(int rows, int cols, int row_data_init=200);
00121 
00123     Sparse_Mat(const Sparse_Mat<T> &m);
00124 
00126     Sparse_Mat(const Mat<T> &m);
00127 
00133     Sparse_Mat(const Mat<T> &m, T epsilon);
00134 
00136     ~Sparse_Mat();
00137   
00148     void set_size(int rows, int cols, int row_data_init=-1);
00149   
00151     int rows() const { return n_rows; }
00152 
00154     int cols() const { return n_cols; }
00155 
00157     int nnz();
00158 
00160     double density();
00161   
00163     void compact();
00164   
00166     void full(Mat<T> &m) const;
00167 
00169     Mat<T> full() const;
00170   
00172     T operator()(int r, int c) const;
00173 
00175     void set(int r, int c, T v);
00176 
00178     void set_new(int r, int c, T v);
00179 
00181     void add_elem(const int r, const int c, const T v);
00182 
00184     void zeros();
00185 
00187     void zero_elem(const int r, const int c);
00188 
00190     void clear();
00191 
00193     void clear_elem(const int r, const int c);
00194 
00196     void set_submatrix(int r1, int r2, int c1, int c2, const Mat<T> &m);
00197 
00199     void set_submatrix(int r, int c, const Mat<T>& m);
00200 
00202     Sparse_Mat<T> get_submatrix(int r1, int r2, int c1, int c2) const;
00203 
00205     Sparse_Mat<T> get_submatrix_cols(int c1, int c2) const;
00206 
00208     void get_col(int c, Sparse_Vec<T> &v) const;
00209 
00211     Sparse_Vec<T> get_col(int c) const;
00212 
00214     void set_col(int c, const Sparse_Vec<T> &v);
00215   
00217     void transpose(Sparse_Mat<T> &m) const;
00218 
00220     Sparse_Mat<T> transpose() const;
00221 
00223     // Sparse_Mat<T> T() const { return this->transpose(); };
00224   
00226     void operator=(const Sparse_Mat<T> &m);
00227 
00229     void operator=(const Mat<T> &m);
00230   
00232     Sparse_Mat<T> operator-() const;
00233   
00235     bool operator==(const Sparse_Mat<T> &m) const;
00236   
00238     void operator+=(const Sparse_Mat<T> &v);
00239   
00241     void operator+=(const Mat<T> &v);
00242   
00244     void operator-=(const Sparse_Mat<T> &v);
00245   
00247     void operator-=(const Mat<T> &v);
00248   
00250     void operator*=(const T &v);
00251   
00253     void operator/=(const T &v);
00254   
00256     friend Sparse_Mat<T> operator+<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00257 
00259     friend Sparse_Mat<T> operator*<>(const T &c, const Sparse_Mat<T> &m);
00260 
00262     friend Sparse_Mat<T> operator*<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00263 
00265     friend Sparse_Vec<T> operator*<>(const Sparse_Mat<T> &m, const Sparse_Vec<T> &v);
00266 
00268     friend Vec<T> operator*<>(const Sparse_Mat<T> &m, const Vec<T> &v);
00269 
00271     friend Vec<T> operator*<>(const Vec<T> &v, const Sparse_Mat<T> &m);
00272 
00274     friend Mat<T> trans_mult <>(const Sparse_Mat<T> &m);
00275 
00277     friend Sparse_Mat<T> trans_mult_s <>(const Sparse_Mat<T> &m);
00278 
00280     friend Sparse_Mat<T> trans_mult <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00281 
00283     friend Vec<T> trans_mult <>(const Sparse_Mat<T> &m, const Vec<T> &v);
00284 
00286     friend Sparse_Mat<T> mult_trans <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
00287 
00288     private:
00289     void init();
00290     void alloc_empty();
00291     void alloc(int row_data_size=200);
00292     void free();
00293   
00294     int n_rows, n_cols;
00295     Sparse_Vec<T> *col;
00296   };
00297 
00302   typedef Sparse_Mat<int> sparse_imat;
00303 
00308   typedef Sparse_Mat<double> sparse_mat;
00309 
00314   typedef Sparse_Mat<std::complex<double> > sparse_cmat;
00315 
00316   //---------------------- Implementation starts here --------------------------------
00317 
00318   template <class T>
00319     void Sparse_Mat<T>::init()
00320     {
00321       n_rows = 0;
00322       n_cols = 0;
00323       col = 0;
00324     }
00325 
00326   template <class T>
00327     void Sparse_Mat<T>::alloc_empty()
00328     {
00329       if (n_cols == 0)
00330         col = 0;
00331       else
00332         col = new Sparse_Vec<T>[n_cols];
00333     }
00334 
00335   template <class T>
00336     void Sparse_Mat<T>::alloc(int row_data_init)
00337     {
00338       if (n_cols == 0)
00339         col = 0;
00340       else
00341         col = new Sparse_Vec<T>[n_cols];
00342       for (int c=0; c<n_cols; c++)
00343         col[c].set_size(n_rows, row_data_init);
00344     }
00345 
00346   template <class T>
00347     void Sparse_Mat<T>::free()
00348     {
00349       delete [] col;
00350       col = 0;
00351     }
00352 
00353   template <class T>
00354     Sparse_Mat<T>::Sparse_Mat()
00355     {
00356       init();
00357     }
00358 
00359   template <class T>
00360     Sparse_Mat<T>::Sparse_Mat(int rows, int cols, int row_data_init)
00361     {
00362       init();
00363       n_rows = rows;
00364       n_cols = cols;
00365       alloc(row_data_init);
00366     }
00367 
00368   template <class T>
00369     Sparse_Mat<T>::Sparse_Mat(const Sparse_Mat<T> &m)
00370     {
00371       init();
00372       n_rows = m.n_rows;
00373       n_cols = m.n_cols;
00374       alloc_empty();
00375     
00376       for (int c=0; c<n_cols; c++)
00377         col[c] = m.col[c];
00378     }
00379 
00380   template <class T>
00381     Sparse_Mat<T>::Sparse_Mat(const Mat<T> &m)
00382     {
00383       init();
00384       n_rows = m.rows();
00385       n_cols = m.cols();
00386       alloc();
00387 
00388       for (int c=0; c<n_cols; c++) {
00389         for (int r=0; r<n_rows; r++) {
00390           //if (abs(m(r,c)) != T(0))
00391           if (m(r,c) != T(0))
00392             col[c].set_new(r, m(r,c));
00393         }
00394         col[c].compact();
00395       }
00396     }
00397 
00398   template <class T>
00399     Sparse_Mat<T>::Sparse_Mat(const Mat<T> &m, T epsilon)
00400     {
00401       init();
00402       n_rows = m.rows();
00403       n_cols = m.cols();
00404       alloc();
00405 
00406       for (int c=0; c<n_cols; c++) {
00407         for (int r=0; r<n_rows; r++) {
00408                 if (std::abs(m(r,c)) > std::abs(epsilon))
00409             col[c].set_new(r, m(r,c));
00410         }
00411         col[c].compact();
00412       }
00413     }
00414 
00415   template <class T>
00416     Sparse_Mat<T>::~Sparse_Mat()
00417     {
00418       free();
00419     }
00420 
00421   template <class T>
00422     void Sparse_Mat<T>::set_size(int rows, int cols, int row_data_init)
00423     {
00424       n_rows = rows;
00425 
00426       //Allocate new memory for data if the number of columns has changed or if row_data_init != -1
00427       if (cols!=n_cols||row_data_init!=-1) {
00428         n_cols = cols;
00429         free();
00430         alloc(row_data_init);
00431       }
00432     }
00433 
00434   template <class T>
00435     int Sparse_Mat<T>::nnz()
00436     {
00437       int n=0;
00438       for (int c=0; c<n_cols; c++)
00439         n += col[c].nnz();
00440 
00441       return n;
00442     }
00443 
00444   template <class T>
00445     double Sparse_Mat<T>::density()
00446     {
00447       //return static_cast<double>(nnz())/(n_rows*n_cols);
00448       return double(nnz())/(n_rows*n_cols);
00449     }
00450 
00451   template <class T>
00452     void Sparse_Mat<T>::compact()
00453     {
00454       for (int c=0; c<n_cols; c++)
00455         col[c].compact();
00456     }
00457 
00458   template <class T>
00459     void Sparse_Mat<T>::full(Mat<T> &m) const
00460     {
00461       m.set_size(n_rows, n_cols);
00462       m = T(0);
00463       for (int c=0; c<n_cols; c++) {
00464         for (int p=0; p<col[c].nnz(); p++)
00465           m(col[c].get_nz_index(p),c) = col[c].get_nz_data(p);
00466       }
00467     }
00468 
00469   template <class T>
00470     Mat<T> Sparse_Mat<T>::full() const
00471     {
00472       Mat<T> r(n_rows, n_cols);
00473       full(r);
00474       return r;
00475     }
00476 
00477   template <class T>
00478     T Sparse_Mat<T>::operator()(int r, int c) const
00479     {
00480       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00481       return col[c](r);
00482     }
00483 
00484   template <class T>
00485     void Sparse_Mat<T>::set(int r, int c, T v)
00486     {
00487       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00488       col[c].set(r, v);
00489     }
00490 
00491   template <class T>
00492     void Sparse_Mat<T>::set_new(int r, int c, T v)
00493     {
00494       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00495       col[c].set_new(r, v);
00496     }
00497 
00498   template <class T>
00499     void Sparse_Mat<T>::add_elem(int r, int c, T v)
00500     {
00501       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00502       col[c].add_elem(r, v);
00503     }
00504 
00505   template <class T>
00506     void Sparse_Mat<T>::zeros()
00507     {
00508       for (int c=0; c<n_cols; c++)
00509         col[c].zeros();      
00510     }
00511 
00512   template <class T>
00513     void Sparse_Mat<T>::zero_elem(const int r, const int c)
00514     {
00515       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00516       col[c].zero_elem(r);
00517     }
00518 
00519   template <class T>
00520     void Sparse_Mat<T>::clear()
00521     {
00522       for (int c=0; c<n_cols; c++)
00523         col[c].clear();      
00524     }
00525 
00526   template <class T>
00527     void Sparse_Mat<T>::clear_elem(const int r, const int c)
00528     {
00529       it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
00530       col[c].clear_elem(r);
00531     }
00532 
00533   template <class T>
00534     void Sparse_Mat<T>::set_submatrix(int r1, int r2, int c1, int c2, const Mat<T>& m) 
00535   {
00536     if (r1 == -1) r1 = n_rows-1;
00537     if (r2 == -1) r2 = n_rows-1;
00538     if (c1 == -1) c1 = n_cols-1;
00539     if (c2 == -1) c2 = n_cols-1;
00540 
00541     it_assert1(r1>=0 && r2>=0 && r1<n_rows && r2<n_rows &&
00542                c1>=0 && c2>=0 && c1<n_cols && c2<n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
00543 
00544     it_assert1(r2>=r1 && c2>=c1, "Sparse_Mat<Num_T>::set_submatrix: r2<r1 or c2<c1");
00545     it_assert1(m.rows() == r2-r1+1 && m.cols() == c2-c1+1, "Mat<Num_T>::set_submatrix(): sizes don't match");
00546 
00547     for (int i=0 ; i<m.rows() ; i++) {
00548       for (int j=0 ; j<m.cols() ; j++) {
00549         set(r1+i, c1+j, m(i,j));
00550       }
00551     }
00552   }
00553 
00554   template <class T>
00555     void Sparse_Mat<T>::set_submatrix(int r, int c, const Mat<T>& m) 
00556   {
00557     it_assert1(r>=0 && r+m.rows()<=n_rows &&
00558                c>=0 && c+m.cols()<=n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
00559 
00560     for (int i=0 ; i<m.rows() ; i++) {
00561       for (int j=0 ; j<m.cols() ; j++) {
00562         set(r+i, c+j, m(i,j));
00563       }
00564     }
00565   }
00566 
00567   template <class T>
00568     Sparse_Mat<T> Sparse_Mat<T>::get_submatrix(int r1, int r2, int c1, int c2) const
00569     {
00570       it_assert0(r1<=r2 && r1>=0 && r1<n_rows && c1<=c2 && c1>=0 && c1<n_cols,
00571                 "Sparse_Mat<T>::get_submatrix(): illegal input variables");
00572     
00573       Sparse_Mat<T> r(r2-r1+1, c2-c1+1);
00574 
00575       for (int c=c1; c<=c2; c++)
00576         r.col[c-c1] = col[c].get_subvector(r1, r2);
00577       r.compact();
00578 
00579       return r;
00580     }
00581 
00582   template <class T>
00583     Sparse_Mat<T> Sparse_Mat<T>::get_submatrix_cols(int c1, int c2) const
00584     {
00585       it_assert0(c1<=c2 && c1>=0 && c1<n_cols, "Sparse_Mat<T>::get_submatrix_cols()");
00586       Sparse_Mat<T> r(n_rows, c2-c1+1, 0);
00587 
00588       for (int c=c1; c<=c2; c++)
00589         r.col[c-c1] = col[c];
00590       r.compact();
00591 
00592       return r;
00593     }
00594 
00595   template <class T>
00596     void Sparse_Mat<T>::get_col(int c, Sparse_Vec<T> &v) const
00597     {
00598       it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::get_col()");
00599       v = col[c];
00600     }
00601 
00602   template <class T>
00603     Sparse_Vec<T> Sparse_Mat<T>::get_col(int c) const
00604     {
00605       it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::get_col()");
00606       return col[c];
00607     }
00608 
00609   template <class T>
00610     void Sparse_Mat<T>::set_col(int c, const Sparse_Vec<T> &v)
00611     {
00612       it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::set_col()");
00613       col[c] = v;
00614     }
00615 
00616   template <class T>
00617     void Sparse_Mat<T>::transpose(Sparse_Mat<T> &m) const
00618   {
00619     m.set_size(n_cols, n_rows);
00620     for (int c=0; c<n_cols; c++) {
00621       for (int p=0; p<col[c].nnz(); p++)
00622         m.col[col[c].get_nz_index(p)].set_new(c, col[c].get_nz_data(p));
00623     }
00624   }
00625 
00626   template <class T>
00627     Sparse_Mat<T> Sparse_Mat<T>::transpose() const
00628   {
00629     Sparse_Mat<T> m;
00630     transpose(m);
00631     return m;
00632   }
00633 
00634   template <class T>
00635     void Sparse_Mat<T>::operator=(const Sparse_Mat<T> &m)
00636     {
00637       free();
00638       n_rows = m.n_rows;
00639       n_cols = m.n_cols;
00640       alloc_empty();
00641     
00642       for (int c=0; c<n_cols; c++)
00643         col[c] = m.col[c];
00644     }
00645 
00646   template <class T>
00647     void Sparse_Mat<T>::operator=(const Mat<T> &m)
00648     {
00649       free();
00650       n_rows = m.rows();
00651       n_cols = m.cols();
00652       alloc();
00653 
00654       for (int c=0; c<n_cols; c++) {
00655         for (int r=0; r<n_rows; r++) {
00656           if (m(r,c) != T(0))
00657             col[c].set_new(r, m(r,c));
00658         }
00659         col[c].compact();
00660       }
00661     }
00662 
00663   template <class T>
00664     Sparse_Mat<T> Sparse_Mat<T>::operator-() const
00665     {
00666       Sparse_Mat r(n_rows, n_cols, 0);
00667 
00668       for (int c=0; c<n_cols; c++) {
00669         r.col[c].resize_data(col[c].nnz());
00670         for (int p=0; p<col[c].nnz(); p++)
00671           r.col[c].set_new(col[c].get_nz_index(p), -col[c].get_nz_data(p));
00672       }
00673 
00674       return r;
00675     }
00676 
00677   template <class T>
00678     bool Sparse_Mat<T>::operator==(const Sparse_Mat<T> &m) const
00679     {
00680       if (n_rows!=m.n_rows || n_cols!=m.n_cols)
00681         return false;
00682       for (int c=0; c<n_cols; c++) {
00683         if (!(col[c] == m.col[c]))
00684           return false;
00685       }
00686       // If they passed all tests, they must be equal
00687       return true;
00688     }
00689 
00690   template <class T>
00691     void Sparse_Mat<T>::operator+=(const Sparse_Mat<T> &m)
00692     {
00693       it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Addition of unequal sized matrices is not allowed");
00694 
00695       Sparse_Vec<T> v;
00696       for (int c=0; c<n_cols; c++) {
00697         m.get_col(c,v);
00698         col[c]+=v;
00699       }
00700     }
00701 
00702   template <class T>
00703     void Sparse_Mat<T>::operator+=(const Mat<T> &m)
00704     {
00705       it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Addition of unequal sized matrices is not allowed");
00706 
00707       for (int c=0; c<n_cols; c++)
00708         col[c]+=(m.get_col(c));
00709     }
00710 
00711   template <class T>
00712     void Sparse_Mat<T>::operator-=(const Sparse_Mat<T> &m)
00713     {
00714       it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Subtraction of unequal sized matrices is not allowed");
00715 
00716       Sparse_Vec<T> v;
00717       for (int c=0; c<n_cols; c++) {
00718         m.get_col(c,v);
00719         col[c]-=v;
00720       }
00721     }
00722 
00723   template <class T>
00724     void Sparse_Mat<T>::operator-=(const Mat<T> &m)
00725     {
00726       it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Subtraction of unequal sized matrices is not allowed");
00727 
00728       for (int c=0; c<n_cols; c++)
00729         col[c]-=(m.get_col(c));
00730     }
00731 
00732   template <class T>
00733     void Sparse_Mat<T>::operator*=(const T &m)
00734     {
00735       for (int c=0; c<n_cols; c++)
00736         col[c]*=m;
00737     }
00738 
00739   template <class T>
00740     void Sparse_Mat<T>::operator/=(const T &m)
00741     {
00742       for (int c=0; c<n_cols; c++)
00743         col[c]/=m;
00744     }
00745 
00746   template <class T>
00747     Sparse_Mat<T> operator+(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
00748     {
00749       it_assert0(m1.n_cols == m2.n_cols && m1.n_rows == m2.n_rows , "Sparse_Mat<T> + Sparse_Mat<T>");
00750     
00751       Sparse_Mat<T> m(m1.n_rows, m1.n_cols, 0);
00752 
00753       for (int c=0; c<m.n_cols; c++)
00754         m.col[c] = m1.col[c] + m2.col[c];
00755 
00756       return m;
00757     }
00758 
00759   // This function added by EGL, May'05
00760   template <class T>
00761     Sparse_Mat<T> operator*(const T &c, const Sparse_Mat<T> &m)
00762     {
00763       int i,j;
00764       Sparse_Mat<T> ret(m.n_rows,m.n_cols);
00765       for (j=0; j<m.n_cols; j++) {
00766         for (i=0; i<m.col[j].nnz(); i++) {
00767           T x = c*m.col[j].get_nz_data(i);
00768           int k = m.col[j].get_nz_index(i);
00769           ret.set_new(k,j,x);
00770         }
00771       }
00772       return ret;
00773     }
00774 
00775   // modified implementation by EGL, May'05
00776   template <class T>
00777     Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
00778     {
00779       it_assert0(m1.n_cols == m2.n_rows, "Sparse_Mat<T> * Sparse_Mat<T>");
00780       
00781       Sparse_Mat<T> ret(m1.n_rows, m2.n_cols);
00782       
00783       for (int c=0; c<m2.n_cols; c++) {
00784         for (int p2=0; p2<m2.col[c].nnz(); p2++) {
00785           Sparse_Vec<T> &mcol = m1.col[m2.col[c].get_nz_index(p2)];
00786           for (int p1=0; p1<mcol.nnz(); p1++) {
00787             int r = mcol.get_nz_index(p1);
00788             T inc = mcol.get_nz_data(p1) * m2.col[c].get_nz_data(p2);
00789             ret.col[c].add_elem(r,inc);
00790           }
00791         }
00792       }
00793       ret.compact();
00794       return ret;
00795     }
00796 
00797 
00798   // This is apparently buggy.
00799 /*   template <class T> */
00800 /*     Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2) */
00801 /*     { */
00802 /*       it_assert0(m1.n_cols == m2.n_rows, "Sparse_Mat<T> * Sparse_Mat<T>"); */
00803     
00804 /*       Sparse_Mat<T> ret(m1.n_rows, m2.n_cols); */
00805 /*       ivec occupied_by(ret.n_rows), pos(ret.n_rows); */
00806 /*       for (int rp=0; rp<m1.n_rows; rp++) */
00807 /*      occupied_by[rp] = -1; */
00808 /*       for (int c=0; c<ret.n_cols; c++) { */
00809 /*      Sparse_Vec<T> &m2col = m2.col[c]; */
00810 /*      for (int p2=0; p2<m2col.nnz(); p2++) { */
00811 /*        Sparse_Vec<T> &m1col = m1.col[m2col.get_nz_index(p2)]; */
00812 /*        for (int p1=0; p1<m1col.nnz(); p1++) { */
00813 /*          int r = m1col.get_nz_index(p1); */
00814 /*          T inc = m1col.get_nz_data(p1) * m2col.get_nz_data(p2); */
00815 /*          if (occupied_by[r] == c) { */
00816 /*            int index=ret.col[c].get_nz_index(pos[r]); */
00817 /*            ret.col[c].add_elem(index,inc); */
00818 /*          } */
00819 /*          else { */
00820 /*            occupied_by[r] = c; */
00821 /*            pos[r] = ret.col[c].nnz(); */
00822 /*            ret.col[c].set_new(r, inc); */
00823 /*          } */
00824 /*        } */
00825 /*      } */
00826 /*       } */
00827 /*       ret.compact(); */
00828 
00829 /*       return ret; */
00830 /*     } */
00831 
00832 
00833   // This function added by EGL, May'05
00834   template <class T>
00835     Sparse_Vec<T> operator*(const Sparse_Mat<T> &m, const Sparse_Vec<T> &v)
00836     {
00837       it_assert0(m.n_cols == v.size(), "Sparse_Mat<T> * Sparse_Vec<T>");
00838 
00839       Sparse_Vec<T> ret(m.n_rows);
00840 
00841       /* The two lines below added because the input parameter "v" is
00842          declared const, but the some functions (e.g., nnz()) change
00843          the vector... Is there a better workaround? */
00844       Sparse_Vec<T> vv;
00845       vv = v;
00846 
00847       for (int p2=0; p2<vv.nnz(); p2++) {
00848         Sparse_Vec<T> &mcol = m.col[vv.get_nz_index(p2)];
00849         for (int p1=0; p1<mcol.nnz(); p1++) {
00850           int r = mcol.get_nz_index(p1);
00851           T inc = mcol.get_nz_data(p1) * vv.get_nz_data(p2);
00852           ret.add_elem(r,inc);
00853         }
00854       }
00855       ret.compact();
00856       return ret;
00857     }
00858 
00859 
00860   template <class T>
00861     Vec<T> operator*(const Sparse_Mat<T> &m, const Vec<T> &v)
00862     {
00863       it_assert0(m.n_cols == v.size(), "Sparse_Mat<T> * Vec<T>");
00864     
00865       Vec<T> r(m.n_rows);
00866       r.clear();
00867 
00868       for (int c=0; c<m.n_cols; c++) {
00869         for (int p=0; p<m.col[c].nnz(); p++)
00870           r(m.col[c].get_nz_index(p)) += m.col[c].get_nz_data(p) * v(c);
00871       }
00872       
00873       return r;
00874     }
00875 
00876   template <class T>
00877     Vec<T> operator*(const Vec<T> &v, const Sparse_Mat<T> &m)
00878     {
00879       it_assert0(v.size() == m.n_rows, "Vec<T> * Sparse_Mat<T>");
00880     
00881       Vec<T> r(m.n_cols);
00882       r.clear();
00883 
00884       for (int c=0; c<m.n_cols; c++)
00885         r[c] = v * m.col[c];
00886 
00887       return r;
00888     }
00889 
00890   template <class T>
00891     Mat<T> trans_mult(const Sparse_Mat<T> &m)
00892     {
00893       Mat<T> ret(m.n_cols, m.n_cols);
00894       Vec<T> col;
00895       for (int c=0; c<ret.cols(); c++) {
00896         m.col[c].full(col);
00897         for (int r=0; r<c; r++) {
00898           T tmp = m.col[r] * col;
00899           ret(r,c) = tmp;
00900           ret(c,r) = tmp;
00901         }
00902         ret(c,c) = m.col[c].sqr();
00903       }
00904 
00905       return ret;
00906     }
00907 
00908   template <class T>
00909     Sparse_Mat<T> trans_mult_s(const Sparse_Mat<T> &m)
00910     {
00911       Sparse_Mat<T> ret(m.n_cols, m.n_cols);
00912       Vec<T> col;
00913       T tmp;
00914       for (int c=0; c<ret.n_cols; c++) {
00915         m.col[c].full(col);
00916         for (int r=0; r<c; r++) {
00917           tmp = m.col[r] * col;
00918           if (tmp != T(0)) {
00919             ret.col[c].set_new(r, tmp);
00920             ret.col[r].set_new(c, tmp);
00921           }
00922         }
00923         tmp = m.col[c].sqr();
00924         if (tmp != T(0))
00925           ret.col[c].set_new(c, tmp);
00926       }
00927 
00928       return ret;
00929     }
00930 
00931   template <class T>
00932     Sparse_Mat<T> trans_mult(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
00933     {
00934       it_assert0(m1.n_rows == m2.n_rows, "trans_mult()");
00935     
00936       Sparse_Mat<T> ret(m1.n_cols, m2.n_cols);
00937       Vec<T> col;
00938       for (int c=0; c<ret.n_cols; c++) {
00939         m2.col[c].full(col);
00940         for (int r=0; r<ret.n_rows; r++)
00941           ret.col[c].set_new(r, m1.col[r] * col);
00942       }
00943 
00944       return ret;
00945     }
00946 
00947   template <class T>
00948     Vec<T> trans_mult(const Sparse_Mat<T> &m, const Vec<T> &v)
00949     {
00950       Vec<T> r(m.n_cols);
00951       for (int c=0; c<m.n_cols; c++)
00952         r(c) = m.col[c] * v;
00953 
00954       return r;
00955     }
00956 
00957   template <class T>
00958     Sparse_Mat<T> mult_trans(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
00959     {
00960       return trans_mult(m1.transpose(),m2.transpose());
00961     }
00962 
00963   template <class T>
00964     inline Sparse_Mat<T> sparse(const Mat<T> &m, T epsilon)
00965     {
00966       Sparse_Mat<T> s(m, epsilon);
00967       return s;
00968     }
00969 
00970   template <class T>
00971     inline Mat<T> full(const Sparse_Mat<T> &s)
00972     {
00973       Mat<T> m;
00974       s.full(m);
00975       return m;
00976     }
00977 
00978   template <class T>
00979     inline Sparse_Mat<T> transpose(const Sparse_Mat<T> &s)
00980     {
00981       Sparse_Mat<T> m;
00982       s.transpose(m);
00983       return m;
00984     }
00985 
00986   /*--------------------------------------------------------------
00987    * Explicit initializations
00988    *--------------------------------------------------------------*/
00989 
00990 #ifndef _MSC_VER
00991 
00993   extern template class Sparse_Mat<int>;
00995   extern template class Sparse_Mat<double>;
00997   extern template class Sparse_Mat<std::complex<double> >;
00998 
00999 
01001   extern template sparse_imat operator+(const sparse_imat &, const sparse_imat &);
01003   extern template sparse_mat operator+(const sparse_mat &, const sparse_mat &);
01005   extern template sparse_cmat operator+(const sparse_cmat &, const sparse_cmat &);
01006 
01008   extern template sparse_imat operator*(const sparse_imat &, const sparse_imat &);
01010   extern template sparse_mat operator*(const sparse_mat &, const sparse_mat &);
01012   extern template sparse_cmat operator*(const sparse_cmat &, const sparse_cmat &);
01013 
01015   extern template ivec operator*(const ivec &, const sparse_imat &);
01017   extern template vec operator*(const vec &, const sparse_mat &);
01019   extern template cvec operator*(const cvec &, const sparse_cmat &);
01020 
01022   extern template ivec operator*(const sparse_imat &, const ivec &);
01024   extern template vec operator*(const sparse_mat &, const vec &);
01026   extern template cvec operator*(const sparse_cmat &, const cvec &);
01027 
01029   extern template imat trans_mult(const sparse_imat &);
01031   extern template mat trans_mult(const sparse_mat &);
01033   extern template cmat trans_mult(const sparse_cmat &);
01034 
01036   extern template sparse_imat trans_mult_s(const sparse_imat &);
01038   extern template sparse_mat trans_mult_s(const sparse_mat &);
01040   extern template sparse_cmat trans_mult_s(const sparse_cmat &);
01041 
01043   extern template sparse_imat trans_mult(const sparse_imat &, const sparse_imat &);
01045   extern template sparse_mat trans_mult(const sparse_mat &, const sparse_mat &);
01047   extern template sparse_cmat trans_mult(const sparse_cmat &, const sparse_cmat &);
01048 
01050   extern template ivec trans_mult(const sparse_imat &, const ivec &);
01052   extern template vec trans_mult(const sparse_mat &, const vec &);
01054   extern template cvec trans_mult(const sparse_cmat &, const cvec &);
01055 
01057   extern template sparse_imat mult_trans(const sparse_imat &, const sparse_imat &);
01059   extern template sparse_mat mult_trans(const sparse_mat &, const sparse_mat &);
01061   extern template sparse_cmat mult_trans(const sparse_cmat &, const sparse_cmat &);
01062 
01063 #endif
01064 
01065 } // namespace itpp
01066 
01067 #endif // #ifndef SMAT_H
01068 
SourceForge Logo

Generated on Wed Mar 21 12:21:37 2007 for IT++ by Doxygen 1.4.7