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