IT++ Logo

smat.h

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

Generated on Sun Dec 9 17:30:25 2007 for IT++ by Doxygen 1.5.4