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

Generated on Sat Jul 9 2011 15:21:30 for IT++ by Doxygen 1.7.4