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

Generated on Fri Aug 14 15:28:12 2009 for IT++ by Doxygen 1.5.9