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
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4