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
Generated on Fri Aug 14 15:28:12 2009 for IT++ by Doxygen 1.5.9