00001 00030 #ifndef MAT_H 00031 #define MAT_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/itassert.h> 00040 #include <itpp/base/math/misc.h> 00041 #include <itpp/base/factory.h> 00042 00043 00044 namespace itpp { 00045 00046 // Declaration of Vec 00047 template<class Num_T> class Vec; 00048 // Declaration of Mat 00049 template<class Num_T> class Mat; 00050 // Declaration of bin 00051 class bin; 00052 00054 template<class Num_T> 00055 Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00057 template<class Num_T> 00058 Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00059 00061 template<class Num_T> 00062 Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00064 template<class Num_T> 00065 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t); 00067 template<class Num_T> 00068 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m); 00069 00071 template<class Num_T> 00072 Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00074 template<class Num_T> 00075 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t); 00077 template<class Num_T> 00078 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m); 00080 template<class Num_T> 00081 Mat<Num_T> operator-(const Mat<Num_T> &m); 00082 00084 template<class Num_T> 00085 Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00087 template<class Num_T> 00088 Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v); 00090 template<class Num_T> 00091 Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m); 00093 template<class Num_T> 00094 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t); 00096 template<class Num_T> 00097 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m); 00098 00100 template<class Num_T> 00101 Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00103 template<class Num_T> 00104 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00105 Mat<Num_T> &out); 00107 template<class Num_T> 00108 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00109 const Mat<Num_T> &m3, Mat<Num_T> &out); 00111 template<class Num_T> 00112 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00113 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 00114 Mat<Num_T> &out); 00116 template<class Num_T> 00117 void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2); 00119 template<class Num_T> 00120 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00121 00123 template<class Num_T> 00124 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t); 00125 00127 template<class Num_T> 00128 Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00130 template<class Num_T> 00131 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00132 Mat<Num_T> &out); 00134 template<class Num_T> 00135 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00136 00137 // ------------------------------------------------------------------------------------- 00138 // Declaration of Mat 00139 // ------------------------------------------------------------------------------------- 00140 00206 template<class Num_T> 00207 class Mat { 00208 public: 00210 typedef Num_T value_type; 00211 00213 explicit Mat(const Factory &f = DEFAULT_FACTORY); 00215 Mat(int rows, int cols, const Factory &f = DEFAULT_FACTORY); 00217 Mat(const Mat<Num_T> &m); 00219 Mat(const Mat<Num_T> &m, const Factory &f); 00221 Mat(const Vec<Num_T> &v, const Factory &f = DEFAULT_FACTORY); 00223 Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY); 00225 Mat(const char *str, const Factory &f = DEFAULT_FACTORY); 00233 Mat(const Num_T *c_array, int rows, int cols, bool row_major = true, 00234 const Factory &f = DEFAULT_FACTORY); 00235 00237 ~Mat(); 00238 00240 int cols() const { return no_cols; } 00242 int rows() const { return no_rows; } 00244 int size() const { return datasize; } 00246 void set_size(int rows, int cols, bool copy = false); 00248 void zeros(); 00250 void clear() { zeros(); } 00252 void ones(); 00254 void set(const char *str); 00256 void set(const std::string &str); 00257 00259 const Num_T &operator()(int r, int c) const; 00261 Num_T &operator()(int r, int c); 00263 const Num_T &operator()(int i) const; 00265 Num_T &operator()(int i); 00267 const Num_T &get(int r, int c) const; 00269 void set(int r, int c, Num_T t); 00270 00276 Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const; 00282 Mat<Num_T> get(int r1, int r2, int c1, int c2) const; 00283 00285 Vec<Num_T> get_row(int r) const; 00287 Mat<Num_T> get_rows(int r1, int r2) const; 00289 Mat<Num_T> get_rows(const Vec<int> &indexlist) const; 00291 Vec<Num_T> get_col(int c) const; 00293 Mat<Num_T> get_cols(int c1, int c2) const; 00295 Mat<Num_T> get_cols(const Vec<int> &indexlist) const; 00297 void set_row(int r, const Vec<Num_T> &v); 00299 void set_col(int c, const Vec<Num_T> &v); 00301 void set_rows(int r, const Mat<Num_T> &m); 00303 void set_cols(int c, const Mat<Num_T> &m); 00305 void copy_row(int to, int from); 00307 void copy_col(int to, int from); 00309 void swap_rows(int r1, int r2); 00311 void swap_cols(int c1, int c2); 00312 00314 void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m); 00316 void set_submatrix(int r, int c, const Mat<Num_T> &m); 00318 void set_submatrix(int r1, int r2, int c1, int c2, Num_T t); 00319 00321 void del_row(int r); 00323 void del_rows(int r1, int r2); 00325 void del_col(int c); 00327 void del_cols(int c1, int c2); 00329 void ins_row(int r, const Vec<Num_T> &v); 00331 void ins_col(int c, const Vec<Num_T> &v); 00333 void append_row(const Vec<Num_T> &v); 00335 void append_col(const Vec<Num_T> &v); 00336 00338 Mat<Num_T> transpose() const; 00340 Mat<Num_T> T() const { return this->transpose(); } 00342 Mat<Num_T> hermitian_transpose() const; 00344 Mat<Num_T> H() const { return this->hermitian_transpose(); } 00345 00347 friend Mat<Num_T> concat_horizontal<>(const Mat<Num_T> &m1, 00348 const Mat<Num_T> &m2); 00350 friend Mat<Num_T> concat_vertical<>(const Mat<Num_T> &m1, 00351 const Mat<Num_T> &m2); 00352 00354 Mat<Num_T>& operator=(Num_T t); 00356 Mat<Num_T>& operator=(const Mat<Num_T> &m); 00358 Mat<Num_T>& operator=(const Vec<Num_T> &v); 00360 Mat<Num_T>& operator=(const char *str); 00361 00363 Mat<Num_T>& operator+=(const Mat<Num_T> &m); 00365 Mat<Num_T>& operator+=(Num_T t); 00367 friend Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00369 friend Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t); 00371 friend Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m); 00372 00374 Mat<Num_T>& operator-=(const Mat<Num_T> &m); 00376 Mat<Num_T>& operator-=(Num_T t); 00378 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00380 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t); 00382 friend Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m); 00384 friend Mat<Num_T> operator-<>(const Mat<Num_T> &m); 00385 00387 Mat<Num_T>& operator*=(const Mat<Num_T> &m); 00389 Mat<Num_T>& operator*=(Num_T t); 00391 friend Mat<Num_T> operator*<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00393 friend Vec<Num_T> operator*<>(const Mat<Num_T> &m, const Vec<Num_T> &v); 00405 friend Mat<Num_T> operator*<>(const Vec<Num_T> &v, const Mat<Num_T> &m); 00407 friend Mat<Num_T> operator*<>(const Mat<Num_T> &m, Num_T t); 00409 friend Mat<Num_T> operator*<>(Num_T t, const Mat<Num_T> &m); 00410 00412 friend Mat<Num_T> elem_mult<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00414 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00415 Mat<Num_T> &out); 00417 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00418 const Mat<Num_T> &m3, Mat<Num_T> &out); 00420 friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00421 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 00422 Mat<Num_T> &out); 00424 friend void elem_mult_inplace<>(const Mat<Num_T> &m1, Mat<Num_T> &m2); 00426 friend Num_T elem_mult_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00427 00429 Mat<Num_T>& operator/=(Num_T t); 00431 friend Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t); 00433 Mat<Num_T>& operator/=(const Mat<Num_T> &m); 00434 00436 friend Mat<Num_T> elem_div<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00438 friend void elem_div_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 00439 Mat<Num_T> &out); 00441 friend Num_T elem_div_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00442 00444 bool operator==(const Mat<Num_T> &m) const; 00446 bool operator!=(const Mat<Num_T> &m) const; 00447 00449 Num_T &_elem(int r, int c) { return data[r+c*no_rows]; } 00451 const Num_T &_elem(int r, int c) const { return data[r+c*no_rows]; } 00453 Num_T &_elem(int i) { return data[i]; } 00455 const Num_T &_elem(int i) const { return data[i]; } 00456 00458 Num_T *_data() { return data; } 00460 const Num_T *_data() const { return data; } 00462 int _datasize() const { return datasize; } 00463 00464 protected: 00466 void alloc(int rows, int cols); 00468 void free(); 00469 00472 int datasize, no_rows, no_cols; 00474 00475 Num_T *data; 00477 const Factory &factory; 00478 00479 private: 00481 bool in_range(int r, int c) const 00482 { 00483 return ((r >=0) && (r < no_rows) && (c >= 0) && (c < no_cols)); 00484 } 00486 bool row_in_range(int r) const { return ((r >= 0) && (r < no_rows)); } 00488 bool col_in_range(int c) const { return ((c >= 0) && (c < no_cols)); } 00490 bool in_range(int i) const { return ((i >= 0) && (i < datasize)); } 00491 }; 00492 00493 // ------------------------------------------------------------------------------------- 00494 // Type definitions of mat, cmat, imat, smat, and bmat 00495 // ------------------------------------------------------------------------------------- 00496 00501 typedef Mat<double> mat; 00502 00507 typedef Mat<std::complex<double> > cmat; 00508 00513 typedef Mat<int> imat; 00514 00519 typedef Mat<short int> smat; 00520 00527 typedef Mat<bin> bmat; 00528 00529 } //namespace itpp 00530 00531 00532 #include <itpp/base/vec.h> 00533 00534 namespace itpp { 00535 00536 // ---------------------------------------------------------------------- 00537 // Declaration of input and output streams for Mat 00538 // ---------------------------------------------------------------------- 00539 00544 template <class Num_T> 00545 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m); 00546 00558 template <class Num_T> 00559 std::istream &operator>>(std::istream &is, Mat<Num_T> &m); 00560 00561 // ---------------------------------------------------------------------- 00562 // Implementation of templated Mat members and friends 00563 // ---------------------------------------------------------------------- 00564 00565 template<class Num_T> inline 00566 void Mat<Num_T>::alloc(int rows, int cols) 00567 { 00568 if ((rows > 0) && (cols > 0)) { 00569 datasize = rows * cols; 00570 no_rows = rows; 00571 no_cols = cols; 00572 create_elements(data, datasize, factory); 00573 } 00574 else { 00575 data = 0; 00576 datasize = 0; 00577 no_rows = 0; 00578 no_cols = 0; 00579 } 00580 } 00581 00582 template<class Num_T> inline 00583 void Mat<Num_T>::free() 00584 { 00585 destroy_elements(data, datasize); 00586 datasize = 0; 00587 no_rows = 0; 00588 no_cols = 0; 00589 } 00590 00591 00592 template<class Num_T> inline 00593 Mat<Num_T>::Mat(const Factory &f) : 00594 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) {} 00595 00596 template<class Num_T> inline 00597 Mat<Num_T>::Mat(int rows, int cols, const Factory &f) : 00598 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00599 { 00600 it_assert_debug((rows >= 0) && (cols >= 0), "Mat<>::Mat(): Wrong size"); 00601 alloc(rows, cols); 00602 } 00603 00604 template<class Num_T> inline 00605 Mat<Num_T>::Mat(const Mat<Num_T> &m) : 00606 datasize(0), no_rows(0), no_cols(0), data(0), factory(m.factory) 00607 { 00608 alloc(m.no_rows, m.no_cols); 00609 copy_vector(m.datasize, m.data, data); 00610 } 00611 00612 template<class Num_T> inline 00613 Mat<Num_T>::Mat(const Mat<Num_T> &m, const Factory &f) : 00614 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00615 { 00616 alloc(m.no_rows, m.no_cols); 00617 copy_vector(m.datasize, m.data, data); 00618 } 00619 00620 template<class Num_T> inline 00621 Mat<Num_T>::Mat(const Vec<Num_T> &v, const Factory &f) : 00622 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00623 { 00624 int size = v.size(); 00625 alloc(size, 1); 00626 copy_vector(size, v._data(), data); 00627 } 00628 00629 template<class Num_T> inline 00630 Mat<Num_T>::Mat(const std::string &str, const Factory &f) : 00631 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00632 { 00633 set(str); 00634 } 00635 00636 template<class Num_T> inline 00637 Mat<Num_T>::Mat(const char *str, const Factory &f) : 00638 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00639 { 00640 set(str); 00641 } 00642 00643 template<class Num_T> 00644 Mat<Num_T>::Mat(const Num_T *c_array, int rows, int cols, bool row_major, 00645 const Factory &f): 00646 datasize(0), no_rows(0), no_cols(0), data(0), factory(f) 00647 { 00648 alloc(rows, cols); 00649 if (!row_major) 00650 copy_vector(datasize, c_array, data); 00651 else 00652 for (int i=0; i<rows; i++) 00653 for (int j=0; j<cols; j++) 00654 data[i+j*no_rows] = c_array[i*no_cols+j]; 00655 } 00656 00657 template<class Num_T> inline 00658 Mat<Num_T>::~Mat() 00659 { 00660 free(); 00661 } 00662 00663 00664 template<class Num_T> 00665 void Mat<Num_T>::set_size(int rows, int cols, bool copy) 00666 { 00667 it_assert_debug((rows >= 0) && (cols >= 0), 00668 "Mat<>::set_size(): Wrong size"); 00669 // check if we have to resize the current matrix 00670 if ((no_rows == rows) && (no_cols == cols)) 00671 return; 00672 // conditionally copy previous matrix content 00673 if (copy) { 00674 // create a temporary pointer to the allocated data 00675 Num_T* tmp = data; 00676 // store the current number of elements and number of rows 00677 int old_datasize = datasize; 00678 int old_rows = no_rows; 00679 // check the boundaries of the copied data 00680 int min_r = (no_rows < rows) ? no_rows : rows; 00681 int min_c = (no_cols < cols) ? no_cols : cols; 00682 // allocate new memory 00683 alloc(rows, cols); 00684 // copy the previous data into the allocated memory 00685 for (int i = 0; i < min_c; ++i) { 00686 copy_vector(min_r, &tmp[i*old_rows], &data[i*no_rows]); 00687 } 00688 // fill-in the rest of matrix with zeros 00689 for (int i = min_r; i < rows; ++i) 00690 for (int j = 0; j < cols; ++j) 00691 data[i+j*rows] = Num_T(0); 00692 for (int j = min_c; j < cols; ++j) 00693 for (int i = 0; i < min_r; ++i) 00694 data[i+j*rows] = Num_T(0); 00695 // delete old elements 00696 destroy_elements(tmp, old_datasize); 00697 } 00698 // if possible, reuse the allocated memory 00699 else if (datasize == rows * cols) { 00700 no_rows = rows; 00701 no_cols = cols; 00702 } 00703 // finally release old memory and allocate a new one 00704 else { 00705 free(); 00706 alloc(rows, cols); 00707 } 00708 } 00709 00710 template<class Num_T> inline 00711 void Mat<Num_T>::zeros() 00712 { 00713 for(int i=0; i<datasize; i++) 00714 data[i] = Num_T(0); 00715 } 00716 00717 template<class Num_T> inline 00718 void Mat<Num_T>::ones() 00719 { 00720 for(int i=0; i<datasize; i++) 00721 data[i] = Num_T(1); 00722 } 00723 00724 template<class Num_T> inline 00725 const Num_T& Mat<Num_T>::operator()(int r, int c) const 00726 { 00727 it_assert_debug(in_range(r, c), 00728 "Mat<>::operator(): Indexing out of range"); 00729 return data[r+c*no_rows]; 00730 } 00731 00732 template<class Num_T> inline 00733 Num_T& Mat<Num_T>::operator()(int r, int c) 00734 { 00735 it_assert_debug(in_range(r, c), 00736 "Mat<>::operator(): Indexing out of range"); 00737 return data[r+c*no_rows]; 00738 } 00739 00740 template<class Num_T> inline 00741 Num_T& Mat<Num_T>::operator()(int i) 00742 { 00743 it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range"); 00744 return data[i]; 00745 } 00746 00747 template<class Num_T> inline 00748 const Num_T& Mat<Num_T>::operator()(int i) const 00749 { 00750 it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range"); 00751 return data[i]; 00752 } 00753 00754 template<class Num_T> inline 00755 const Num_T& Mat<Num_T>::get(int r, int c) const 00756 { 00757 it_assert_debug(in_range(r, c), "Mat<>::get(): Indexing out of range"); 00758 return data[r+c*no_rows]; 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> 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 r2, int c1, int c2, 01028 const Mat<Num_T> &m) 01029 { 01030 01031 if (r1 == -1) r1 = no_rows-1; 01032 if (r2 == -1) r2 = no_rows-1; 01033 if (c1 == -1) c1 = no_cols-1; 01034 if (c2 == -1) c2 = no_cols-1; 01035 01036 it_assert_debug(r1>=0 && r2>=0 && r1<no_rows && r2<no_rows && 01037 c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, "Mat<Num_T>::set_submatrix(): index out of range"); 01038 01039 it_assert_debug(r2>=r1 && c2>=c1, "Mat<Num_T>::set_submatrix: r2<r1 or c2<c1"); 01040 it_assert_debug(m.no_rows == r2-r1+1 && m.no_cols == c2-c1+1, "Mat<Num_T>::set_submatrix(): sizes don't match"); 01041 01042 for (int i=0; i<m.no_cols; i++) 01043 copy_vector(m.no_rows, m.data+i*m.no_rows, data+(c1+i)*no_rows+r1); 01044 } 01045 01046 01047 01048 template<class Num_T> inline 01049 void Mat<Num_T>::set_submatrix(int r, int c, const Mat<Num_T> &m) 01050 { 01051 it_assert_debug((r >= 0) && (r + m.no_rows <= no_rows) && 01052 (c >= 0) && (c + m.no_cols <= no_cols), 01053 "Mat<>::set_submatrix(): Indexing out of range " 01054 "or wrong input matrix"); 01055 for (int i=0; i<m.no_cols; i++) 01056 copy_vector(m.no_rows, m.data+i*m.no_rows, data+(c+i)*no_rows+r); 01057 } 01058 01059 01060 01061 template<class Num_T> inline 01062 void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, Num_T t) 01063 { 01064 01065 if (r1 == -1) r1 = no_rows-1; 01066 if (r2 == -1) r2 = no_rows-1; 01067 if (c1 == -1) c1 = no_cols-1; 01068 if (c2 == -1) c2 = no_cols-1; 01069 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) && 01070 (c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 01071 "Mat<>::set_submatrix(): Wrong indexing"); 01072 01073 int i, j, pos, rows = r2-r1+1; 01074 01075 for (i=c1; i<=c2; i++) { 01076 pos = i*no_rows+r1; 01077 for (j=0; j<rows; j++) { 01078 data[pos++] = t; 01079 } 01080 } 01081 } 01082 01083 template<class Num_T> 01084 void Mat<Num_T>::del_row(int r) 01085 { 01086 it_assert_debug(row_in_range(r), "Mat<>::del_row(): Index out of range"); 01087 Mat<Num_T> Temp(*this); 01088 set_size(no_rows-1, no_cols, false); 01089 for (int i=0 ; i < r ; i++) { 01090 copy_vector(no_cols, &Temp.data[i], no_rows+1, &data[i], no_rows); 01091 } 01092 for (int i=r ; i < no_rows ; i++) { 01093 copy_vector(no_cols, &Temp.data[i+1], no_rows+1, &data[i], no_rows); 01094 } 01095 01096 } 01097 01098 template<class Num_T> 01099 void Mat<Num_T>::del_rows(int r1, int r2) 01100 { 01101 it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows), 01102 "Mat<>::del_rows(): Indexing out of range"); 01103 Mat<Num_T> Temp(*this); 01104 int no_del_rows = r2-r1+1; 01105 set_size(no_rows-no_del_rows, no_cols, false); 01106 for (int i = 0; i < r1 ; ++i) { 01107 copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows); 01108 } 01109 for (int i = r2+1; i < Temp.no_rows; ++i) { 01110 copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i-no_del_rows], 01111 no_rows); 01112 } 01113 } 01114 01115 template<class Num_T> 01116 void Mat<Num_T>::del_col(int c) 01117 { 01118 it_assert_debug(col_in_range(c), "Mat<>::del_col(): Index out of range"); 01119 Mat<Num_T> Temp(*this); 01120 01121 set_size(no_rows, no_cols-1, false); 01122 copy_vector(c*no_rows, Temp.data, data); 01123 copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]); 01124 } 01125 01126 template<class Num_T> 01127 void Mat<Num_T>::del_cols(int c1, int c2) 01128 { 01129 it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols), 01130 "Mat<>::del_cols(): Indexing out of range"); 01131 Mat<Num_T> Temp(*this); 01132 int n_deleted_cols = c2-c1+1; 01133 set_size(no_rows, no_cols-n_deleted_cols, false); 01134 copy_vector(c1*no_rows, Temp.data, data); 01135 copy_vector((no_cols-c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]); 01136 } 01137 01138 template<class Num_T> 01139 void Mat<Num_T>::ins_row(int r, const Vec<Num_T> &v) 01140 { 01141 it_assert_debug((r >= 0) && (r <= no_rows), 01142 "Mat<>::ins_row(): Index out of range"); 01143 it_assert_debug((v.size() == no_cols) || (no_rows == 0), 01144 "Mat<>::ins_row(): Wrong size of the input vector"); 01145 01146 if (no_cols==0) { 01147 no_cols = v.size(); 01148 } 01149 01150 Mat<Num_T> Temp(*this); 01151 set_size(no_rows+1, no_cols, false); 01152 01153 for (int i=0 ; i < r ; i++) { 01154 copy_vector(no_cols, &Temp.data[i], no_rows-1, &data[i], no_rows); 01155 } 01156 copy_vector(no_cols, v._data(), 1, &data[r], no_rows); 01157 for (int i=r+1 ; i < no_rows ; i++) { 01158 copy_vector(no_cols, &Temp.data[i-1], no_rows-1, &data[i], no_rows); 01159 } 01160 } 01161 01162 template<class Num_T> 01163 void Mat<Num_T>::ins_col(int c, const Vec<Num_T> &v) 01164 { 01165 it_assert_debug((c >= 0) && (c <= no_cols), 01166 "Mat<>::ins_col(): Index out of range"); 01167 it_assert_debug((v.size() == no_rows) || (no_cols == 0), 01168 "Mat<>::ins_col(): Wrong size of the input vector"); 01169 01170 if (no_rows==0) { 01171 no_rows = v.size(); 01172 } 01173 01174 Mat<Num_T> Temp(*this); 01175 set_size(no_rows, no_cols+1, false); 01176 01177 copy_vector(c*no_rows, Temp.data, data); 01178 copy_vector(no_rows, v._data(), &data[c*no_rows]); 01179 copy_vector((no_cols-c-1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]); 01180 } 01181 01182 template<class Num_T> inline 01183 void Mat<Num_T>::append_row(const Vec<Num_T> &v) 01184 { 01185 ins_row(no_rows, v); 01186 } 01187 01188 template<class Num_T> inline 01189 void Mat<Num_T>::append_col(const Vec<Num_T> &v) 01190 { 01191 ins_col(no_cols, v); 01192 } 01193 01194 template<class Num_T> 01195 Mat<Num_T> Mat<Num_T>::transpose() const 01196 { 01197 Mat<Num_T> temp(no_cols, no_rows); 01198 for (int i = 0; i < no_rows; ++i) { 01199 copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1); 01200 } 01201 return temp; 01202 } 01203 01205 template<> 01206 cmat cmat::hermitian_transpose() const; 01208 01209 template<class Num_T> 01210 Mat<Num_T> Mat<Num_T>::hermitian_transpose() const 01211 { 01212 Mat<Num_T> temp(no_cols, no_rows); 01213 for (int i = 0; i < no_rows; ++i) { 01214 copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1); 01215 } 01216 return temp; 01217 } 01218 01219 template<class Num_T> 01220 Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01221 { 01222 it_assert_debug(m1.no_rows == m2.no_rows, 01223 "Mat<>::concat_horizontal(): Wrong sizes"); 01224 int no_rows = m1.no_rows; 01225 Mat<Num_T> temp(no_rows, m1.no_cols + m2.no_cols); 01226 for (int i = 0; i < m1.no_cols; ++i) { 01227 copy_vector(no_rows, &m1.data[i * no_rows], &temp.data[i * no_rows]); 01228 } 01229 for (int i = 0; i < m2.no_cols; ++i) { 01230 copy_vector(no_rows, &m2.data[i * no_rows], &temp.data[(m1.no_cols + i) 01231 * no_rows]); 01232 } 01233 return temp; 01234 } 01235 01236 template<class Num_T> 01237 Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01238 { 01239 it_assert_debug(m1.no_cols == m2.no_cols, 01240 "Mat<>::concat_vertical(): Wrong sizes"); 01241 int no_cols = m1.no_cols; 01242 Mat<Num_T> temp(m1.no_rows + m2.no_rows, no_cols); 01243 for (int i = 0; i < no_cols; ++i) { 01244 copy_vector(m1.no_rows, &m1.data[i * m1.no_rows], 01245 &temp.data[i * temp.no_rows]); 01246 copy_vector(m2.no_rows, &m2.data[i * m2.no_rows], 01247 &temp.data[i * temp.no_rows + m1.no_rows]); 01248 } 01249 return temp; 01250 } 01251 01252 template<class Num_T> inline 01253 Mat<Num_T>& Mat<Num_T>::operator=(Num_T t) 01254 { 01255 for (int i=0; i<datasize; i++) 01256 data[i] = t; 01257 return *this; 01258 } 01259 01260 template<class Num_T> inline 01261 Mat<Num_T>& Mat<Num_T>::operator=(const Mat<Num_T> &m) 01262 { 01263 if (this != &m) { 01264 set_size(m.no_rows,m.no_cols, false); 01265 if (m.datasize != 0) 01266 copy_vector(m.datasize, m.data, data); 01267 } 01268 return *this; 01269 } 01270 01271 template<class Num_T> inline 01272 Mat<Num_T>& Mat<Num_T>::operator=(const Vec<Num_T> &v) 01273 { 01274 it_assert_debug(((no_rows == 1) && (no_cols == v.size())) 01275 || ((no_cols == 1) && (no_rows == v.size())), 01276 "Mat<>::operator=(): Wrong size of the input vector"); 01277 set_size(v.size(), 1, false); 01278 copy_vector(v.size(), v._data(), data); 01279 return *this; 01280 } 01281 01282 template<class Num_T> inline 01283 Mat<Num_T>& Mat<Num_T>::operator=(const char *str) 01284 { 01285 set(str); 01286 return *this; 01287 } 01288 01289 //-------------------- Templated friend functions -------------------------- 01290 01291 template<class Num_T> 01292 Mat<Num_T>& Mat<Num_T>::operator+=(const Mat<Num_T> &m) 01293 { 01294 if (datasize == 0) 01295 operator=(m); 01296 else { 01297 int i, j, m_pos=0, pos=0; 01298 it_assert_debug(m.no_rows==no_rows && m.no_cols==no_cols,"Mat<Num_T>::operator+=: wrong sizes"); 01299 for (i=0; i<no_cols; i++) { 01300 for (j=0; j<no_rows; j++) 01301 data[pos+j] += m.data[m_pos+j]; 01302 pos += no_rows; 01303 m_pos += m.no_rows; 01304 } 01305 } 01306 return *this; 01307 } 01308 01309 template<class Num_T> inline 01310 Mat<Num_T>& Mat<Num_T>::operator+=(Num_T t) 01311 { 01312 for (int i=0; i<datasize; i++) 01313 data[i] += t; 01314 return *this; 01315 } 01316 01317 template<class Num_T> 01318 Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01319 { 01320 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01321 int i, j, m1_pos=0, m2_pos=0, r_pos=0; 01322 01323 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01324 "Mat<>::operator+(): Wrong sizes"); 01325 01326 for (i=0; i<r.no_cols; i++) { 01327 for (j=0; j<r.no_rows; j++) 01328 r.data[r_pos+j] = m1.data[m1_pos+j] + m2.data[m2_pos+j]; 01329 // next column 01330 m1_pos += m1.no_rows; 01331 m2_pos += m2.no_rows; 01332 r_pos += r.no_rows; 01333 } 01334 01335 return r; 01336 } 01337 01338 01339 template<class Num_T> 01340 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t) 01341 { 01342 Mat<Num_T> r(m.no_rows, m.no_cols); 01343 01344 for (int i=0; i<r.datasize; i++) 01345 r.data[i] = m.data[i] + t; 01346 01347 return r; 01348 } 01349 01350 template<class Num_T> 01351 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m) 01352 { 01353 Mat<Num_T> r(m.no_rows, m.no_cols); 01354 01355 for (int i=0; i<r.datasize; i++) 01356 r.data[i] = t + m.data[i]; 01357 01358 return r; 01359 } 01360 01361 template<class Num_T> 01362 Mat<Num_T>& Mat<Num_T>::operator-=(const Mat<Num_T> &m) 01363 { 01364 int i,j, m_pos=0, pos=0; 01365 01366 if (datasize == 0) { 01367 set_size(m.no_rows, m.no_cols, false); 01368 for (i=0; i<no_cols; i++) { 01369 for (j=0; j<no_rows; j++) 01370 data[pos+j] = -m.data[m_pos+j]; 01371 // next column 01372 m_pos += m.no_rows; 01373 pos += no_rows; 01374 } 01375 } 01376 else { 01377 it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols), 01378 "Mat<>::operator-=(): Wrong sizes"); 01379 for (i=0; i<no_cols; i++) { 01380 for (j=0; j<no_rows; j++) 01381 data[pos+j] -= m.data[m_pos+j]; 01382 // next column 01383 m_pos += m.no_rows; 01384 pos += no_rows; 01385 } 01386 } 01387 return *this; 01388 } 01389 01390 template<class Num_T> 01391 Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01392 { 01393 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01394 int i, j, m1_pos=0, m2_pos=0, r_pos=0; 01395 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01396 "Mat<>::operator-(): Wrong sizes"); 01397 01398 for (i=0; i<r.no_cols; i++) { 01399 for (j=0; j<r.no_rows; j++) 01400 r.data[r_pos+j] = m1.data[m1_pos+j] - m2.data[m2_pos+j]; 01401 // next column 01402 m1_pos += m1.no_rows; 01403 m2_pos += m2.no_rows; 01404 r_pos += r.no_rows; 01405 } 01406 01407 return r; 01408 } 01409 01410 template<class Num_T> inline 01411 Mat<Num_T>& Mat<Num_T>::operator-=(Num_T t) 01412 { 01413 for (int i=0; i<datasize; i++) 01414 data[i] -= t; 01415 return *this; 01416 } 01417 01418 template<class Num_T> 01419 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t) 01420 { 01421 Mat<Num_T> r(m.no_rows, m.no_cols); 01422 int i, j, m_pos=0, r_pos=0; 01423 01424 for (i=0; i<r.no_cols; i++) { 01425 for (j=0; j<r.no_rows; j++) 01426 r.data[r_pos+j] = m.data[m_pos+j] - t; 01427 // next column 01428 m_pos += m.no_rows; 01429 r_pos += r.no_rows; 01430 } 01431 01432 return r; 01433 } 01434 01435 template<class Num_T> 01436 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m) 01437 { 01438 Mat<Num_T> r(m.no_rows, m.no_cols); 01439 int i, j, m_pos=0, r_pos=0; 01440 01441 for (i=0; i<r.no_cols; i++) { 01442 for (j=0; j<r.no_rows; j++) 01443 r.data[r_pos+j] = t - m.data[m_pos+j]; 01444 // next column 01445 m_pos += m.no_rows; 01446 r_pos += r.no_rows; 01447 } 01448 01449 return r; 01450 } 01451 01452 template<class Num_T> 01453 Mat<Num_T> operator-(const Mat<Num_T> &m) 01454 { 01455 Mat<Num_T> r(m.no_rows, m.no_cols); 01456 int i, j, m_pos=0, r_pos=0; 01457 01458 for (i=0; i<r.no_cols; i++) { 01459 for (j=0; j<r.no_rows; j++) 01460 r.data[r_pos+j] = -m.data[m_pos+j]; 01461 // next column 01462 m_pos += m.no_rows; 01463 r_pos += r.no_rows; 01464 } 01465 01466 return r; 01467 } 01468 01469 #if defined(HAVE_BLAS) 01470 template<> mat& mat::operator*=(const mat &m); 01471 template<> cmat& cmat::operator*=(const cmat &m); 01472 #endif 01473 01474 template<class Num_T> 01475 Mat<Num_T>& Mat<Num_T>::operator*=(const Mat<Num_T> &m) 01476 { 01477 it_assert_debug(no_cols == m.no_rows,"Mat<>::operator*=(): Wrong sizes"); 01478 Mat<Num_T> r(no_rows, m.no_cols); 01479 01480 Num_T tmp; 01481 01482 int i,j,k, r_pos=0, pos=0, m_pos=0; 01483 01484 for (i=0; i<r.no_cols; i++) { 01485 for (j=0; j<r.no_rows; j++) { 01486 tmp = Num_T(0); 01487 pos = 0; 01488 for (k=0; k<no_cols; k++) { 01489 tmp += data[pos+j] * m.data[m_pos+k]; 01490 pos += no_rows; 01491 } 01492 r.data[r_pos+j] = tmp; 01493 } 01494 r_pos += r.no_rows; 01495 m_pos += m.no_rows; 01496 } 01497 operator=(r); // time consuming 01498 return *this; 01499 } 01500 01501 template<class Num_T> inline 01502 Mat<Num_T>& Mat<Num_T>::operator*=(Num_T t) 01503 { 01504 scal_vector(datasize, t, data); 01505 return *this; 01506 } 01507 01508 #if defined(HAVE_BLAS) 01509 template<> mat operator*(const mat &m1, const mat &m2); 01510 template<> cmat operator*(const cmat &m1, const cmat &m2); 01511 #endif 01512 01513 01514 template<class Num_T> 01515 Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01516 { 01517 it_assert_debug(m1.no_cols == m2.no_rows, 01518 "Mat<>::operator*(): Wrong sizes"); 01519 Mat<Num_T> r(m1.no_rows, m2.no_cols); 01520 01521 Num_T tmp; 01522 int i, j, k; 01523 Num_T *tr=r.data, *t1, *t2=m2.data; 01524 01525 for (i=0; i<r.no_cols; i++) { 01526 for (j=0; j<r.no_rows; j++) { 01527 tmp = Num_T(0); t1 = m1.data+j; 01528 for (k=m1.no_cols; k>0; k--) { 01529 tmp += *(t1) * *(t2++); 01530 t1 += m1.no_rows; 01531 } 01532 *(tr++) = tmp; t2 -= m2.no_rows; 01533 } 01534 t2 += m2.no_rows; 01535 } 01536 01537 return r; 01538 } 01539 01540 #if defined(HAVE_BLAS) 01541 template<> vec operator*(const mat &m, const vec &v); 01542 template<> cvec operator*(const cmat &m, const cvec &v); 01543 #endif 01544 01545 template<class Num_T> 01546 Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v) 01547 { 01548 it_assert_debug(m.no_cols == v.size(), 01549 "Mat<>::operator*(): Wrong sizes"); 01550 Vec<Num_T> r(m.no_rows); 01551 int i, k, m_pos; 01552 01553 for (i=0; i<m.no_rows; i++) { 01554 r(i) = Num_T(0); 01555 m_pos = 0; 01556 for (k=0; k<m.no_cols; k++) { 01557 r(i) += m.data[m_pos+i] * v(k); 01558 m_pos += m.no_rows; 01559 } 01560 } 01561 01562 return r; 01563 } 01564 01565 template<class Num_T> 01566 Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m) 01567 { 01568 it_assert((m.no_rows == 1),"Mat<Num_T>::operator*(): wrong sizes"); 01569 it_warning("Mat<Num_T>::operator*(v, m): This operator is deprecated. " 01570 "Please use outer_product(v, m.get_row(0)) instead."); 01571 return outer_product(v, m.get_row(0)); 01572 } 01573 01574 template<class Num_T> 01575 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t) 01576 { 01577 Mat<Num_T> r(m.no_rows, m.no_cols); 01578 01579 for (int i=0; i<r.datasize; i++) 01580 r.data[i] = m.data[i] * t; 01581 01582 return r; 01583 } 01584 01585 template<class Num_T> inline 01586 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m) 01587 { 01588 return operator*(m, t); 01589 } 01590 01591 template<class Num_T> inline 01592 Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01593 { 01594 Mat<Num_T> out; 01595 elem_mult_out(m1,m2,out); 01596 return out; 01597 } 01598 01599 template<class Num_T> 01600 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01601 Mat<Num_T> &out) 01602 { 01603 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.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]; 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, Mat<Num_T> &out) 01613 { 01614 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows) 01615 && (m1.no_cols == m2.no_cols) && (m1.no_cols == m3.no_cols), 01616 "Mat<>::elem_mult_out(): Wrong sizes"); 01617 out.set_size(m1.no_rows, m1.no_cols); 01618 for(int i=0; i<out.datasize; i++) 01619 out.data[i] = m1.data[i] * m2.data[i] * m3.data[i]; 01620 } 01621 01622 template<class Num_T> 01623 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01624 const Mat<Num_T> &m3, const Mat<Num_T> &m4, 01625 Mat<Num_T> &out) 01626 { 01627 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows) 01628 && (m1.no_rows == m4.no_rows) && (m1.no_cols == m2.no_cols) 01629 && (m1.no_cols == m3.no_cols) && (m1.no_cols == m4.no_cols), 01630 "Mat<>::elem_mult_out(): Wrong sizes"); 01631 out.set_size(m1.no_rows, m1.no_cols); 01632 for(int i=0; i<out.datasize; i++) 01633 out.data[i] = m1.data[i] * m2.data[i] * m3.data[i] * m4.data[i]; 01634 } 01635 01636 template<class Num_T> inline 01637 void elem_mult_inplace(const Mat<Num_T> &m1, 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_inplace(): Wrong sizes"); 01641 for(int i=0; i<m2.datasize; i++) 01642 m2.data[i] *= m1.data[i]; 01643 } 01644 01645 template<class Num_T> inline 01646 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01647 { 01648 it_assert_debug( (m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01649 "Mat<>::elem_mult_sum(): Wrong sizes"); 01650 Num_T acc = 0; 01651 01652 for(int i=0; i<m1.datasize; i++) 01653 acc += m1.data[i] * m2.data[i]; 01654 01655 return acc; 01656 } 01657 01658 template<class Num_T> inline 01659 Mat<Num_T>& Mat<Num_T>::operator/=(Num_T t) 01660 { 01661 for (int i=0; i<datasize; i++) 01662 data[i] /= t; 01663 return *this; 01664 } 01665 01666 template<class Num_T> 01667 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t) 01668 { 01669 Mat<Num_T> r(m.no_rows, m.no_cols); 01670 01671 for (int i=0; i<r.datasize; i++) 01672 r.data[i] = m.data[i] / t; 01673 01674 return r; 01675 } 01676 01677 template<class Num_T> inline 01678 Mat<Num_T>& Mat<Num_T>::operator/=(const Mat<Num_T> &m) 01679 { 01680 it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols), 01681 "Mat<>::operator/=(): Wrong sizes"); 01682 for (int i=0; i<datasize; i++) 01683 data[i] /= m.data[i]; 01684 return *this; 01685 } 01686 01687 template<class Num_T> inline 01688 Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01689 { 01690 Mat<Num_T> out; 01691 elem_div_out(m1,m2,out); 01692 return out; 01693 } 01694 01695 template<class Num_T> 01696 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2, 01697 Mat<Num_T> &out) 01698 { 01699 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01700 "Mat<>::elem_div_out(): Wrong sizes"); 01701 01702 if( (out.no_rows != m1.no_rows) || (out.no_cols != m1.no_cols) ) 01703 out.set_size(m1.no_rows, m1.no_cols); 01704 01705 for(int i=0; i<out.datasize; i++) 01706 out.data[i] = m1.data[i] / m2.data[i]; 01707 } 01708 01709 template<class Num_T> inline 01710 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01711 { 01712 it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols), 01713 "Mat<>::elem_div_sum(): Wrong sizes"); 01714 Num_T acc = 0; 01715 01716 for(int i=0; i<m1.datasize; i++) 01717 acc += m1.data[i] / m2.data[i]; 01718 01719 return acc; 01720 } 01721 01722 template<class Num_T> 01723 bool Mat<Num_T>::operator==(const Mat<Num_T> &m) const 01724 { 01725 if (no_rows!=m.no_rows || no_cols != m.no_cols) return false; 01726 for (int i=0;i<datasize;i++) { 01727 if (data[i]!=m.data[i]) return false; 01728 } 01729 return true; 01730 } 01731 01732 template<class Num_T> 01733 bool Mat<Num_T>::operator!=(const Mat<Num_T> &m) const 01734 { 01735 if (no_rows != m.no_rows || no_cols != m.no_cols) return true; 01736 for (int i=0;i<datasize;i++) { 01737 if (data[i]!=m.data[i]) return true; 01738 } 01739 return false; 01740 } 01741 01742 template <class Num_T> 01743 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m) 01744 { 01745 int i; 01746 01747 switch (m.rows()) { 01748 case 0 : 01749 os << "[]"; 01750 break; 01751 case 1 : 01752 os << '[' << m.get_row(0) << ']'; 01753 break; 01754 default: 01755 os << '[' << m.get_row(0) << std::endl; 01756 for (i=1; i<m.rows()-1; i++) 01757 os << ' ' << m.get_row(i) << std::endl; 01758 os << ' ' << m.get_row(m.rows()-1) << ']'; 01759 } 01760 01761 return os; 01762 } 01763 01764 template <class Num_T> 01765 std::istream &operator>>(std::istream &is, Mat<Num_T> &m) 01766 { 01767 std::ostringstream buffer; 01768 bool started = false; 01769 bool finished = false; 01770 bool brackets = false; 01771 bool within_double_brackets = false; 01772 char c; 01773 01774 while (!finished) { 01775 if (is.eof()) { 01776 finished = true; 01777 } else { 01778 c = is.get(); 01779 01780 if (is.eof() || (c == '\n')) { 01781 if (brackets) { 01782 // Right bracket missing 01783 is.setstate(std::ios_base::failbit); 01784 finished = true; 01785 } else if (!((c == '\n') && !started)) { 01786 finished = true; 01787 } 01788 } else if ((c == ' ') || (c == '\t')) { 01789 if (started) { 01790 buffer << ' '; 01791 } 01792 } else if (c == '[') { 01793 if ((started && !brackets) || within_double_brackets) { 01794 // Unexpected left bracket 01795 is.setstate(std::ios_base::failbit); 01796 finished = true; 01797 } else if (!started) { 01798 started = true; 01799 brackets = true; 01800 } else { 01801 within_double_brackets = true; 01802 } 01803 } else if (c == ']') { 01804 if (!started || !brackets) { 01805 // Unexpected right bracket 01806 is.setstate(std::ios_base::failbit); 01807 finished = true; 01808 } else if (within_double_brackets) { 01809 within_double_brackets = false; 01810 buffer << ';'; 01811 } else { 01812 finished = true; 01813 } 01814 while (!is.eof() && (((c = is.peek()) == ' ') || (c == '\t'))) { 01815 is.get(); 01816 } 01817 if (!is.eof() && (c == '\n')) { 01818 is.get(); 01819 } 01820 } else { 01821 started = true; 01822 buffer << c; 01823 } 01824 } 01825 } 01826 01827 if (!started) { 01828 m.set_size(0, false); 01829 } else { 01830 m.set(buffer.str()); 01831 } 01832 01833 return is; 01834 } 01835 01837 01838 // --------------------------------------------------------------------- 01839 // Instantiations 01840 // --------------------------------------------------------------------- 01841 01842 #ifdef HAVE_EXTERN_TEMPLATE 01843 01844 // class instantiations 01845 01846 extern template class Mat<double>; 01847 extern template class Mat<std::complex<double> >; 01848 extern template class Mat<int>; 01849 extern template class Mat<short int>; 01850 extern template class Mat<bin>; 01851 01852 // addition operators 01853 01854 extern template mat operator+(const mat &m1, const mat &m2); 01855 extern template cmat operator+(const cmat &m1, const cmat &m2); 01856 extern template imat operator+(const imat &m1, const imat &m2); 01857 extern template smat operator+(const smat &m1, const smat &m2); 01858 extern template bmat operator+(const bmat &m1, const bmat &m2); 01859 01860 extern template mat operator+(const mat &m, double t); 01861 extern template cmat operator+(const cmat &m, std::complex<double> t); 01862 extern template imat operator+(const imat &m, int t); 01863 extern template smat operator+(const smat &m, short t); 01864 extern template bmat operator+(const bmat &m, bin t); 01865 01866 extern template mat operator+(double t, const mat &m); 01867 extern template cmat operator+(std::complex<double> t, const cmat &m); 01868 extern template imat operator+(int t, const imat &m); 01869 extern template smat operator+(short t, const smat &m); 01870 extern template bmat operator+(bin t, const bmat &m); 01871 01872 // subtraction operators 01873 01874 extern template mat operator-(const mat &m1, const mat &m2); 01875 extern template cmat operator-(const cmat &m1, const cmat &m2); 01876 extern template imat operator-(const imat &m1, const imat &m2); 01877 extern template smat operator-(const smat &m1, const smat &m2); 01878 extern template bmat operator-(const bmat &m1, const bmat &m2); 01879 01880 extern template mat operator-(const mat &m, double t); 01881 extern template cmat operator-(const cmat &m, std::complex<double> t); 01882 extern template imat operator-(const imat &m, int t); 01883 extern template smat operator-(const smat &m, short t); 01884 extern template bmat operator-(const bmat &m, bin t); 01885 01886 extern template mat operator-(double t, const mat &m); 01887 extern template cmat operator-(std::complex<double> t, const cmat &m); 01888 extern template imat operator-(int t, const imat &m); 01889 extern template smat operator-(short t, const smat &m); 01890 extern template bmat operator-(bin t, const bmat &m); 01891 01892 // unary minus 01893 01894 extern template mat operator-(const mat &m); 01895 extern template cmat operator-(const cmat &m); 01896 extern template imat operator-(const imat &m); 01897 extern template smat operator-(const smat &m); 01898 extern template bmat operator-(const bmat &m); 01899 01900 // multiplication operators 01901 01902 #if !defined(HAVE_BLAS) 01903 extern template mat operator*(const mat &m1, const mat &m2); 01904 extern template cmat operator*(const cmat &m1, const cmat &m2); 01905 #endif 01906 extern template imat operator*(const imat &m1, const imat &m2); 01907 extern template smat operator*(const smat &m1, const smat &m2); 01908 extern template bmat operator*(const bmat &m1, const bmat &m2); 01909 01910 #if !defined(HAVE_BLAS) 01911 extern template vec operator*(const mat &m, const vec &v); 01912 extern template cvec operator*(const cmat &m, const cvec &v); 01913 #endif 01914 extern template ivec operator*(const imat &m, const ivec &v); 01915 extern template svec operator*(const smat &m, const svec &v); 01916 extern template bvec operator*(const bmat &m, const bvec &v); 01917 01918 extern template mat operator*(const vec &v, const mat &m); 01919 extern template cmat operator*(const cvec &v, const cmat &m); 01920 extern template imat operator*(const ivec &v, const imat &m); 01921 extern template smat operator*(const svec &v, const smat &m); 01922 extern template bmat operator*(const bvec &v, const bmat &m); 01923 01924 extern template mat operator*(const mat &m, double t); 01925 extern template cmat operator*(const cmat &m, std::complex<double> t); 01926 extern template imat operator*(const imat &m, int t); 01927 extern template smat operator*(const smat &m, short t); 01928 extern template bmat operator*(const bmat &m, bin t); 01929 01930 extern template mat operator*(double t, const mat &m); 01931 extern template cmat operator*(std::complex<double> t, const cmat &m); 01932 extern template imat operator*(int t, const imat &m); 01933 extern template smat operator*(short t, const smat &m); 01934 extern template bmat operator*(bin t, const bmat &m); 01935 01936 // element-wise multiplication 01937 01938 extern template mat elem_mult(const mat &m1, const mat &m2); 01939 extern template cmat elem_mult(const cmat &m1, const cmat &m2); 01940 extern template imat elem_mult(const imat &m1, const imat &m2); 01941 extern template smat elem_mult(const smat &m1, const smat &m2); 01942 extern template bmat elem_mult(const bmat &m1, const bmat &m2); 01943 01944 extern template void elem_mult_out(const mat &m1, const mat &m2, mat &out); 01945 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01946 cmat &out); 01947 extern template void elem_mult_out(const imat &m1, const imat &m2, 01948 imat &out); 01949 extern template void elem_mult_out(const smat &m1, const smat &m2, 01950 smat &out); 01951 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01952 bmat &out); 01953 01954 extern template void elem_mult_out(const mat &m1, const mat &m2, 01955 const mat &m3, mat &out); 01956 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01957 const cmat &m3, cmat &out); 01958 extern template void elem_mult_out(const imat &m1, const imat &m2, 01959 const imat &m3, imat &out); 01960 extern template void elem_mult_out(const smat &m1, const smat &m2, 01961 const smat &m3, smat &out); 01962 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01963 const bmat &m3, bmat &out); 01964 01965 extern template void elem_mult_out(const mat &m1, const mat &m2, 01966 const mat &m3, const mat &m4, mat &out); 01967 extern template void elem_mult_out(const cmat &m1, const cmat &m2, 01968 const cmat &m3, const cmat &m4, 01969 cmat &out); 01970 extern template void elem_mult_out(const imat &m1, const imat &m2, 01971 const imat &m3, const imat &m4, 01972 imat &out); 01973 extern template void elem_mult_out(const smat &m1, const smat &m2, 01974 const smat &m3, const smat &m4, 01975 smat &out); 01976 extern template void elem_mult_out(const bmat &m1, const bmat &m2, 01977 const bmat &m3, const bmat &m4, 01978 bmat &out); 01979 01980 extern template void elem_mult_inplace(const mat &m1, mat &m2); 01981 extern template void elem_mult_inplace(const cmat &m1, cmat &m2); 01982 extern template void elem_mult_inplace(const imat &m1, imat &m2); 01983 extern template void elem_mult_inplace(const smat &m1, smat &m2); 01984 extern template void elem_mult_inplace(const bmat &m1, bmat &m2); 01985 01986 extern template double elem_mult_sum(const mat &m1, const mat &m2); 01987 extern template std::complex<double> elem_mult_sum(const cmat &m1, 01988 const cmat &m2); 01989 extern template int elem_mult_sum(const imat &m1, const imat &m2); 01990 extern template short elem_mult_sum(const smat &m1, const smat &m2); 01991 extern template bin elem_mult_sum(const bmat &m1, const bmat &m2); 01992 01993 // division operator 01994 01995 extern template mat operator/(const mat &m, double t); 01996 extern template cmat operator/(const cmat &m, std::complex<double> t); 01997 extern template imat operator/(const imat &m, int t); 01998 extern template smat operator/(const smat &m, short t); 01999 extern template bmat operator/(const bmat &m, bin t); 02000 02001 // element-wise division 02002 02003 extern template mat elem_div(const mat &m1, const mat &m2); 02004 extern template cmat elem_div(const cmat &m1, const cmat &m2); 02005 extern template imat elem_div(const imat &m1, const imat &m2); 02006 extern template smat elem_div(const smat &m1, const smat &m2); 02007 extern template bmat elem_div(const bmat &m1, const bmat &m2); 02008 02009 extern template void elem_div_out(const mat &m1, const mat &m2, mat &out); 02010 extern template void elem_div_out(const cmat &m1, const cmat &m2, cmat &out); 02011 extern template void elem_div_out(const imat &m1, const imat &m2, imat &out); 02012 extern template void elem_div_out(const smat &m1, const smat &m2, smat &out); 02013 extern template void elem_div_out(const bmat &m1, const bmat &m2, bmat &out); 02014 02015 extern template double elem_div_sum(const mat &m1, const mat &m2); 02016 extern template std::complex<double> elem_div_sum(const cmat &m1, 02017 const cmat &m2); 02018 extern template int elem_div_sum(const imat &m1, const imat &m2); 02019 extern template short elem_div_sum(const smat &m1, const smat &m2); 02020 extern template bin elem_div_sum(const bmat &m1, const bmat &m2); 02021 02022 // concatenation 02023 02024 extern template mat concat_horizontal(const mat &m1, const mat &m2); 02025 extern template cmat concat_horizontal(const cmat &m1, const cmat &m2); 02026 extern template imat concat_horizontal(const imat &m1, const imat &m2); 02027 extern template smat concat_horizontal(const smat &m1, const smat &m2); 02028 extern template bmat concat_horizontal(const bmat &m1, const bmat &m2); 02029 02030 extern template mat concat_vertical(const mat &m1, const mat &m2); 02031 extern template cmat concat_vertical(const cmat &m1, const cmat &m2); 02032 extern template imat concat_vertical(const imat &m1, const imat &m2); 02033 extern template smat concat_vertical(const smat &m1, const smat &m2); 02034 extern template bmat concat_vertical(const bmat &m1, const bmat &m2); 02035 02036 // I/O streams 02037 02038 extern template std::ostream &operator<<(std::ostream &os, const mat &m); 02039 extern template std::ostream &operator<<(std::ostream &os, const cmat &m); 02040 extern template std::ostream &operator<<(std::ostream &os, const imat &m); 02041 extern template std::ostream &operator<<(std::ostream &os, const smat &m); 02042 extern template std::ostream &operator<<(std::ostream &os, const bmat &m); 02043 02044 extern template std::istream &operator>>(std::istream &is, mat &m); 02045 extern template std::istream &operator>>(std::istream &is, cmat &m); 02046 extern template std::istream &operator>>(std::istream &is, imat &m); 02047 extern template std::istream &operator>>(std::istream &is, smat &m); 02048 extern template std::istream &operator>>(std::istream &is, bmat &m); 02049 02050 #endif // HAVE_EXTERN_TEMPLATE 02051 02053 02054 } // namespace itpp 02055 02056 #endif // #ifndef MAT_H
Generated on Sat Apr 19 10:41:55 2008 for IT++ by Doxygen 1.5.5