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