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