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