00001 00029 #ifndef MAT_H 00030 #define MAT_H 00031 00032 #include <itpp/base/itassert.h> 00033 #include <itpp/base/math/misc.h> 00034 #include <itpp/base/factory.h> 00035 00036 00037 namespace itpp 00038 { 00039 00040 // Declaration of Vec 00041 template<class Num_T> class Vec; 00042 // Declaration of Mat 00043 template<class Num_T> class Mat; 00044 // Declaration of bin 00045 class bin; 00046 00048 template<class Num_T> 00049 Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00051 template<class Num_T> 00052 Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00053 00055 template<class Num_T> 00056 Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00058 template<class Num_T> 00059 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t); 00061 template<class Num_T> 00062 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m); 00063 00065 template<class Num_T> 00066 Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00068 template<class Num_T> 00069 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t); 00071 template<class Num_T> 00072 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m); 00074 template<class Num_T> 00075 Mat<Num_T> operator-(const Mat<Num_T> &m); 00076 00078 template<class Num_T> 00079 Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00081 template<class Num_T> 00082 Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v); 00084 template<class Num_T> 00085 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t); 00087 template<class Num_T> 00088 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m); 00089 00091 template<class Num_T> 00092 Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00094 template<class Num_T> 00095 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00096 Mat<Num_T> &out); 00098 template<class Num_T> 00099 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00100 const Mat<Num_T> &m3, Mat<Num_T> &out); 00102 template<class Num_T> 00103 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00104 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 00105 Mat<Num_T> &out); 00107 template<class Num_T> 00108 void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2); 00110 template<class Num_T> 00111 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00112 00114 template<class Num_T> 00115 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t); 00117 template<class Num_T> 00118 Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m); 00119 00121 template<class Num_T> 00122 Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00124 template<class Num_T> 00125 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00126 Mat<Num_T> &out); 00128 template<class Num_T> 00129 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00130 00131 // ------------------------------------------------------------------------------------- 00132 // Declaration of Mat 00133 // ------------------------------------------------------------------------------------- 00134 00200 template<class Num_T> 00201 class Mat 00202 { 00203 public: 00205 typedef Num_T value_type; 00206 00208 explicit Mat(const Factory &f = DEFAULT_FACTORY); 00210 Mat(int rows, int cols, const Factory &f = DEFAULT_FACTORY); 00212 Mat(const Mat<Num_T> &m); 00214 Mat(const Mat<Num_T> &m, const Factory &f); 00216 Mat(const Vec<Num_T> &v, const Factory &f = DEFAULT_FACTORY); 00218 Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY); 00220 Mat(const char *str, const Factory &f = DEFAULT_FACTORY); 00228 Mat(const Num_T *c_array, int rows, int cols, bool row_major = true, 00229 const Factory &f = DEFAULT_FACTORY); 00230 00232 ~Mat(); 00233 00235 int cols() const { return no_cols; } 00237 int rows() const { return no_rows; } 00239 int size() const { return datasize; } 00241 void set_size(int rows, int cols, bool copy = false); 00243 void zeros(); 00245 void clear() { zeros(); } 00247 void ones(); 00249 void set(const std::string &str); 00251 void set(const char *str); 00252 00254 const Num_T &operator()(int r, int c) const; 00256 Num_T &operator()(int r, int c); 00258 const Num_T &operator()(int i) const; 00260 Num_T &operator()(int i); 00262 const Num_T &get(int r, int c) const; 00264 const Num_T &get(int i) const; 00266 void set(int r, int c, Num_T t); 00267 00273 Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const; 00279 Mat<Num_T> get(int r1, int r2, int c1, int c2) const; 00280 00282 Vec<Num_T> get_row(int r) const; 00284 Mat<Num_T> get_rows(int r1, int r2) const; 00286 Mat<Num_T> get_rows(const Vec<int> &indexlist) const; 00288 Vec<Num_T> get_col(int c) const; 00290 Mat<Num_T> get_cols(int c1, int c2) const; 00292 Mat<Num_T> get_cols(const Vec<int> &indexlist) const; 00294 void set_row(int r, const Vec<Num_T> &v); 00296 void set_col(int c, const Vec<Num_T> &v); 00298 void set_rows(int r, const Mat<Num_T> &m); 00300 void set_cols(int c, const Mat<Num_T> &m); 00302 void copy_row(int to, int from); 00304 void copy_col(int to, int from); 00306 void swap_rows(int r1, int r2); 00308 void swap_cols(int c1, int c2); 00309 00311 void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m); 00313 void set_submatrix(int r, int c, const Mat<Num_T> &m); 00315 void set_submatrix(int r1, int r2, int c1, int c2, Num_T t); 00316 00318 void del_row(int r); 00320 void del_rows(int r1, int r2); 00322 void del_col(int c); 00324 void del_cols(int c1, int c2); 00326 void ins_row(int r, const Vec<Num_T> &v); 00328 void ins_col(int c, const Vec<Num_T> &v); 00330 void append_row(const Vec<Num_T> &v); 00332 void append_col(const Vec<Num_T> &v); 00333 00335 Mat<Num_T> transpose() const; 00337 Mat<Num_T> T() const { return this->transpose(); } 00339 Mat<Num_T> hermitian_transpose() const; 00341 Mat<Num_T> H() const { return this->hermitian_transpose(); } 00342 00344 friend Mat<Num_T> concat_horizontal<>(const Mat<Num_T> &m1, 00345 const Mat<Num_T> &m2); 00347 friend Mat<Num_T> concat_vertical<>(const Mat<Num_T> &m1, 00348 const Mat<Num_T> &m2); 00349 00351 Mat<Num_T>& operator=(Num_T t); 00353 Mat<Num_T>& operator=(const Mat<Num_T> &m); 00355 Mat<Num_T>& operator=(const Vec<Num_T> &v); 00357 Mat<Num_T>& operator=(const std::string &str); 00359 Mat<Num_T>& operator=(const char *str); 00360 00362 Mat<Num_T>& operator+=(const Mat<Num_T> &m); 00364 Mat<Num_T>& operator+=(Num_T t); 00366 friend Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00368 friend Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t); 00370 friend Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m); 00371 00373 Mat<Num_T>& operator-=(const Mat<Num_T> &m); 00375 Mat<Num_T>& operator-=(Num_T t); 00377 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00379 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t); 00381 friend Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m); 00383 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m); 00384 00386 Mat<Num_T>& operator*=(const Mat<Num_T> &m); 00388 Mat<Num_T>& operator*=(Num_T t); 00390 friend Mat<Num_T> operator*<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00392 friend Vec<Num_T> operator*<>(const Mat<Num_T> &m, const Vec<Num_T> &v); 00394 friend Mat<Num_T> operator*<>(const Mat<Num_T> &m, Num_T t); 00396 friend Mat<Num_T> operator*<>(Num_T t, const Mat<Num_T> &m); 00397 00399 friend Mat<Num_T> elem_mult<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00401 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00402 Mat<Num_T> &out); 00404 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00405 const Mat<Num_T> &m3, Mat<Num_T> &out); 00407 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00408 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 00409 Mat<Num_T> &out); 00411 friend void elem_mult_inplace<>(const Mat<Num_T> &m1, Mat<Num_T> &m2); 00413 friend Num_T elem_mult_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00414 00416 Mat<Num_T>& operator/=(Num_T t); 00418 Mat<Num_T>& operator/=(const Mat<Num_T> &m); 00419 00421 friend Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t); 00423 friend Mat<Num_T> operator/<>(Num_T t, const Mat<Num_T> &m); 00424 00426 friend Mat<Num_T> elem_div<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00428 friend void elem_div_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00429 Mat<Num_T> &out); 00431 friend Num_T elem_div_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00432 00434 bool operator==(const Mat<Num_T> &m) const; 00436 bool operator!=(const Mat<Num_T> &m) const; 00437 00439 Num_T &_elem(int r, int c) { return data[r+c*no_rows]; } 00441 const Num_T &_elem(int r, int c) const { return data[r+c*no_rows]; } 00443 Num_T &_elem(int i) { return data[i]; } 00445 const Num_T &_elem(int i) const { return data[i]; } 00446 00448 Num_T *_data() { return data; } 00450 const Num_T *_data() const { return data; } 00452 int _datasize() const { return datasize; } 00453 00454 protected: 00456 void alloc(int rows, int cols); 00458 void free(); 00459 00462 int datasize, no_rows, no_cols; 00464 00465 Num_T *data; 00467 const Factory &factory; 00468 00469 private: 00471 bool in_range(int r, int c) const { 00472 return ((r >= 0) && (r < no_rows) && (c >= 0) && (c < no_cols)); 00473 } 00475 bool row_in_range(int r) const { return ((r >= 0) && (r < no_rows)); } 00477 bool col_in_range(int c) const { return ((c >= 0) && (c < no_cols)); } 00479 bool in_range(int i) const { return ((i >= 0) && (i < datasize)); } 00480 }; 00481 00482 // ------------------------------------------------------------------------------------- 00483 // Type definitions of mat, cmat, imat, smat, and bmat 00484 // ------------------------------------------------------------------------------------- 00485 00490 typedef Mat<double> mat; 00491 00496 typedef Mat<std::complex<double> > cmat; 00497 00502 typedef Mat<int> imat; 00503 00508 typedef Mat<short int> smat; 00509 00516 typedef Mat<bin> bmat; 00517 00518 } //namespace itpp 00519 00520 00521 #include <itpp/base/vec.h> 00522 00523 namespace itpp 00524 { 00525 00526 // ---------------------------------------------------------------------- 00527 // Declaration of input and output streams for Mat 00528 // ---------------------------------------------------------------------- 00529 00534 template <class Num_T> 00535 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m); 00536 00548 template <class Num_T> 00549 std::istream &operator>>(std::istream &is, Mat<Num_T> &m); 00550 00551 // ---------------------------------------------------------------------- 00552 // Implementation of templated Mat members and friends 00553 // ---------------------------------------------------------------------- 00554 00555 template<class Num_T> inline 00556 void Mat<Num_T>::alloc(int rows, int cols) 00557 { 00558 if ((rows > 0) && (cols > 0)) { 00559 datasize = rows * cols; 00560 no_rows = rows; 00561 no_cols = cols; 00562 create_elements(data, datasize, factory); 00563 } 00564 else { 00565 data = 0; 00566 datasize = 0; 00567 no_rows = 0; 00568 no_cols = 0; 00569 } 00570 } 00571 00572 template<class Num_T> inline 00573 void Mat<Num_T>::free() 00574 { 00575 destroy_elements(data, datasize); 00576 datasize = 0; 00577 no_rows = 0; 00578 no_cols = 0; 00579 } 00580 00581 00582 template<class Num_T> inline 00583 Mat<Num_T>::Mat(const Factory &f) : 00584 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) {} 00585 00586 template<class Num_T> inline 00587 Mat<Num_T>::Mat(int rows, int cols, const Factory &f) : 00588 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00589 { 00590 it_assert_debug((rows >= 0) && (cols >= 0), "Mat<>::Mat(): Wrong size"); 00591 alloc(rows, cols); 00592 } 00593 00594 template<class Num_T> inline 00595 Mat<Num_T>::Mat(const Mat<Num_T> &m) : 00596 datasize(0), no_rows(0), no_cols(0), data(0), factory(m.factory) 00597 { 00598 alloc(m.no_rows, m.no_cols); 00599 copy_vector(m.datasize, m.data, data); 00600 } 00601 00602 template<class Num_T> inline 00603 Mat<Num_T>::Mat(const Mat<Num_T> &m, const Factory &f) : 00604 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00605 { 00606 alloc(m.no_rows, m.no_cols); 00607 copy_vector(m.datasize, m.data, data); 00608 } 00609 00610 template<class Num_T> inline 00611 Mat<Num_T>::Mat(const Vec<Num_T> &v, const Factory &f) : 00612 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00613 { 00614 int size = v.size(); 00615 alloc(size, 1); 00616 copy_vector(size, v._data(), data); 00617 } 00618 00619 template<class Num_T> inline 00620 Mat<Num_T>::Mat(const std::string &str, const Factory &f) : 00621 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00622 { 00623 set(str); 00624 } 00625 00626 template<class Num_T> inline 00627 Mat<Num_T>::Mat(const char *str, const Factory &f) : 00628 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00629 { 00630 set(std::string(str)); 00631 } 00632 00633 template<class Num_T> 00634 Mat<Num_T>::Mat(const Num_T *c_array, int rows, int cols, bool row_major, 00635 const Factory &f): 00636 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00637 { 00638 alloc(rows, cols); 00639 if (!row_major) 00640 copy_vector(datasize, c_array, data); 00641 else 00642 for (int i = 0; i < rows; i++) 00643 for (int j = 0; j < cols; j++) 00644 data[i+j*no_rows] = c_array[i*no_cols+j]; 00645 } 00646 00647 template<class Num_T> inline 00648 Mat<Num_T>::~Mat() 00649 { 00650 free(); 00651 } 00652 00653 00654 template<class Num_T> 00655 void Mat<Num_T>::set_size(int rows, int cols, bool copy) 00656 { 00657 it_assert_debug((rows >= 0) && (cols >= 0), 00658 "Mat<>::set_size(): Wrong size"); 00659 // check if we have to resize the current matrix 00660 if ((no_rows == rows) && (no_cols == cols)) 00661 return; 00662 // check if one of dimensions is zero 00663 if ((rows == 0) || (cols == 0)) { 00664 free(); 00665 return; 00666 } 00667 // conditionally copy previous matrix content 00668 if (copy) { 00669 // create a temporary pointer to the allocated data 00670 Num_T* tmp = data; 00671 // store the current number of elements and number of rows 00672 int old_datasize = datasize; 00673 int old_rows = no_rows; 00674 // check the boundaries of the copied data 00675 int min_r = (no_rows < rows) ? no_rows : rows; 00676 int min_c = (no_cols < cols) ? no_cols : cols; 00677 // allocate new memory 00678 alloc(rows, cols); 00679 // copy the previous data into the allocated memory 00680 for (int i = 0; i < min_c; ++i) { 00681 copy_vector(min_r, &tmp[i*old_rows], &data[i*no_rows]); 00682 } 00683 // fill-in the rest of matrix with zeros 00684 for (int i = min_r; i < rows; ++i) 00685 for (int j = 0; j < cols; ++j) 00686 data[i+j*rows] = Num_T(0); 00687 for (int j = min_c; j < cols; ++j) 00688 for (int i = 0; i < min_r; ++i) 00689 data[i+j*rows] = Num_T(0); 00690 // delete old elements 00691 destroy_elements(tmp, old_datasize); 00692 } 00693 // if possible, reuse the allocated memory 00694 else if (datasize == rows * cols) { 00695 no_rows = rows; 00696 no_cols = cols; 00697 } 00698 // finally release old memory and allocate a new one 00699 else { 00700 free(); 00701 alloc(rows, cols); 00702 } 00703 } 00704 00705 template<class Num_T> inline 00706 void Mat<Num_T>::zeros() 00707 { 00708 for (int i = 0; i < datasize; i++) 00709 data[i] = Num_T(0); 00710 } 00711 00712 template<class Num_T> inline 00713 void Mat<Num_T>::ones() 00714 { 00715 for (int i = 0; i < datasize; i++) 00716 data[i] = Num_T(1); 00717 } 00718 00719 template<class Num_T> inline 00720 const Num_T& Mat<Num_T>::operator()(int r, int c) const 00721 { 00722 it_assert_debug(in_range(r, c), 00723 "Mat<>::operator(): Indexing out of range"); 00724 return data[r+c*no_rows]; 00725 } 00726 00727 template<class Num_T> inline 00728 Num_T& Mat<Num_T>::operator()(int r, int c) 00729 { 00730 it_assert_debug(in_range(r, c), 00731 "Mat<>::operator(): Indexing out of range"); 00732 return data[r+c*no_rows]; 00733 } 00734 00735 template<class Num_T> inline 00736 Num_T& Mat<Num_T>::operator()(int i) 00737 { 00738 it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range"); 00739 return data[i]; 00740 } 00741 00742 template<class Num_T> inline 00743 const Num_T& Mat<Num_T>::operator()(int i) const 00744 { 00745 it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range"); 00746 return data[i]; 00747 } 00748 00749 template<class Num_T> inline 00750 const Num_T& Mat<Num_T>::get(int r, int c) const 00751 { 00752 return (*this)(r, c); 00753 } 00754 00755 template<class Num_T> inline 00756 const Num_T& Mat<Num_T>::get(int i) const 00757 { 00758 return (*this)(i); 00759 } 00760 00761 template<class Num_T> inline 00762 void Mat<Num_T>::set(int r, int c, Num_T t) 00763 { 00764 it_assert_debug(in_range(r, c), "Mat<>::set(): Indexing out of range"); 00765 data[r+c*no_rows] = t; 00766 } 00767 00768 00769 template<class Num_T> 00770 void Mat<Num_T>::set(const std::string &str) 00771 { 00772 // actual row counter 00773 int rows = 0; 00774 // number of rows to allocate next time (8, 16, 32, 64, etc.) 00775 int maxrows = 8; 00776 00777 // clean the current matrix content 00778 free(); 00779 00780 // variable to store the start of a current vector 00781 std::string::size_type beg = 0; 00782 std::string::size_type end = 0; 00783 while (end != std::string::npos) { 00784 // find next occurrence of a semicolon in string str 00785 end = str.find(';', beg); 00786 // parse first row into a vector v 00787 Vec<Num_T> v(str.substr(beg, end - beg)); 00788 int v_size = v.size(); 00789 00790 // this check is necessary to parse the following two strings as the 00791 // same matrix: "1 0 1; ; 1 1; " and "1 0 1; 0 0 0; 1 1 0" 00792 if ((end != std::string::npos) || (v_size > 0)) { 00793 // matrix empty -> insert v as a first row and allocate maxrows 00794 if (rows == 0) { 00795 set_size(maxrows, v_size, true); 00796 set_row(rows++, v); 00797 } 00798 else { 00799 // check if we need to resize the matrix 00800 if ((rows == maxrows) || (v_size != no_cols)) { 00801 // we need to add new rows 00802 if (rows == maxrows) { 00803 maxrows *= 2; 00804 } 00805 // check if we need to add new columns 00806 if (v_size > no_cols) { 00807 set_size(maxrows, v_size, true); 00808 } 00809 else { 00810 set_size(maxrows, no_cols, true); 00811 // set the size of the parsed vector to the number of columns 00812 v.set_size(no_cols, true); 00813 } 00814 } 00815 // set the parsed vector as the next row 00816 set_row(rows++, v); 00817 } 00818 } 00819 // update the starting position of the next vector in the parsed 00820 // string 00821 beg = end + 1; 00822 } // if ((end != std::string::npos) || (v.size > 0)) 00823 00824 set_size(rows, no_cols, true); 00825 } 00826 00827 template<class Num_T> inline 00828 void Mat<Num_T>::set(const char *str) 00829 { 00830 set(std::string(str)); 00831 } 00832 00833 template<class Num_T> inline 00834 Mat<Num_T> Mat<Num_T>::operator()(int r1, int r2, int c1, int c2) const 00835 { 00836 if (r1 == -1) r1 = no_rows - 1; 00837 if (r2 == -1) r2 = no_rows - 1; 00838 if (c1 == -1) c1 = no_cols - 1; 00839 if (c2 == -1) c2 = no_cols - 1; 00840 00841 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) && 00842 (c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 00843 "Mat<>::operator()(r1, r2, c1, c2): Wrong indexing"); 00844 00845 Mat<Num_T> s(r2 - r1 + 1, c2 - c1 + 1); 00846 00847 for (int i = 0;i < s.no_cols;i++) 00848 copy_vector(s.no_rows, data + r1 + (c1 + i)*no_rows, s.data + i*s.no_rows); 00849 00850 return s; 00851 } 00852 00853 template<class Num_T> inline 00854 Mat<Num_T> Mat<Num_T>::get(int r1, int r2, int c1, int c2) const 00855 { 00856 return (*this)(r1, r2, c1, c2); 00857 } 00858 00859 template<class Num_T> inline 00860 Vec<Num_T> Mat<Num_T>::get_row(int r) const 00861 { 00862 it_assert_debug(row_in_range(r), "Mat<>::get_row(): Index out of range"); 00863 Vec<Num_T> a(no_cols); 00864 00865 copy_vector(no_cols, data + r, no_rows, a._data(), 1); 00866 return a; 00867 } 00868 00869 template<class Num_T> 00870 Mat<Num_T> Mat<Num_T>::get_rows(int r1, int r2) const 00871 { 00872 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows), 00873 "Mat<>::get_rows(): Wrong indexing"); 00874 Mat<Num_T> m(r2 - r1 + 1, no_cols); 00875 00876 for (int i = 0; i < m.rows(); i++) 00877 copy_vector(no_cols, data + i + r1, no_rows, m.data + i, m.no_rows); 00878 00879 return m; 00880 } 00881 00882 template<class Num_T> 00883 Mat<Num_T> Mat<Num_T>::get_rows(const Vec<int> &indexlist) const 00884 { 00885 Mat<Num_T> m(indexlist.size(), no_cols); 00886 00887 for (int i = 0;i < indexlist.size();i++) { 00888 it_assert_debug(row_in_range(indexlist(i)), 00889 "Mat<>::get_rows(indexlist): Indexing out of range"); 00890 copy_vector(no_cols, data + indexlist(i), no_rows, m.data + i, m.no_rows); 00891 } 00892 00893 return m; 00894 } 00895 00896 template<class Num_T> inline 00897 Vec<Num_T> Mat<Num_T>::get_col(int c) const 00898 { 00899 it_assert_debug(col_in_range(c), "Mat<>::get_col(): Index out of range"); 00900 Vec<Num_T> a(no_rows); 00901 00902 copy_vector(no_rows, data + c*no_rows, a._data()); 00903 00904 return a; 00905 } 00906 00907 template<class Num_T> 00908 Mat<Num_T> Mat<Num_T>::get_cols(int c1, int c2) const 00909 { 00910 it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 00911 "Mat<>::get_cols(): Wrong indexing"); 00912 Mat<Num_T> m(no_rows, c2 - c1 + 1); 00913 00914 for (int i = 0; i < m.cols(); i++) 00915 copy_vector(no_rows, data + (i + c1)*no_rows, m.data + i*m.no_rows); 00916 00917 return m; 00918 } 00919 00920 template<class Num_T> 00921 Mat<Num_T> Mat<Num_T>::get_cols(const Vec<int> &indexlist) const 00922 { 00923 Mat<Num_T> m(no_rows, indexlist.size()); 00924 00925 for (int i = 0; i < indexlist.size(); i++) { 00926 it_assert_debug(col_in_range(indexlist(i)), 00927 "Mat<>::get_cols(indexlist): Indexing out of range"); 00928 copy_vector(no_rows, data + indexlist(i)*no_rows, m.data + i*m.no_rows); 00929 } 00930 00931 return m; 00932 } 00933 00934 template<class Num_T> inline 00935 void Mat<Num_T>::set_row(int r, const Vec<Num_T> &v) 00936 { 00937 it_assert_debug(row_in_range(r), "Mat<>::set_row(): Index out of range"); 00938 it_assert_debug(v.size() == no_cols, 00939 "Mat<>::set_row(): Wrong size of input vector"); 00940 copy_vector(v.size(), v._data(), 1, data + r, no_rows); 00941 } 00942 00943 template<class Num_T> inline 00944 void Mat<Num_T>::set_col(int c, const Vec<Num_T> &v) 00945 { 00946 it_assert_debug(col_in_range(c), "Mat<>::set_col(): Index out of range"); 00947 it_assert_debug(v.size() == no_rows, 00948 "Mat<>::set_col(): Wrong size of input vector"); 00949 copy_vector(v.size(), v._data(), data + c*no_rows); 00950 } 00951 00952 00953 template<class Num_T> 00954 void Mat<Num_T>::set_rows(int r, const Mat<Num_T> &m) 00955 { 00956 it_assert_debug(row_in_range(r), "Mat<>::set_rows(): Index out of range"); 00957 it_assert_debug(no_cols == m.cols(), 00958 "Mat<>::set_rows(): Column sizes do not match"); 00959 it_assert_debug(m.rows() + r <= no_rows, 00960 "Mat<>::set_rows(): Not enough rows"); 00961 00962 for (int i = 0; i < m.rows(); ++i) { 00963 copy_vector(no_cols, m.data + i, m.no_rows, data + i + r, no_rows); 00964 } 00965 } 00966 00967 template<class Num_T> 00968 void Mat<Num_T>::set_cols(int c, const Mat<Num_T> &m) 00969 { 00970 it_assert_debug(col_in_range(c), "Mat<>::set_cols(): Index out of range"); 00971 it_assert_debug(no_rows == m.rows(), 00972 "Mat<>::set_cols(): Row sizes do not match"); 00973 it_assert_debug(m.cols() + c <= no_cols, 00974 "Mat<>::set_cols(): Not enough colums"); 00975 00976 for (int i = 0; i < m.cols(); ++i) { 00977 copy_vector(no_rows, m.data + i*no_rows, data + (i + c)*no_rows); 00978 } 00979 } 00980 00981 00982 template<class Num_T> inline 00983 void Mat<Num_T>::copy_row(int to, int from) 00984 { 00985 it_assert_debug(row_in_range(to) && row_in_range(from), 00986 "Mat<>::copy_row(): Indexing out of range"); 00987 if (from == to) 00988 return; 00989 00990 copy_vector(no_cols, data + from, no_rows, data + to, no_rows); 00991 } 00992 00993 template<class Num_T> inline 00994 void Mat<Num_T>::copy_col(int to, int from) 00995 { 00996 it_assert_debug(col_in_range(to) && col_in_range(from), 00997 "Mat<>::copy_col(): Indexing out of range"); 00998 if (from == to) 00999 return; 01000 01001 copy_vector(no_rows, data + from*no_rows, data + to*no_rows); 01002 } 01003 01004 template<class Num_T> inline 01005 void Mat<Num_T>::swap_rows(int r1, int r2) 01006 { 01007 it_assert_debug(row_in_range(r1) && row_in_range(r2), 01008 "Mat<>::swap_rows(): Indexing out of range"); 01009 if (r1 == r2) 01010 return; 01011 01012 swap_vector(no_cols, data + r1, no_rows, data + r2, no_rows); 01013 } 01014 01015 template<class Num_T> inline 01016 void Mat<Num_T>::swap_cols(int c1, int c2) 01017 { 01018 it_assert_debug(col_in_range(c1) && col_in_range(c2), 01019 "Mat<>::swap_cols(): Indexing out of range"); 01020 if (c1 == c2) 01021 return; 01022 01023 swap_vector(no_rows, data + c1*no_rows, data + c2*no_rows); 01024 } 01025 01026 template<class Num_T> 01027 void Mat<Num_T>::set_submatrix(int r1, int, int c1, int, const Mat<Num_T> &m) 01028 { 01029 it_warning("Mat<>::set_submatrix(r1, r2, r3, r4, m): This function is " 01030 "deprecated and might be removed from future IT++ releases. " 01031 "Please use Mat<>::set_submatrix(r, c, m) function instead."); 01032 set_submatrix(r1, c1, m); 01033 } 01034 01035 template<class Num_T> inline 01036 void Mat<Num_T>::set_submatrix(int r, int c, const Mat<Num_T> &m) 01037 { 01038 it_assert_debug((r >= 0) && (r + m.no_rows <= no_rows) && 01039 (c >= 0) && (c + m.no_cols <= no_cols), 01040 "Mat<>::set_submatrix(): Indexing out of range " 01041 "or wrong input matrix"); 01042 for (int i = 0; i < m.no_cols; i++) 01043 copy_vector(m.no_rows, m.data + i*m.no_rows, data + (c + i)*no_rows + r); 01044 } 01045 01046 01047 01048 template<class Num_T> inline 01049 void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, Num_T t) 01050 { 01051 if (r1 == -1) r1 = no_rows - 1; 01052 if (r2 == -1) r2 = no_rows - 1; 01053 if (c1 == -1) c1 = no_cols - 1; 01054 if (c2 == -1) c2 = no_cols - 1; 01055 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) && 01056 (c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 01057 "Mat<>::set_submatrix(): Wrong indexing"); 01058 for (int i = c1; i <= c2; i++) { 01059 int pos = i * no_rows + r1; 01060 for (int j = r1; j <= r2; j++) 01061 data[pos++] = t; 01062 } 01063 } 01064 01065 template<class Num_T> 01066 void Mat<Num_T>::del_row(int r) 01067 { 01068 it_assert_debug(row_in_range(r), "Mat<>::del_row(): Index out of range"); 01069 Mat<Num_T> Temp(*this); 01070 set_size(no_rows - 1, no_cols, false); 01071 for (int i = 0 ; i < r ; i++) { 01072 copy_vector(no_cols, &Temp.data[i], no_rows + 1, &data[i], no_rows); 01073 } 01074 for (int i = r ; i < no_rows ; i++) { 01075 copy_vector(no_cols, &Temp.data[i+1], no_rows + 1, &data[i], no_rows); 01076 } 01077 01078 } 01079 01080 template<class Num_T> 01081 void Mat<Num_T>::del_rows(int r1, int r2) 01082 { 01083 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows), 01084 "Mat<>::del_rows(): Indexing out of range"); 01085 Mat<Num_T> Temp(*this); 01086 int no_del_rows = r2 - r1 + 1; 01087 set_size(no_rows - no_del_rows, no_cols, false); 01088 for (int i = 0; i < r1 ; ++i) { 01089 copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows); 01090 } 01091 for (int i = r2 + 1; i < Temp.no_rows; ++i) { 01092 copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i-no_del_rows], 01093 no_rows); 01094 } 01095 } 01096 01097 template<class Num_T> 01098 void Mat<Num_T>::del_col(int c) 01099 { 01100 it_assert_debug(col_in_range(c), "Mat<>::del_col(): Index out of range"); 01101 Mat<Num_T> Temp(*this); 01102 01103 set_size(no_rows, no_cols - 1, false); 01104 copy_vector(c*no_rows, Temp.data, data); 01105 copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]); 01106 } 01107 01108 template<class Num_T> 01109 void Mat<Num_T>::del_cols(int c1, int c2) 01110 { 01111 it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 01112 "Mat<>::del_cols(): Indexing out of range"); 01113 Mat<Num_T> Temp(*this); 01114 int n_deleted_cols = c2 - c1 + 1; 01115 set_size(no_rows, no_cols - n_deleted_cols, false); 01116 copy_vector(c1*no_rows, Temp.data, data); 01117 copy_vector((no_cols - c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]); 01118 } 01119 01120 template<class Num_T> 01121 void Mat<Num_T>::ins_row(int r, const Vec<Num_T> &v) 01122 { 01123 it_assert_debug((r >= 0) && (r <= no_rows), 01124 "Mat<>::ins_row(): Index out of range"); 01125 it_assert_debug((v.size() == no_cols) || (no_rows == 0), 01126 "Mat<>::ins_row(): Wrong size of the input vector"); 01127 01128 if (no_cols == 0) { 01129 no_cols = v.size(); 01130 } 01131 01132 Mat<Num_T> Temp(*this); 01133 set_size(no_rows + 1, no_cols, false); 01134 01135 for (int i = 0 ; i < r ; i++) { 01136 copy_vector(no_cols, &Temp.data[i], no_rows - 1, &data[i], no_rows); 01137 } 01138 copy_vector(no_cols, v._data(), 1, &data[r], no_rows); 01139 for (int i = r + 1 ; i < no_rows ; i++) { 01140 copy_vector(no_cols, &Temp.data[i-1], no_rows - 1, &data[i], no_rows); 01141 } 01142 } 01143 01144 template<class Num_T> 01145 void Mat<Num_T>::ins_col(int c, const Vec<Num_T> &v) 01146 { 01147 it_assert_debug((c >= 0) && (c <= no_cols), 01148 "Mat<>::ins_col(): Index out of range"); 01149 it_assert_debug((v.size() == no_rows) || (no_cols == 0), 01150 "Mat<>::ins_col(): Wrong size of the input vector"); 01151 01152 if (no_rows == 0) { 01153 no_rows = v.size(); 01154 } 01155 01156 Mat<Num_T> Temp(*this); 01157 set_size(no_rows, no_cols + 1, false); 01158 01159 copy_vector(c*no_rows, Temp.data, data); 01160 copy_vector(no_rows, v._data(), &data[c*no_rows]); 01161 copy_vector((no_cols - c - 1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]); 01162 } 01163 01164 template<class Num_T> inline 01165 void Mat<Num_T>::append_row(const Vec<Num_T> &v) 01166 { 01167 ins_row(no_rows, v); 01168 } 01169 01170 template<class Num_T> inline 01171 void Mat<Num_T>::append_col(const Vec<Num_T> &v) 01172 { 01173 ins_col(no_cols, v); 01174 } 01175 01176 template<class Num_T> 01177 Mat<Num_T> Mat<Num_T>::transpose() const 01178 { 01179 Mat<Num_T> temp(no_cols, no_rows); 01180 for (int i = 0; i < no_rows; ++i) { 01181 copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1); 01182 } 01183 return temp; 01184 } 01185 01187 template<> 01188 cmat cmat::hermitian_transpose() const; 01190 01191 template<class Num_T> 01192 Mat<Num_T> Mat<Num_T>::hermitian_transpose() const 01193 { 01194 Mat<Num_T> temp(no_cols, no_rows); 01195 for (int i = 0; i < no_rows; ++i) { 01196 copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1); 01197 } 01198 return temp; 01199 } 01200 01201 template<class Num_T> 01202 Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01203 { 01204 // if one of the input matrix is empty just copy the other one as a result 01205 if (m1.no_cols == 0) 01206 return m2; 01207 if (m2.no_cols == 0) 01208 return m1; 01209 it_assert_debug(m1.no_rows == m2.no_rows, 01210 "Mat<>::concat_horizontal(): Wrong sizes"); 01211 int no_rows = m1.no_rows; 01212 Mat<Num_T> temp(no_rows, m1.no_cols + m2.no_cols); 01213 for (int i = 0; i < m1.no_cols; ++i) { 01214 copy_vector(no_rows, &m1.data[i * no_rows], &temp.data[i * no_rows]); 01215 } 01216 for (int i = 0; i < m2.no_cols; ++i) { 01217 copy_vector(no_rows, &m2.data[i * no_rows], &temp.data[(m1.no_cols + i) 01218 * no_rows]); 01219 } 01220 return temp; 01221 } 01222 01223 template<class Num_T> 01224 Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01225 { 01226 // if one of the input matrix is empty just copy the other one as a result 01227 if (m1.no_rows == 0) 01228 return m2; 01229 if (m2.no_rows == 0) 01230 return m1; 01231 it_assert_debug(m1.no_cols == m2.no_cols, 01232 "Mat<>::concat_vertical(): Wrong sizes"); 01233 int no_cols = m1.no_cols; 01234 Mat<Num_T> temp(m1.no_rows + m2.no_rows, no_cols); 01235 for (int i = 0; i < no_cols; ++i) { 01236 copy_vector(m1.no_rows, &m1.data[i * m1.no_rows], 01237 &temp.data[i * temp.no_rows]); 01238 copy_vector(m2.no_rows, &m2.data[i * m2.no_rows], 01239 &temp.data[i * temp.no_rows + m1.no_rows]); 01240 } 01241 return temp; 01242 } 01243 01244 template<class Num_T> inline 01245 Mat<Num_T>& Mat<Num_T>::operator=(Num_T t) 01246 { 01247 for (int i = 0; i < datasize; i++) 01248 data[i] = t; 01249 return *this; 01250 } 01251 01252 template<class Num_T> inline 01253 Mat<Num_T>& Mat<Num_T>::operator=(const Mat<Num_T> &m) 01254 { 01255 if (this != &m) { 01256 set_size(m.no_rows, m.no_cols, false); 01257 if (m.datasize != 0) 01258 copy_vector(m.datasize, m.data, data); 01259 } 01260 return *this; 01261 } 01262 01263 template<class Num_T> inline 01264 Mat<Num_T>& Mat<Num_T>::operator=(const Vec<Num_T> &v) 01265 { 01266 it_assert_debug(((no_rows == 1) && (no_cols == v.size())) 01267 || ((no_cols == 1) && (no_rows == v.size())), 01268 "Mat<>::operator=(): Wrong size of the input vector"); 01269 set_size(v.size(), 1, false); 01270 copy_vector(v.size(), v._data(), data); 01271 return *this; 01272 } 01273 01274 template<class Num_T> inline 01275 Mat<Num_T>& Mat<Num_T>::operator=(const std::string &str) 01276 { 01277 set(str); 01278 return *this; 01279 } 01280 01281 template<class Num_T> inline 01282 Mat<Num_T>& Mat<Num_T>::operator=(const char *str) 01283 { 01284 set(std::string(str)); 01285 return *this; 01286 } 01287 01288 //-------------------- Templated friend functions -------------------------- 01289 01290 template<class Num_T> 01291 Mat<Num_T>& Mat<Num_T>::operator+=(const Mat<Num_T> &m) 01292 { 01293 if (datasize == 0) 01294 operator=(m); 01295 else { 01296 int i, j, m_pos = 0, pos = 0; 01297 it_assert_debug(m.no_rows == no_rows && m.no_cols == no_cols, "Mat<Num_T>::operator+=: wrong sizes"); 01298 for (i = 0; i < no_cols; i++) { 01299 for (j = 0; j < no_rows; j++) 01300 data[pos+j] += m.data[m_pos+j]; 01301 pos += no_rows; 01302 m_pos += m.no_rows; 01303 } 01304 } 01305 return *this; 01306 } 01307 01308 template<class Num_T> inline 01309 Mat<Num_T>& Mat<Num_T>::operator+=(Num_T t) 01310 { 01311 for (int i = 0; i < datasize; i++) 01312 data[i] += t; 01313 return *this; 01314 } 01315 01316 template<class Num_T> 01317 Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01318 { 01319 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01320 int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0; 01321 01322 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01323 "Mat<>::operator+(): Wrong sizes"); 01324 01325 for (i = 0; i < r.no_cols; i++) { 01326 for (j = 0; j < r.no_rows; j++) 01327 r.data[r_pos+j] = m1.data[m1_pos+j] + m2.data[m2_pos+j]; 01328 // next column 01329 m1_pos += m1.no_rows; 01330 m2_pos += m2.no_rows; 01331 r_pos += r.no_rows; 01332 } 01333 01334 return r; 01335 } 01336 01337 01338 template<class Num_T> 01339 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t) 01340 { 01341 Mat<Num_T> r(m.no_rows, m.no_cols); 01342 01343 for (int i = 0; i < r.datasize; i++) 01344 r.data[i] = m.data[i] + t; 01345 01346 return r; 01347 } 01348 01349 template<class Num_T> 01350 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m) 01351 { 01352 Mat<Num_T> r(m.no_rows, m.no_cols); 01353 01354 for (int i = 0; i < r.datasize; i++) 01355 r.data[i] = t + m.data[i]; 01356 01357 return r; 01358 } 01359 01360 template<class Num_T> 01361 Mat<Num_T>& Mat<Num_T>::operator-=(const Mat<Num_T> &m) 01362 { 01363 int i, j, m_pos = 0, pos = 0; 01364 01365 if (datasize == 0) { 01366 set_size(m.no_rows, m.no_cols, false); 01367 for (i = 0; i < no_cols; i++) { 01368 for (j = 0; j < no_rows; j++) 01369 data[pos+j] = -m.data[m_pos+j]; 01370 // next column 01371 m_pos += m.no_rows; 01372 pos += no_rows; 01373 } 01374 } 01375 else { 01376 it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols), 01377 "Mat<>::operator-=(): Wrong sizes"); 01378 for (i = 0; i < no_cols; i++) { 01379 for (j = 0; j < no_rows; j++) 01380 data[pos+j] -= m.data[m_pos+j]; 01381 // next column 01382 m_pos += m.no_rows; 01383 pos += no_rows; 01384 } 01385 } 01386 return *this; 01387 } 01388 01389 template<class Num_T> 01390 Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01391 { 01392 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01393 int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0; 01394 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01395 "Mat<>::operator-(): Wrong sizes"); 01396 01397 for (i = 0; i < r.no_cols; i++) { 01398 for (j = 0; j < r.no_rows; j++) 01399 r.data[r_pos+j] = m1.data[m1_pos+j] - m2.data[m2_pos+j]; 01400 // next column 01401 m1_pos += m1.no_rows; 01402 m2_pos += m2.no_rows; 01403 r_pos += r.no_rows; 01404 } 01405 01406 return r; 01407 } 01408 01409 template<class Num_T> inline 01410 Mat<Num_T>& Mat<Num_T>::operator-=(Num_T t) 01411 { 01412 for (int i = 0; i < datasize; i++) 01413 data[i] -= t; 01414 return *this; 01415 } 01416 01417 template<class Num_T> 01418 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t) 01419 { 01420 Mat<Num_T> r(m.no_rows, m.no_cols); 01421 int i, j, m_pos = 0, r_pos = 0; 01422 01423 for (i = 0; i < r.no_cols; i++) { 01424 for (j = 0; j < r.no_rows; j++) 01425 r.data[r_pos+j] = m.data[m_pos+j] - t; 01426 // next column 01427 m_pos += m.no_rows; 01428 r_pos += r.no_rows; 01429 } 01430 01431 return r; 01432 } 01433 01434 template<class Num_T> 01435 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m) 01436 { 01437 Mat<Num_T> r(m.no_rows, m.no_cols); 01438 int i, j, m_pos = 0, r_pos = 0; 01439 01440 for (i = 0; i < r.no_cols; i++) { 01441 for (j = 0; j < r.no_rows; j++) 01442 r.data[r_pos+j] = t - m.data[m_pos+j]; 01443 // next column 01444 m_pos += m.no_rows; 01445 r_pos += r.no_rows; 01446 } 01447 01448 return r; 01449 } 01450 01451 template<class Num_T> 01452 Mat<Num_T> operator-(const Mat<Num_T> &m) 01453 { 01454 Mat<Num_T> r(m.no_rows, m.no_cols); 01455 int i, j, m_pos = 0, r_pos = 0; 01456 01457 for (i = 0; i < r.no_cols; i++) { 01458 for (j = 0; j < r.no_rows; j++) 01459 r.data[r_pos+j] = -m.data[m_pos+j]; 01460 // next column 01461 m_pos += m.no_rows; 01462 r_pos += r.no_rows; 01463 } 01464 01465 return r; 01466 } 01467 01468 01469 template<> mat& mat::operator*=(const mat &m); 01470 template<> cmat& cmat::operator*=(const cmat &m); 01471 01472 template<class Num_T> 01473 Mat<Num_T>& Mat<Num_T>::operator*=(const Mat<Num_T> &m) 01474 { 01475 it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes"); 01476 Mat<Num_T> r(no_rows, m.no_cols); 01477 01478 Num_T tmp; 01479 01480 int i, j, k, r_pos = 0, pos = 0, m_pos = 0; 01481 01482 for (i = 0; i < r.no_cols; i++) { 01483 for (j = 0; j < r.no_rows; j++) { 01484 tmp = Num_T(0); 01485 pos = 0; 01486 for (k = 0; k < no_cols; k++) { 01487 tmp += data[pos+j] * m.data[m_pos+k]; 01488 pos += no_rows; 01489 } 01490 r.data[r_pos+j] = tmp; 01491 } 01492 r_pos += r.no_rows; 01493 m_pos += m.no_rows; 01494 } 01495 operator=(r); // time consuming 01496 return *this; 01497 } 01498 01499 template<class Num_T> inline 01500 Mat<Num_T>& Mat<Num_T>::operator*=(Num_T t) 01501 { 01502 scal_vector(datasize, t, data); 01503 return *this; 01504 } 01505 01506 01507 template<> mat operator*(const mat &m1, const mat &m2); 01508 template<> cmat operator*(const cmat &m1, const cmat &m2); 01509 01510 template<class Num_T> 01511 Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01512 { 01513 it_assert_debug(m1.no_cols == m2.no_rows, 01514 "Mat<>::operator*(): Wrong sizes"); 01515 Mat<Num_T> r(m1.no_rows, m2.no_cols); 01516 01517 Num_T tmp; 01518 int i, j, k; 01519 Num_T *tr = r.data, *t1, *t2 = m2.data; 01520 01521 for (i = 0; i < r.no_cols; i++) { 01522 for (j = 0; j < r.no_rows; j++) { 01523 tmp = Num_T(0); 01524 t1 = m1.data + j; 01525 for (k = m1.no_cols; k > 0; k--) { 01526 tmp += *(t1) * *(t2++); 01527 t1 += m1.no_rows; 01528 } 01529 *(tr++) = tmp; 01530 t2 -= m2.no_rows; 01531 } 01532 t2 += m2.no_rows; 01533 } 01534 01535 return r; 01536 } 01537 01538 01539 template<> vec operator*(const mat &m, const vec &v); 01540 template<> cvec operator*(const cmat &m, const cvec &v); 01541 01542 template<class Num_T> 01543 Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v) 01544 { 01545 it_assert_debug(m.no_cols == v.size(), 01546 "Mat<>::operator*(): Wrong sizes"); 01547 Vec<Num_T> r(m.no_rows); 01548 int i, k, m_pos; 01549 01550 for (i = 0; i < m.no_rows; i++) { 01551 r(i) = Num_T(0); 01552 m_pos = 0; 01553 for (k = 0; k < m.no_cols; k++) { 01554 r(i) += m.data[m_pos+i] * v(k); 01555 m_pos += m.no_rows; 01556 } 01557 } 01558 01559 return r; 01560 } 01561 01562 template<class Num_T> 01563 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t) 01564 { 01565 Mat<Num_T> r(m.no_rows, m.no_cols); 01566 01567 for (int i = 0; i < r.datasize; i++) 01568 r.data[i] = m.data[i] * t; 01569 01570 return r; 01571 } 01572 01573 template<class Num_T> inline 01574 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m) 01575 { 01576 return operator*(m, t); 01577 } 01578 01579 template<class Num_T> inline 01580 Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01581 { 01582 Mat<Num_T> out; 01583 elem_mult_out(m1, m2, out); 01584 return out; 01585 } 01586 01587 template<class Num_T> 01588 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01589 Mat<Num_T> &out) 01590 { 01591 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01592 "Mat<>::elem_mult_out(): Wrong sizes"); 01593 out.set_size(m1.no_rows, m1.no_cols); 01594 for (int i = 0; i < out.datasize; i++) 01595 out.data[i] = m1.data[i] * m2.data[i]; 01596 } 01597 01598 template<class Num_T> 01599 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01600 const Mat<Num_T> &m3, Mat<Num_T> &out) 01601 { 01602 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows) 01603 && (m1.no_cols == m2.no_cols) && (m1.no_cols == m3.no_cols), 01604 "Mat<>::elem_mult_out(): Wrong sizes"); 01605 out.set_size(m1.no_rows, m1.no_cols); 01606 for (int i = 0; i < out.datasize; i++) 01607 out.data[i] = m1.data[i] * m2.data[i] * m3.data[i]; 01608 } 01609 01610 template<class Num_T> 01611 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01612 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 01613 Mat<Num_T> &out) 01614 { 01615 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows) 01616 && (m1.no_rows == m4.no_rows) && (m1.no_cols == m2.no_cols) 01617 && (m1.no_cols == m3.no_cols) && (m1.no_cols == m4.no_cols), 01618 "Mat<>::elem_mult_out(): Wrong sizes"); 01619 out.set_size(m1.no_rows, m1.no_cols); 01620 for (int i = 0; i < out.datasize; i++) 01621 out.data[i] = m1.data[i] * m2.data[i] * m3.data[i] * m4.data[i]; 01622 } 01623 01624 template<class Num_T> 01625 #ifndef _MSC_VER 01626 inline 01627 #endif 01628 void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2) 01629 { 01630 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01631 "Mat<>::elem_mult_inplace(): Wrong sizes"); 01632 for (int i = 0; i < m2.datasize; i++) 01633 m2.data[i] *= m1.data[i]; 01634 } 01635 01636 template<class Num_T> inline 01637 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01638 { 01639 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01640 "Mat<>::elem_mult_sum(): Wrong sizes"); 01641 Num_T acc = 0; 01642 01643 for (int i = 0; i < m1.datasize; i++) 01644 acc += m1.data[i] * m2.data[i]; 01645 01646 return acc; 01647 } 01648 01649 template<class Num_T> inline 01650 Mat<Num_T>& Mat<Num_T>::operator/=(Num_T t) 01651 { 01652 for (int i = 0; i < datasize; i++) 01653 data[i] /= t; 01654 return *this; 01655 } 01656 01657 template<class Num_T> inline 01658 Mat<Num_T>& Mat<Num_T>::operator/=(const Mat<Num_T> &m) 01659 { 01660 it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols), 01661 "Mat<>::operator/=(): Wrong sizes"); 01662 for (int i = 0; i < datasize; i++) 01663 data[i] /= m.data[i]; 01664 return *this; 01665 } 01666 01667 template<class Num_T> 01668 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t) 01669 { 01670 Mat<Num_T> r(m.no_rows, m.no_cols); 01671 for (int i = 0; i < r.datasize; ++i) 01672 r.data[i] = m.data[i] / t; 01673 return r; 01674 } 01675 01676 template<class Num_T> 01677 Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m) 01678 { 01679 Mat<Num_T> r(m.no_rows, m.no_cols); 01680 for (int i = 0; i < r.datasize; ++i) 01681 r.data[i] = t / m.data[i]; 01682 return r; 01683 } 01684 01685 template<class Num_T> inline 01686 Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01687 { 01688 Mat<Num_T> out; 01689 elem_div_out(m1, m2, out); 01690 return out; 01691 } 01692 01693 template<class Num_T> 01694 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01695 Mat<Num_T> &out) 01696 { 01697 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01698 "Mat<>::elem_div_out(): Wrong sizes"); 01699 01700 if ((out.no_rows != m1.no_rows) || (out.no_cols != m1.no_cols)) 01701 out.set_size(m1.no_rows, m1.no_cols); 01702 01703 for (int i = 0; i < out.datasize; i++) 01704 out.data[i] = m1.data[i] / m2.data[i]; 01705 } 01706 01707 template<class Num_T> inline 01708 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01709 { 01710 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01711 "Mat<>::elem_div_sum(): Wrong sizes"); 01712 Num_T acc = 0; 01713 01714 for (int i = 0; i < m1.datasize; i++) 01715 acc += m1.data[i] / m2.data[i]; 01716 01717 return acc; 01718 } 01719 01720 template<class Num_T> 01721 bool Mat<Num_T>::operator==(const Mat<Num_T> &m) const 01722 { 01723 if (no_rows != m.no_rows || no_cols != m.no_cols) return false; 01724 for (int i = 0;i < datasize;i++) { 01725 if (data[i] != m.data[i]) return false; 01726 } 01727 return true; 01728 } 01729 01730 template<class Num_T> 01731 bool Mat<Num_T>::operator!=(const Mat<Num_T> &m) const 01732 { 01733 if (no_rows != m.no_rows || no_cols != m.no_cols) return true; 01734 for (int i = 0;i < datasize;i++) { 01735 if (data[i] != m.data[i]) return true; 01736 } 01737 return false; 01738 } 01739 01740 template <class Num_T> 01741 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m) 01742 { 01743 int i; 01744 01745 switch (m.rows()) { 01746 case 0 : 01747 os << "[]"; 01748 break; 01749 case 1 : 01750 os << '[' << m.get_row(0) << ']'; 01751 break; 01752 default: 01753 os << '[' << m.get_row(0) << std::endl; 01754 for (i = 1; i < m.rows() - 1; i++) 01755 os << ' ' << m.get_row(i) << std::endl; 01756 os << ' ' << m.get_row(m.rows() - 1) << ']'; 01757 } 01758 01759 return os; 01760 } 01761 01762 template <class Num_T> 01763 std::istream &operator>>(std::istream &is, Mat<Num_T> &m) 01764 { 01765 std::ostringstream buffer; 01766 bool started = false; 01767 bool finished = false; 01768 bool brackets = false; 01769 bool within_double_brackets = false; 01770 char c; 01771 01772 while (!finished) { 01773 if (is.eof()) { 01774 finished = true; 01775 } 01776 else { 01777 is.get(c); 01778 01779 if (is.eof() || (c == '\n')) { 01780 if (brackets) { 01781 // Right bracket missing 01782 is.setstate(std::ios_base::failbit); 01783 finished = true; 01784 } 01785 else if (!((c == '\n') && !started)) { 01786 finished = true; 01787 } 01788 } 01789 else if ((c == ' ') || (c == '\t')) { 01790 if (started) { 01791 buffer << ' '; 01792 } 01793 } 01794 else if (c == '[') { 01795 if ((started && !brackets) || within_double_brackets) { 01796 // Unexpected left bracket 01797 is.setstate(std::ios_base::failbit); 01798 finished = true; 01799 } 01800 else if (!started) { 01801 started = true; 01802 brackets = true; 01803 } 01804 else { 01805 within_double_brackets = true; 01806 } 01807 } 01808 else if (c == ']') { 01809 if (!started || !brackets) { 01810 // Unexpected right bracket 01811 is.setstate(std::ios_base::failbit); 01812 finished = true; 01813 } 01814 else if (within_double_brackets) { 01815 within_double_brackets = false; 01816 buffer << ';'; 01817 } 01818 else { 01819 finished = true; 01820 } 01821 while (!is.eof() && (((c = static_cast<char>(is.peek())) == ' ') 01822 || (c == '\t'))) { 01823 is.get(); 01824 } 01825 if (!is.eof() && (c == '\n')) { 01826 is.get(); 01827 } 01828 } 01829 else { 01830 started = true; 01831 buffer << c; 01832 } 01833 } 01834 } 01835 01836 if (!started) { 01837 m.set_size(0, false); 01838 } 01839 else { 01840 m.set(buffer.str()); 01841 } 01842 01843 return is; 01844 } 01845 01847 01848 // --------------------------------------------------------------------- 01849 // Instantiations 01850 // --------------------------------------------------------------------- 01851 01852 #ifndef _MSC_VER 01853 01854 // class instantiations 01855 01856 extern template class Mat<double>; 01857 extern template class Mat<std::complex<double> >; 01858 extern template class Mat<int>; 01859 extern template class Mat<short int>; 01860 extern template class Mat<bin>; 01861 01862 // addition operators 01863 01864 extern template mat operator+(const mat &m1, const mat &m2); 01865 extern template cmat operator+(const cmat &m1, const cmat &m2); 01866 extern template imat operator+(const imat &m1, const imat &m2); 01867 extern template smat operator+(const smat &m1, const smat &m2); 01868 extern template bmat operator+(const bmat &m1, const bmat &m2); 01869 01870 extern template mat operator+(const mat &m, double t); 01871 extern template cmat operator+(const cmat &m, std::complex<double> t); 01872 extern template imat operator+(const imat &m, int t); 01873 extern template smat operator+(const smat &m, short t); 01874 extern template bmat operator+(const bmat &m, bin t); 01875 01876 extern template mat operator+(double t, const mat &m); 01877 extern template cmat operator+(std::complex<double> t, const cmat &m); 01878 extern template imat operator+(int t, const imat &m); 01879 extern template smat operator+(short t, const smat &m); 01880 extern template bmat operator+(bin t, const bmat &m); 01881 01882 // subtraction operators 01883 01884 extern template mat operator-(const mat &m1, const mat &m2); 01885 extern template cmat operator-(const cmat &m1, const cmat &m2); 01886 extern template imat operator-(const imat &m1, const imat &m2); 01887 extern template smat operator-(const smat &m1, const smat &m2); 01888 extern template bmat operator-(const bmat &m1, const bmat &m2); 01889 01890 extern template mat operator-(const mat &m, double t); 01891 extern template cmat operator-(const cmat &m, std::complex<double> t); 01892 extern template imat operator-(const imat &m, int t); 01893 extern template smat operator-(const smat &m, short t); 01894 extern template bmat operator-(const bmat &m, bin t); 01895 01896 extern template mat operator-(double t, const mat &m); 01897 extern template cmat operator-(std::complex<double> t, const cmat &m); 01898 extern template imat operator-(int t, const imat &m); 01899 extern template smat operator-(short t, const smat &m); 01900 extern template bmat operator-(bin t, const bmat &m); 01901 01902 // unary minus 01903 01904 extern template mat operator-(const mat &m); 01905 extern template cmat operator-(const cmat &m); 01906 extern template imat operator-(const imat &m); 01907 extern template smat operator-(const smat &m); 01908 extern template bmat operator-(const bmat &m); 01909 01910 // multiplication operators 01911 01912 extern template imat operator*(const imat &m1, const imat &m2); 01913 extern template smat operator*(const smat &m1, const smat &m2); 01914 extern template bmat operator*(const bmat &m1, const bmat &m2); 01915 01916 extern template ivec operator*(const imat &m, const ivec &v); 01917 extern template svec operator*(const smat &m, const svec &v); 01918 extern template bvec operator*(const bmat &m, const bvec &v); 01919 01920 extern template mat operator*(const mat &m, double t); 01921 extern template cmat operator*(const cmat &m, std::complex<double> t); 01922 extern template imat operator*(const imat &m, int t); 01923 extern template smat operator*(const smat &m, short t); 01924 extern template bmat operator*(const bmat &m, bin t); 01925 01926 extern template mat operator*(double t, const mat &m); 01927 extern template cmat operator*(std::complex<double> t, const cmat &m); 01928 extern template imat operator*(int t, const imat &m); 01929 extern template smat operator*(short t, const smat &m); 01930 extern template bmat operator*(bin t, const bmat &m); 01931 01932 // element-wise multiplication 01933 01934 extern template mat elem_mult(const mat &m1, const mat &m2); 01935 extern template cmat elem_mult(const cmat &m1, const cmat &m2); 01936 extern template imat elem_mult(const imat &m1, const imat &m2); 01937 extern template smat elem_mult(const smat &m1, const smat &m2); 01938 extern template bmat elem_mult(const bmat &m1, const bmat &m2); 01939 01940 extern template void elem_mult_out(const mat &m1, const mat &m2, mat &out); 01941 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01942 cmat &out); 01943 extern template void elem_mult_out(const imat &m1, const imat &m2, 01944 imat &out); 01945 extern template void elem_mult_out(const smat &m1, const smat &m2, 01946 smat &out); 01947 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01948 bmat &out); 01949 01950 extern template void elem_mult_out(const mat &m1, const mat &m2, 01951 const mat &m3, mat &out); 01952 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01953 const cmat &m3, cmat &out); 01954 extern template void elem_mult_out(const imat &m1, const imat &m2, 01955 const imat &m3, imat &out); 01956 extern template void elem_mult_out(const smat &m1, const smat &m2, 01957 const smat &m3, smat &out); 01958 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01959 const bmat &m3, bmat &out); 01960 01961 extern template void elem_mult_out(const mat &m1, const mat &m2, 01962 const mat &m3, const mat &m4, mat &out); 01963 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01964 const cmat &m3, const cmat &m4, 01965 cmat &out); 01966 extern template void elem_mult_out(const imat &m1, const imat &m2, 01967 const imat &m3, const imat &m4, 01968 imat &out); 01969 extern template void elem_mult_out(const smat &m1, const smat &m2, 01970 const smat &m3, const smat &m4, 01971 smat &out); 01972 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01973 const bmat &m3, const bmat &m4, 01974 bmat &out); 01975 01976 extern template void elem_mult_inplace(const mat &m1, mat &m2); 01977 extern template void elem_mult_inplace(const cmat &m1, cmat &m2); 01978 extern template void elem_mult_inplace(const imat &m1, imat &m2); 01979 extern template void elem_mult_inplace(const smat &m1, smat &m2); 01980 extern template void elem_mult_inplace(const bmat &m1, bmat &m2); 01981 01982 extern template double elem_mult_sum(const mat &m1, const mat &m2); 01983 extern template std::complex<double> elem_mult_sum(const cmat &m1, 01984 const cmat &m2); 01985 extern template int elem_mult_sum(const imat &m1, const imat &m2); 01986 extern template short elem_mult_sum(const smat &m1, const smat &m2); 01987 extern template bin elem_mult_sum(const bmat &m1, const bmat &m2); 01988 01989 // division operator 01990 01991 extern template mat operator/(double t, const mat &m); 01992 extern template cmat operator/(std::complex<double> t, const cmat &m); 01993 extern template imat operator/(int t, const imat &m); 01994 extern template smat operator/(short t, const smat &m); 01995 extern template bmat operator/(bin t, const bmat &m); 01996 01997 extern template mat operator/(const mat &m, double t); 01998 extern template cmat operator/(const cmat &m, std::complex<double> t); 01999 extern template imat operator/(const imat &m, int t); 02000 extern template smat operator/(const smat &m, short t); 02001 extern template bmat operator/(const bmat &m, bin t); 02002 02003 // element-wise division 02004 02005 extern template mat elem_div(const mat &m1, const mat &m2); 02006 extern template cmat elem_div(const cmat &m1, const cmat &m2); 02007 extern template imat elem_div(const imat &m1, const imat &m2); 02008 extern template smat elem_div(const smat &m1, const smat &m2); 02009 extern template bmat elem_div(const bmat &m1, const bmat &m2); 02010 02011 extern template void elem_div_out(const mat &m1, const mat &m2, mat &out); 02012 extern template void elem_div_out(const cmat &m1, const cmat &m2, cmat &out); 02013 extern template void elem_div_out(const imat &m1, const imat &m2, imat &out); 02014 extern template void elem_div_out(const smat &m1, const smat &m2, smat &out); 02015 extern template void elem_div_out(const bmat &m1, const bmat &m2, bmat &out); 02016 02017 extern template double elem_div_sum(const mat &m1, const mat &m2); 02018 extern template std::complex<double> elem_div_sum(const cmat &m1, 02019 const cmat &m2); 02020 extern template int elem_div_sum(const imat &m1, const imat &m2); 02021 extern template short elem_div_sum(const smat &m1, const smat &m2); 02022 extern template bin elem_div_sum(const bmat &m1, const bmat &m2); 02023 02024 // concatenation 02025 02026 extern template mat concat_horizontal(const mat &m1, const mat &m2); 02027 extern template cmat concat_horizontal(const cmat &m1, const cmat &m2); 02028 extern template imat concat_horizontal(const imat &m1, const imat &m2); 02029 extern template smat concat_horizontal(const smat &m1, const smat &m2); 02030 extern template bmat concat_horizontal(const bmat &m1, const bmat &m2); 02031 02032 extern template mat concat_vertical(const mat &m1, const mat &m2); 02033 extern template cmat concat_vertical(const cmat &m1, const cmat &m2); 02034 extern template imat concat_vertical(const imat &m1, const imat &m2); 02035 extern template smat concat_vertical(const smat &m1, const smat &m2); 02036 extern template bmat concat_vertical(const bmat &m1, const bmat &m2); 02037 02038 // I/O streams 02039 02040 extern template std::ostream &operator<<(std::ostream &os, const mat &m); 02041 extern template std::ostream &operator<<(std::ostream &os, const cmat &m); 02042 extern template std::ostream &operator<<(std::ostream &os, const imat &m); 02043 extern template std::ostream &operator<<(std::ostream &os, const smat &m); 02044 extern template std::ostream &operator<<(std::ostream &os, const bmat &m); 02045 02046 extern template std::istream &operator>>(std::istream &is, mat &m); 02047 extern template std::istream &operator>>(std::istream &is, cmat &m); 02048 extern template std::istream &operator>>(std::istream &is, imat &m); 02049 extern template std::istream &operator>>(std::istream &is, smat &m); 02050 extern template std::istream &operator>>(std::istream &is, bmat &m); 02051 02052 #endif // _MSC_VER 02053 02055 02056 } // namespace itpp 02057 02058 #endif // #ifndef MAT_H
Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4