00001 00033 #ifndef MAT_H 00034 #define MAT_H 00035 00036 #ifndef _MSC_VER 00037 # include <itpp/config.h> 00038 #else 00039 # include <itpp/config_msvc.h> 00040 #endif 00041 00042 #include <string> 00043 #include <complex> 00044 00045 #include <itpp/itconfig.h> 00046 #include <itpp/base/itassert.h> 00047 #include <itpp/base/factory.h> 00048 00049 00050 namespace itpp { 00051 00052 // Declaration of Vec 00053 template<class Num_T> class Vec; 00054 // Declaration of Mat 00055 template<class Num_T> class Mat; 00056 // Declaration of bin 00057 class bin; 00058 00059 // ------------------------------------------------------------------------------------- 00060 // Declaration of Mat Friends 00061 // ------------------------------------------------------------------------------------- 00062 00064 template<class Num_T> const Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00066 template<class Num_T> const Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00067 00069 template<class Num_T> const Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00071 template<class Num_T> const Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t); 00073 template<class Num_T> const Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m); 00074 00076 template<class Num_T> const Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00078 template<class Num_T> const Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t); 00080 template<class Num_T> const Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m); 00082 template<class Num_T> const Mat<Num_T> operator-(const Mat<Num_T> &m); 00083 00085 template<class Num_T> const Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00087 template<class Num_T> const Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v); 00089 template<class Num_T> const Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m); 00091 template<class Num_T> const Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t); 00093 template<class Num_T> const Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m); 00095 template<class Num_T> const Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00096 00098 template<class Num_T> const Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t); 00100 template<class Num_T> const Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00101 00102 // ------------------------------------------------------------------------------------- 00103 // Declaration of Mat 00104 // ------------------------------------------------------------------------------------- 00105 00170 template<class Num_T> 00171 class Mat { 00172 public: 00174 typedef Num_T value_type; 00176 explicit Mat(const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); } 00178 Mat(int inrow, int incol, const Factory &f = DEFAULT_FACTORY) : factory(f) { 00179 init(); 00180 it_assert1(inrow>=0 && incol>=0, "The rows and columns must be >= 0"); 00181 alloc(inrow, incol); } 00183 Mat(const Mat<Num_T> &m); 00185 Mat(const Mat<Num_T> &m, const Factory &f); 00187 Mat(const Vec<Num_T> &invector, const Factory &f = DEFAULT_FACTORY); 00189 Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(str); } 00191 Mat(const char *str, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(str); } 00198 Mat(Num_T *c_array, int rows, int cols, bool RowMajor = true, const Factory &f = DEFAULT_FACTORY); 00199 00201 ~Mat() { free(); } 00202 00204 int cols() const { return no_cols; } 00206 int rows() const { return no_rows; } 00208 int size() const { return datasize; } 00210 void set_size(int inrow, int incol, bool copy=false); 00212 void zeros(); 00214 void clear() { zeros(); } 00216 void ones(); 00218 bool set(const char *values); 00220 bool set(const std::string &str); 00221 00223 const Num_T &operator()(int R,int C) const 00224 { it_assert0(R>=0 && R<no_rows && C>=0 && C<no_cols, 00225 "Mat<Num_T>::operator(): index out of range"); return data[R+C*no_rows]; } 00227 Num_T &operator()(int R,int C) 00228 { it_assert0(R>=0 && R<no_rows && C>=0 && C<no_cols, 00229 "Mat<Num_T>::operator(): index out of range"); return data[R+C*no_rows]; } 00231 Num_T &operator()(int index) 00232 { it_assert0(index<no_rows*no_cols && index>=0,"Mat<Num_T>::operator(): index out of range"); 00233 return data[index]; } 00235 const Num_T &operator()(int index) const 00236 { it_assert0(index<no_rows*no_cols && index>=0,"Mat<Num_T>::operator(): index out of range"); 00237 return data[index]; } 00239 const Num_T &get(int R,int C) const 00240 { it_assert0(R>=0 && R<no_rows && C>=0 && C<no_cols, 00241 "Mat<Num_T>::get(): index out of range"); return data[R+C*no_rows]; } 00243 void set(int R,int C, const Num_T &v) 00244 { it_assert0(R>=0 && R<no_rows && C>=0 && C<no_cols, 00245 "Mat<Num_T>::set(): index out of range"); data[R+C*no_rows] = v; } 00246 00252 const Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const; 00258 const Mat<Num_T> get(int r1, int r2, int c1, int c2) const; 00259 00261 const Vec<Num_T> get_row(int Index) const ; 00263 const Mat<Num_T> get_rows(int r1, int r2) const; 00265 const Mat<Num_T> get_rows(const Vec<int> &indexlist) const; 00267 const Vec<Num_T> get_col(int Index) const ; 00269 const Mat<Num_T> get_cols(int c1, int c2) const; 00271 const Mat<Num_T> get_cols(const Vec<int> &indexlist) const; 00273 void set_row(int Index, const Vec<Num_T> &invector); 00275 void set_col(int Index, const Vec<Num_T> &invector); 00277 void copy_row(int to, int from); 00279 void copy_col(int to, int from); 00281 void swap_rows(int r1, int r2); 00283 void swap_cols(int c1, int c2); 00284 00286 void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m); 00288 void set_submatrix(int r, int c, const Mat<Num_T> &m); 00290 void set_submatrix(int r1, int r2, int c1, int c2, const Num_T t); 00291 00293 void del_row(int r); 00295 void del_rows(int r1, int r2); 00297 void del_col(int c); 00299 void del_cols(int c1, int c2); 00301 void ins_row(int r, const Vec<Num_T> &in); 00303 void ins_col(int c, const Vec<Num_T> &in); 00305 void append_row(const Vec<Num_T> &in); 00307 void append_col(const Vec<Num_T> &in); 00308 00310 const Mat<Num_T> transpose() const; 00312 const Mat<Num_T> T() const { return this->transpose(); } 00314 const Mat<Num_T> hermitian_transpose() const; 00316 const Mat<Num_T> H() const { return this->hermitian_transpose(); } 00317 00319 friend const Mat<Num_T> concat_horizontal <>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00321 friend const Mat<Num_T> concat_vertical <>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00322 00324 Mat<Num_T>& operator=(Num_T t); 00326 Mat<Num_T>& operator=(const Mat<Num_T> &m); 00328 Mat<Num_T>& operator=(const Vec<Num_T> &v); 00330 Mat<Num_T>& operator=(const char *values); 00331 00333 Mat<Num_T>& operator+=(const Mat<Num_T> &m); 00335 Mat<Num_T>& operator+=(Num_T t); 00337 friend const Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00339 friend const Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t); 00341 friend const Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m); 00342 00344 Mat<Num_T>& operator-=(const Mat<Num_T> &m); 00346 Mat<Num_T>& operator-=(Num_T t); 00348 friend const Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00350 friend const Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t); 00352 friend const Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m); 00354 friend const Mat<Num_T> operator-<>(const Mat<Num_T> &m); 00355 00357 Mat<Num_T>& operator*=(const Mat<Num_T> &m); 00359 Mat<Num_T>& operator*=(Num_T t); 00361 friend const Mat<Num_T> operator*<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00363 friend const Vec<Num_T> operator*<>(const Mat<Num_T> &m, const Vec<Num_T> &v); 00365 friend const Mat<Num_T> operator*<>(const Vec<Num_T> &v, const Mat<Num_T> &m); 00367 friend const Mat<Num_T> operator*<>(const Mat<Num_T> &m, Num_T t); 00369 friend const Mat<Num_T> operator*<>(Num_T t, const Mat<Num_T> &m); 00371 friend const Mat<Num_T> elem_mult <>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00372 00374 Mat<Num_T>& operator/=(Num_T t); 00376 friend const Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t); 00378 Mat<Num_T>& operator/=(const Mat<Num_T> &m); 00380 friend const Mat<Num_T> elem_div <>(const Mat<Num_T> &m1, const Mat<Num_T> &m2); 00381 00383 bool operator==(const Mat<Num_T> &m) const; 00385 bool operator!=(const Mat<Num_T> &m) const; 00386 00388 Num_T &_elem(int R,int C) { return data[R+C*no_rows]; } 00390 const Num_T &_elem(int R,int C) const { return data[R+C*no_rows]; } 00392 Num_T &_elem(int index) { return data[index]; } 00394 const Num_T &_elem(int index) const { return data[index]; } 00395 00397 Num_T *_data() { return data; } 00399 const Num_T *_data() const { return data; } 00401 int _datasize() const { return datasize; } 00402 00403 protected: 00405 void alloc(int rows, int cols) 00406 { 00407 if ( datasize == rows * cols ) { // Reuse the memory 00408 no_rows = rows; no_cols = cols; 00409 return; 00410 } 00411 free(); // Free memory (if any allocated) 00412 if (rows == 0 || cols == 0) 00413 return; 00414 00415 datasize = rows * cols; 00416 create_elements(data, datasize, factory); 00417 no_rows = rows; no_cols = cols; 00418 00419 it_assert1(data, "Mat<Num_T>::alloc: Out of memory"); 00420 } 00421 00423 void free() { delete [] data; init(); } 00424 00426 int datasize, no_rows, no_cols; 00427 00429 Num_T *data; 00430 00432 const Factory &factory; 00433 00434 private: 00435 void init() 00436 { 00437 data = 0; 00438 datasize = no_rows = no_cols = 0; 00439 } 00440 }; 00441 00442 // ------------------------------------------------------------------------------------- 00443 // Type definitions of mat, cmat, imat, smat, and bmat 00444 // ------------------------------------------------------------------------------------- 00445 00450 typedef Mat<double> mat; 00451 00456 typedef Mat<std::complex<double> > cmat; 00457 00462 typedef Mat<int> imat; 00463 00468 typedef Mat<short int> smat; 00469 00476 typedef Mat<bin> bmat; 00477 00478 } //namespace itpp 00479 00480 #include <itpp/base/vec.h> 00481 00482 namespace itpp { 00483 00484 // ------------------------------------------------------------------------------------- 00485 // Declaration of input and output streams for Mat 00486 // ------------------------------------------------------------------------------------- 00487 00492 template <class Num_T> 00493 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m); 00494 00506 template <class Num_T> 00507 std::istream &operator>>(std::istream &is, Mat<Num_T> &m); 00508 00509 // ------------------------------------------------------------------------------------- 00510 // Implementation of templated Mat members and friends 00511 // ------------------------------------------------------------------------------------- 00512 00513 template<class Num_T> inline 00514 Mat<Num_T>::Mat(const Mat<Num_T> &m) : factory(m.factory) 00515 { 00516 init(); 00517 alloc(m.no_rows, m.no_cols); 00518 copy_vector(m.datasize, m.data, data); 00519 } 00520 00521 template<class Num_T> inline 00522 Mat<Num_T>::Mat(const Mat<Num_T> &m, const Factory &f) : factory(f) 00523 { 00524 init(); 00525 alloc(m.no_rows, m.no_cols); 00526 copy_vector(m.datasize, m.data, data); 00527 } 00528 00529 template<class Num_T> inline 00530 Mat<Num_T>::Mat(const Vec<Num_T> &invector, const Factory &f) : factory(f) 00531 { 00532 init(); 00533 set_size(invector.length(), 1, false); 00534 set_col(0,invector); 00535 } 00536 00537 template<class Num_T> inline 00538 Mat<Num_T>::Mat(Num_T *c_array, int rows, int cols, bool RowMajor, const Factory &f) : factory(f) 00539 { 00540 init(); 00541 alloc(rows, cols); 00542 00543 if (!RowMajor) 00544 copy_vector(datasize, c_array, data); 00545 else { // Row Major 00546 Num_T *ptr = c_array; 00547 for (int i=0; i<rows; i++) { 00548 for (int j=0; j<cols; j++) 00549 data[i+j*no_rows] = *(ptr++); 00550 } 00551 } 00552 } 00553 00554 template<class Num_T> inline 00555 void Mat<Num_T>::set_size(int inrow, int incol, bool copy) 00556 { 00557 it_assert1(inrow>=0 && incol>=0, "Mat<Num_T>::set_size: The rows and columns must be >= 0"); 00558 if (no_rows!=inrow || no_cols!=incol) { 00559 if (copy) { 00560 Mat<Num_T> temp(*this); 00561 int i, j; 00562 00563 alloc(inrow, incol); 00564 for (i=0;i<inrow;i++) 00565 for (j=0;j<incol;j++) 00566 data[i+j*inrow] = (i<temp.no_rows && j<temp.no_cols) ? temp(i,j) : Num_T(0); 00567 } else 00568 alloc(inrow, incol); 00569 } 00570 } 00571 00572 template<class Num_T> inline 00573 void Mat<Num_T>::zeros() 00574 { 00575 for(int i=0; i<datasize; i++) 00576 data[i] = Num_T(0); 00577 } 00578 00579 template<class Num_T> inline 00580 void Mat<Num_T>::ones() 00581 { 00582 for(int i=0; i<datasize; i++) 00583 data[i] = Num_T(1); 00584 } 00585 00586 template<> bool Mat<std::complex<double> >::set(const char *values); 00587 template<> bool Mat<bin>::set(const char *values); 00588 00589 template<class Num_T> 00590 bool Mat<Num_T>::set(const char *values) 00591 { 00592 std::istringstream buffer(values); 00593 int rows=0, maxrows=10, cols=0, nocols=0, maxcols=10; 00594 00595 alloc(maxrows, maxcols); 00596 00597 while (buffer.peek()!=EOF) { 00598 rows++; 00599 if (rows > maxrows) { 00600 maxrows=maxrows*2; 00601 set_size(maxrows, maxcols, true); 00602 } 00603 00604 cols=0; 00605 while ( (buffer.peek() != ';') && (buffer.peek() != EOF) ) { 00606 if (buffer.peek()==',') { 00607 buffer.get(); 00608 } else { 00609 cols++; 00610 if (cols > nocols) { 00611 nocols=cols; 00612 if (cols > maxcols) { 00613 maxcols=maxcols*2; 00614 set_size(maxrows, maxcols, true); 00615 } 00616 } 00617 buffer >> this->operator()(rows-1,cols-1); 00618 while (buffer.peek()==' ') { buffer.get(); } 00619 } 00620 } 00621 00622 if (!buffer.eof()) 00623 buffer.get(); 00624 } 00625 set_size(rows, nocols, true); 00626 00627 return true; 00628 } 00629 00630 template<class Num_T> 00631 bool Mat<Num_T>::set(const std::string &str) 00632 { 00633 return set(str.c_str()); 00634 } 00635 00636 template<class Num_T> inline 00637 const Mat<Num_T> Mat<Num_T>::operator()(int r1, int r2, int c1, int c2) const 00638 { 00639 if (r1 == -1) r1 = no_rows-1; 00640 if (r2 == -1) r2 = no_rows-1; 00641 if (c1 == -1) c1 = no_cols-1; 00642 if (c2 == -1) c2 = no_cols-1; 00643 00644 it_assert1(r1>=0 && r2>=0 && r1<no_rows && r2<no_rows && 00645 c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, "operator()(r1,r2,c1,c2)"); 00646 00647 it_assert1(r2>=r1 && c2>=c1, "Mat<Num_T>::op(): r2>=r1 or c2>=c1 not fulfilled"); 00648 00649 Mat<Num_T> s(r2-r1+1, c2-c1+1); 00650 00651 for (int i=0;i<s.no_cols;i++) 00652 copy_vector(s.no_rows, data+r1+(c1+i)*no_rows, s.data+i*s.no_rows); 00653 00654 return s; 00655 } 00656 00657 template<class Num_T> inline 00658 const Mat<Num_T> Mat<Num_T>::get(int r1, int r2, int c1, int c2) const 00659 { 00660 return (*this)(r1, r2, c1, c2); 00661 } 00662 00663 template<class Num_T> inline 00664 const Vec<Num_T> Mat<Num_T>::get_row(int Index) const 00665 { 00666 it_assert1(Index>=0 && Index<no_rows, "Mat<Num_T>::get_row: index out of range"); 00667 Vec<Num_T> a(no_cols); 00668 00669 copy_vector(no_cols, data+Index, no_rows, a._data(), 1); 00670 return a; 00671 } 00672 00673 template<class Num_T> 00674 const Mat<Num_T> Mat<Num_T>::get_rows(int r1, int r2) const 00675 { 00676 it_assert1(r1>=0 && r2<no_rows && r1<=r2, "Mat<Num_T>::get_rows: index out of range"); 00677 Mat<Num_T> m(r2-r1+1, no_cols); 00678 00679 for (int i=0; i<m.rows(); i++) 00680 copy_vector(no_cols, data+i+r1, no_rows, m.data+i, m.no_rows); 00681 00682 return m; 00683 } 00684 00685 template<class Num_T> 00686 const Mat<Num_T> Mat<Num_T>::get_rows(const Vec<int> &indexlist) const 00687 { 00688 Mat<Num_T> m(indexlist.size(),no_cols); 00689 00690 for (int i=0;i<indexlist.size();i++) { 00691 it_assert1(indexlist(i)>=0 && indexlist(i)<no_rows, "Mat<Num_T>::get_rows: index out of range"); 00692 copy_vector(no_cols, data+indexlist(i), no_rows, m.data+i, m.no_rows); 00693 } 00694 00695 return m; 00696 } 00697 00698 template<class Num_T> inline 00699 const Vec<Num_T> Mat<Num_T>::get_col(int Index) const 00700 { 00701 it_assert1(Index>=0 && Index<no_cols, "Mat<Num_T>::get_col: index out of range"); 00702 Vec<Num_T> a(no_rows); 00703 00704 copy_vector(no_rows, data+Index*no_rows, a._data()); 00705 00706 return a; 00707 } 00708 00709 template<class Num_T> inline 00710 const Mat<Num_T> Mat<Num_T>::get_cols(int c1, int c2) const 00711 { 00712 it_assert1(c1>=0 && c2<no_cols && c1<=c2, "Mat<Num_T>::get_cols: index out of range"); 00713 Mat<Num_T> m(no_rows, c2-c1+1); 00714 00715 for (int i=0; i<m.cols(); i++) 00716 copy_vector(no_rows, data+(i+c1)*no_rows, m.data+i*m.no_rows); 00717 00718 return m; 00719 } 00720 00721 template<class Num_T> inline 00722 const Mat<Num_T> Mat<Num_T>::get_cols(const Vec<int> &indexlist) const 00723 { 00724 Mat<Num_T> m(no_rows,indexlist.size()); 00725 00726 for (int i=0; i<indexlist.size(); i++) { 00727 it_assert1(indexlist(i)>=0 && indexlist(i)<no_cols, "Mat<Num_T>::get_cols: index out of range"); 00728 copy_vector(no_rows, data+indexlist(i)*no_rows, m.data+i*m.no_rows); 00729 } 00730 00731 return m; 00732 } 00733 00734 template<class Num_T> inline 00735 void Mat<Num_T>::set_row(int Index, const Vec<Num_T> &v) 00736 { 00737 it_assert1(Index>=0 && Index<no_rows, "Mat<Num_T>::set_row: index out of range"); 00738 it_assert1(v.length() == no_cols, "Mat<Num_T>::set_row: lengths doesn't match"); 00739 00740 copy_vector(v.size(), v._data(), 1, data+Index, no_rows); 00741 } 00742 00743 template<class Num_T> inline 00744 void Mat<Num_T>::set_col(int Index, const Vec<Num_T> &v) 00745 { 00746 it_assert1(Index>=0 && Index<no_cols, "Mat<Num_T>::set_col: index out of range"); 00747 it_assert1(v.length() == no_rows, "Mat<Num_T>::set_col: lengths doesn't match"); 00748 00749 copy_vector(v.size(), v._data(), data+Index*no_rows); 00750 } 00751 00752 template<class Num_T> inline 00753 void Mat<Num_T>::copy_row(int to, int from) 00754 { 00755 it_assert1(to>=0 && from>=0 && to<no_rows && from<no_rows, 00756 "Mat<Num_T>::copy_row: index out of range"); 00757 00758 if (from == to) 00759 return; 00760 00761 copy_vector(no_cols, data+from, no_rows, data+to, no_rows); 00762 } 00763 00764 template<class Num_T> inline 00765 void Mat<Num_T>::copy_col(int to, int from) 00766 { 00767 it_assert1(to>=0 && from>=0 && to<no_cols && from<no_cols, 00768 "Mat<Num_T>::copy_col: index out of range"); 00769 00770 if (from == to) 00771 return; 00772 00773 copy_vector(no_rows, data+from*no_rows, data+to*no_rows); 00774 } 00775 00776 template<class Num_T> inline 00777 void Mat<Num_T>::swap_rows(int r1, int r2) 00778 { 00779 it_assert1(r1>=0 && r2>=0 && r1<no_rows && r2<no_rows, 00780 "Mat<Num_T>::swap_rows: index out of range"); 00781 00782 if (r1 == r2) 00783 return; 00784 00785 swap_vector(no_cols, data+r1, no_rows, data+r2, no_rows); 00786 } 00787 00788 template<class Num_T> inline 00789 void Mat<Num_T>::swap_cols(int c1, int c2) 00790 { 00791 it_assert1(c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, 00792 "Mat<Num_T>::swap_cols: index out of range"); 00793 00794 if (c1 == c2) 00795 return; 00796 00797 swap_vector(no_rows, data+c1*no_rows, data+c2*no_rows); 00798 } 00799 00800 template<class Num_T> inline 00801 void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m) 00802 { 00803 00804 if (r1 == -1) r1 = no_rows-1; 00805 if (r2 == -1) r2 = no_rows-1; 00806 if (c1 == -1) c1 = no_cols-1; 00807 if (c2 == -1) c2 = no_cols-1; 00808 00809 it_assert1(r1>=0 && r2>=0 && r1<no_rows && r2<no_rows && 00810 c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, "Mat<Num_T>::set_submatrix(): index out of range"); 00811 00812 it_assert1(r2>=r1 && c2>=c1, "Mat<Num_T>::set_submatrix: r2<r1 or c2<c1"); 00813 it_assert1(m.no_rows == r2-r1+1 && m.no_cols == c2-c1+1, "Mat<Num_T>::set_submatrix(): sizes don't match"); 00814 00815 for (int i=0; i<m.no_cols; i++) 00816 copy_vector(m.no_rows, m.data+i*m.no_rows, data+(c1+i)*no_rows+r1); 00817 } 00818 00819 00820 00821 template<class Num_T> inline 00822 void Mat<Num_T>::set_submatrix(int r, int c, const Mat<Num_T> &m) 00823 { 00824 00825 it_assert1(r>=0 && r+m.no_rows<=no_rows && 00826 c>=0 && c+m.no_cols<=no_cols, "Mat<Num_T>::set_submatrix(): index out of range"); 00827 00828 for (int i=0; i<m.no_cols; i++) 00829 copy_vector(m.no_rows, m.data+i*m.no_rows, data+(c+i)*no_rows+r); 00830 } 00831 00832 00833 00834 template<class Num_T> inline 00835 void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, const Num_T t) 00836 { 00837 00838 if (r1 == -1) r1 = no_rows-1; 00839 if (r2 == -1) r2 = no_rows-1; 00840 if (c1 == -1) c1 = no_cols-1; 00841 if (c2 == -1) c2 = no_cols-1; 00842 00843 it_assert1(r1>=0 && r2>=0 && r1<no_rows && r2<no_rows && 00844 c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, "Mat<Num_T>::set_submatrix(): i\ 00845 ndex out of range"); 00846 00847 it_assert1(r2>=r1 && c2>=c1, "Mat<Num_T>::set_submatrix: r2<r1 or c2<c1"); 00848 00849 int i, j, pos, rows = r2-r1+1; 00850 00851 for (i=c1; i<=c2; i++) { 00852 pos = i*no_rows+r1; 00853 for (j=0; j<rows; j++) { 00854 data[pos++] = t; 00855 } 00856 } 00857 } 00858 00859 template<class Num_T> inline 00860 void Mat<Num_T>::del_row(int r) 00861 { 00862 it_assert1(r>=0 && r<no_rows, "Mat<Num_T>::del_row(): index out of range"); 00863 Mat<Num_T> Temp(*this); 00864 set_size(no_rows-1, no_cols, false); 00865 for (int i=0 ; i < r ; i++) { 00866 copy_vector(no_cols, &Temp.data[i], no_rows+1, &data[i], no_rows); 00867 } 00868 for (int i=r ; i < no_rows ; i++) { 00869 copy_vector(no_cols, &Temp.data[i+1], no_rows+1, &data[i], no_rows); 00870 } 00871 00872 } 00873 00874 template<class Num_T> inline 00875 void Mat<Num_T>::del_rows(int r1, int r2) 00876 { 00877 it_assert1(r1>=0 && r2<no_rows && r1<=r2, "Mat<Num_T>::del_rows(): index out of range"); 00878 00879 for (int i=r1 ; i<=r2 ; i++) { 00880 del_row(r1); 00881 } 00882 00883 // this should work, I don't understand why it does not. 00884 /* Mat<Num_T> Temp(*this); */ 00885 /* int n_deleted_rows = r2-r1+1; */ 00886 /* set_size(no_rows-n_deleted_rows, no_cols, false); */ 00887 /* for (int i=0 ; i < r1 ; i++) { */ 00888 /* copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows); */ 00889 /* } */ 00890 /* for (int i=r2 ; i < no_rows ; i++) { */ 00891 /* copy_vector(no_cols, &Temp.data[i+1], Temp.no_rows, &data[i-n_deleted_rows+1], no_rows); */ 00892 /* } */ 00893 } 00894 00895 template<class Num_T> inline 00896 void Mat<Num_T>::del_col(int c) 00897 { 00898 it_assert1(c>=0 && c<no_cols, "Mat<Num_T>::del_col(): index out of range"); 00899 Mat<Num_T> Temp(*this); 00900 00901 set_size(no_rows, no_cols-1, false); 00902 copy_vector(c*no_rows, Temp.data, data); 00903 copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]); 00904 } 00905 00906 template<class Num_T> inline 00907 void Mat<Num_T>::del_cols(int c1, int c2) 00908 { 00909 it_assert1(c1>=0 && c2<no_cols && c1<=c2, "Mat<Num_T>::del_cols(): index out of range"); 00910 Mat<Num_T> Temp(*this); 00911 int n_deleted_cols = c2-c1+1; 00912 set_size(no_rows, no_cols-n_deleted_cols, false); 00913 copy_vector(c1*no_rows, Temp.data, data); 00914 copy_vector((no_cols-c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]); 00915 } 00916 00917 template<class Num_T> inline 00918 void Mat<Num_T>::ins_row(int r, const Vec<Num_T> &in) 00919 { 00920 it_assert1(r>=0 && r<=no_rows, "Mat<Num_T>::ins_row(): index out of range"); 00921 it_assert1((in.size() == no_cols) || no_rows==0, "Mat<Num_T>::ins_row(): vector size does not match"); 00922 00923 if (no_cols==0) { 00924 no_cols = in.size(); 00925 } 00926 00927 Mat<Num_T> Temp(*this); 00928 set_size(no_rows+1, no_cols, false); 00929 00930 for (int i=0 ; i < r ; i++) { 00931 copy_vector(no_cols, &Temp.data[i], no_rows-1, &data[i], no_rows); 00932 } 00933 copy_vector(no_cols, in._data(), 1, &data[r], no_rows); 00934 for (int i=r+1 ; i < no_rows ; i++) { 00935 copy_vector(no_cols, &Temp.data[i-1], no_rows-1, &data[i], no_rows); 00936 } 00937 } 00938 00939 template<class Num_T> inline 00940 void Mat<Num_T>::ins_col(int c, const Vec<Num_T> &in) 00941 { 00942 it_assert1(c>=0 && c<=no_cols, "Mat<Num_T>::ins_col(): index out of range"); 00943 it_assert1(in.size() == no_rows || no_cols==0, "Mat<Num_T>::ins_col(): vector size does not match"); 00944 00945 if (no_rows==0) { 00946 no_rows = in.size(); 00947 } 00948 00949 Mat<Num_T> Temp(*this); 00950 set_size(no_rows, no_cols+1, false); 00951 00952 copy_vector(c*no_rows, Temp.data, data); 00953 copy_vector(no_rows, in._data(), &data[c*no_rows]); 00954 copy_vector((no_cols-c-1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]); 00955 } 00956 00957 template<class Num_T> inline 00958 void Mat<Num_T>::append_row(const Vec<Num_T> &in) 00959 { 00960 ins_row(no_rows, in); 00961 } 00962 00963 template<class Num_T> inline 00964 void Mat<Num_T>::append_col(const Vec<Num_T> &in) 00965 { 00966 ins_col(no_cols, in); 00967 } 00968 00969 template<class Num_T> 00970 const Mat<Num_T> Mat<Num_T>::transpose() const 00971 { 00972 Mat<Num_T> temp(no_cols, no_rows); 00973 for (int i=0; i<no_rows; i++) 00974 for (int j=0; j<no_cols; j++) 00975 temp(j,i) = operator()(i,j); 00976 00977 return temp; 00978 } 00979 00980 template<> const cmat Mat<std::complex<double> >::hermitian_transpose() const; 00981 00982 template<class Num_T> 00983 const Mat<Num_T> Mat<Num_T>::hermitian_transpose() const 00984 { 00985 Mat<Num_T> temp(no_cols, no_rows); 00986 for (int i=0; i<no_rows; i++) 00987 for (int j=0; j<no_cols; j++) 00988 temp(j,i) = operator()(i,j); 00989 00990 return temp; 00991 } 00992 00993 template<class Num_T> 00994 const Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 00995 { 00996 it_assert1(m1.no_rows == m2.no_rows, "Mat<Num_T>::concat_horizontal; wrong sizes"); 00997 00998 Mat<Num_T> temp(m1.no_rows, m1.no_cols+m2.no_cols); 00999 int i; 01000 01001 for (i=0; i<m1.no_cols; i++) { 01002 temp.set_col(i, m1.get_col(i)); 01003 } 01004 for (i=0; i<m2.no_cols; i++) { 01005 temp.set_col(i+m1.no_cols, m2.get_col(i)); 01006 } 01007 01008 return temp; 01009 } 01010 01011 template<class Num_T> 01012 const Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01013 { 01014 it_assert1(m1.no_cols == m2.no_cols, "Mat<Num_T>::concat_vertical; wrong sizes"); 01015 01016 Mat<Num_T> temp(m1.no_rows+m2.no_rows, m1.no_cols); 01017 int i; 01018 01019 for (i=0; i<m1.no_rows; i++) { 01020 temp.set_row(i, m1.get_row(i)); 01021 } 01022 for (i=0; i<m2.no_rows; i++) { 01023 temp.set_row(i+m1.no_rows, m2.get_row(i)); 01024 } 01025 01026 return temp; 01027 } 01028 01029 template<class Num_T> inline 01030 Mat<Num_T>& Mat<Num_T>::operator=(Num_T t) 01031 { 01032 for (int i=0; i<datasize; i++) 01033 data[i] = t; 01034 return *this; 01035 } 01036 01037 template<class Num_T> inline 01038 Mat<Num_T>& Mat<Num_T>::operator=(const Mat<Num_T> &m) 01039 { 01040 if (this != &m) { 01041 set_size(m.no_rows,m.no_cols, false); 01042 if (m.datasize != 0) 01043 copy_vector(m.datasize, m.data, data); 01044 } 01045 return *this; 01046 } 01047 01048 template<class Num_T> inline 01049 Mat<Num_T>& Mat<Num_T>::operator=(const Vec<Num_T> &v) 01050 { 01051 it_assert1( (no_rows == 1 && no_cols == v.size()) || (no_cols == 1 && no_rows == v.size()), 01052 "Mat<Num_T>::operator=(vec); wrong size"); 01053 01054 // construct a 1-d column of the vector 01055 set_size(v.size(), 1, false); 01056 set_col(0, v); 01057 01058 return *this; 01059 } 01060 01061 template<class Num_T> inline 01062 Mat<Num_T>& Mat<Num_T>::operator=(const char *values) 01063 { 01064 set(values); 01065 return *this; 01066 } 01067 01068 //-------------------- Templated friend functions -------------------------- 01069 01070 template<class Num_T> inline 01071 Mat<Num_T>& Mat<Num_T>::operator+=(const Mat<Num_T> &m) 01072 { 01073 if (datasize == 0) 01074 operator=(m); 01075 else { 01076 int i, j, m_pos=0, pos=0; 01077 it_assert1(m.no_rows==no_rows && m.no_cols==no_cols,"Mat<Num_T>::operator+=: wrong sizes"); 01078 for (i=0; i<no_cols; i++) { 01079 for (j=0; j<no_rows; j++) 01080 data[pos+j] += m.data[m_pos+j]; 01081 pos += no_rows; 01082 m_pos += m.no_rows; 01083 } 01084 } 01085 return *this; 01086 } 01087 01088 template<class Num_T> inline 01089 Mat<Num_T>& Mat<Num_T>::operator+=(Num_T t) 01090 { 01091 for (int i=0; i<datasize; i++) 01092 data[i] += t; 01093 return *this; 01094 } 01095 01096 template<class Num_T> inline 01097 const Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01098 { 01099 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01100 int i, j, m1_pos=0, m2_pos=0, r_pos=0; 01101 01102 it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat<Num_T>::operator+: wrong sizes"); 01103 01104 for (i=0; i<r.no_cols; i++) { 01105 for (j=0; j<r.no_rows; j++) 01106 r.data[r_pos+j] = m1.data[m1_pos+j] + m2.data[m2_pos+j]; 01107 // next column 01108 m1_pos += m1.no_rows; 01109 m2_pos += m2.no_rows; 01110 r_pos += r.no_rows; 01111 } 01112 01113 return r; 01114 } 01115 01116 01117 template<class Num_T> inline 01118 const Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t) 01119 { 01120 Mat<Num_T> r(m.no_rows, m.no_cols); 01121 01122 for (int i=0; i<r.datasize; i++) 01123 r.data[i] = m.data[i] + t; 01124 01125 return r; 01126 } 01127 01128 template<class Num_T> inline 01129 const Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m) 01130 { 01131 Mat<Num_T> r(m.no_rows, m.no_cols); 01132 01133 for (int i=0; i<r.datasize; i++) 01134 r.data[i] = t + m.data[i]; 01135 01136 return r; 01137 } 01138 01139 template<class Num_T> inline 01140 Mat<Num_T>& Mat<Num_T>::operator-=(const Mat<Num_T> &m) 01141 { 01142 int i,j, m_pos=0, pos=0; 01143 01144 if (datasize == 0) { 01145 set_size(m.no_rows, m.no_cols, false); 01146 for (i=0; i<no_cols; i++) { 01147 for (j=0; j<no_rows; j++) 01148 data[pos+j] = -m.data[m_pos+j]; 01149 // next column 01150 m_pos += m.no_rows; 01151 pos += no_rows; 01152 } 01153 } else { 01154 it_assert1(m.no_rows==no_rows && m.no_cols==no_cols,"Mat<Num_T>::operator-=: wrong sizes"); 01155 for (i=0; i<no_cols; i++) { 01156 for (j=0; j<no_rows; j++) 01157 data[pos+j] -= m.data[m_pos+j]; 01158 // next column 01159 m_pos += m.no_rows; 01160 pos += no_rows; 01161 } 01162 } 01163 return *this; 01164 } 01165 01166 template<class Num_T> inline 01167 const Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01168 { 01169 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01170 int i, j, m1_pos=0, m2_pos=0, r_pos=0; 01171 01172 it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat<Num_T>::operator-: wrong sizes"); 01173 01174 for (i=0; i<r.no_cols; i++) { 01175 for (j=0; j<r.no_rows; j++) 01176 r.data[r_pos+j] = m1.data[m1_pos+j] - m2.data[m2_pos+j]; 01177 // next column 01178 m1_pos += m1.no_rows; 01179 m2_pos += m2.no_rows; 01180 r_pos += r.no_rows; 01181 } 01182 01183 return r; 01184 } 01185 01186 template<class Num_T> inline 01187 Mat<Num_T>& Mat<Num_T>::operator-=(Num_T t) 01188 { 01189 for (int i=0; i<datasize; i++) 01190 data[i] -= t; 01191 return *this; 01192 } 01193 01194 template<class Num_T> inline 01195 const Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t) 01196 { 01197 Mat<Num_T> r(m.no_rows, m.no_cols); 01198 int i, j, m_pos=0, r_pos=0; 01199 01200 for (i=0; i<r.no_cols; i++) { 01201 for (j=0; j<r.no_rows; j++) 01202 r.data[r_pos+j] = m.data[m_pos+j] - t; 01203 // next column 01204 m_pos += m.no_rows; 01205 r_pos += r.no_rows; 01206 } 01207 01208 return r; 01209 } 01210 01211 template<class Num_T> inline 01212 const Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m) 01213 { 01214 Mat<Num_T> r(m.no_rows, m.no_cols); 01215 int i, j, m_pos=0, r_pos=0; 01216 01217 for (i=0; i<r.no_cols; i++) { 01218 for (j=0; j<r.no_rows; j++) 01219 r.data[r_pos+j] = t - m.data[m_pos+j]; 01220 // next column 01221 m_pos += m.no_rows; 01222 r_pos += r.no_rows; 01223 } 01224 01225 return r; 01226 } 01227 01228 template<class Num_T> inline 01229 const Mat<Num_T> operator-(const Mat<Num_T> &m) 01230 { 01231 Mat<Num_T> r(m.no_rows, m.no_cols); 01232 int i, j, m_pos=0, r_pos=0; 01233 01234 for (i=0; i<r.no_cols; i++) { 01235 for (j=0; j<r.no_rows; j++) 01236 r.data[r_pos+j] = -m.data[m_pos+j]; 01237 // next column 01238 m_pos += m.no_rows; 01239 r_pos += r.no_rows; 01240 } 01241 01242 return r; 01243 } 01244 01245 #if defined(HAVE_CBLAS) 01246 template<> mat& Mat<double>::operator*=(const mat &m); 01247 template<> cmat& Mat<std::complex<double> >::operator*=(const cmat &m); 01248 #endif 01249 01250 template<class Num_T> inline 01251 Mat<Num_T>& Mat<Num_T>::operator*=(const Mat<Num_T> &m) 01252 { 01253 it_assert1(no_cols == m.no_rows,"Mat<Num_T>::operator*=: wrong sizes"); 01254 Mat<Num_T> r(no_rows, m.no_cols); 01255 01256 Num_T tmp; 01257 01258 int i,j,k, r_pos=0, pos=0, m_pos=0; 01259 01260 for (i=0; i<r.no_cols; i++) { 01261 for (j=0; j<r.no_rows; j++) { 01262 tmp = Num_T(0); 01263 pos = 0; 01264 for (k=0; k<no_cols; k++) { 01265 tmp += data[pos+j] * m.data[m_pos+k]; 01266 pos += no_rows; 01267 } 01268 r.data[r_pos+j] = tmp; 01269 } 01270 r_pos += r.no_rows; 01271 m_pos += m.no_rows; 01272 } 01273 operator=(r); // time consuming 01274 return *this; 01275 } 01276 01277 template<class Num_T> inline 01278 Mat<Num_T>& Mat<Num_T>::operator*=(Num_T t) 01279 { 01280 for (int i=0; i<datasize; i++) 01281 data[i] *= t; 01282 return *this; 01283 } 01284 01285 #if defined(HAVE_CBLAS) 01286 template<> const mat operator*(const mat &m1, const mat &m2); 01287 template<> const cmat operator*(const cmat &m1, const cmat &m2); 01288 #endif 01289 01290 01291 template<class Num_T> 01292 const Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01293 { 01294 it_assert1(m1.no_cols == m2.no_rows,"Mat<Num_T>::operator*: wrong sizes"); 01295 Mat<Num_T> r(m1.no_rows, m2.no_cols); 01296 01297 Num_T tmp; 01298 int i, j, k; 01299 Num_T *tr=r.data, *t1, *t2=m2.data; 01300 01301 for (i=0; i<r.no_cols; i++) { 01302 for (j=0; j<r.no_rows; j++) { 01303 tmp = Num_T(0); t1 = m1.data+j; 01304 for (k=m1.no_cols; k>0; k--) { 01305 tmp += *(t1) * *(t2++); 01306 t1 += m1.no_rows; 01307 } 01308 *(tr++) = tmp; t2 -= m2.no_rows; 01309 } 01310 t2 += m2.no_rows; 01311 } 01312 01313 return r; 01314 } 01315 01316 #if defined(HAVE_CBLAS) 01317 template<> const vec operator*(const mat &m, const vec &v); 01318 template<> const cvec operator*(const cmat &m, const cvec &v); 01319 #endif 01320 01321 template<class Num_T> 01322 const Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v) 01323 { 01324 it_assert1(m.no_cols == v.size(),"Mat<Num_T>::operator*: wrong sizes"); 01325 Vec<Num_T> r(m.no_rows); 01326 int i, k, m_pos; 01327 01328 for (i=0; i<m.no_rows; i++) { 01329 r(i) = Num_T(0); 01330 m_pos = 0; 01331 for (k=0; k<m.no_cols; k++) { 01332 r(i) += m.data[m_pos+i] * v(k); 01333 m_pos += m.no_rows; 01334 } 01335 } 01336 01337 return r; 01338 } 01339 01340 template<class Num_T> inline 01341 const Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m) 01342 { 01343 it_assert1(m.no_rows == 1,"Mat<Num_T>::operator*: wrong sizes"); 01344 Mat<Num_T> r(v.size(), m.no_cols); 01345 01346 for (int i = 0; i < v.size(); ++i) 01347 for (int j = 0; j < m.no_cols; ++j) 01348 r(i,j) = v(i) * m.data[j]; 01349 01350 return r; 01351 } 01352 01353 template<class Num_T> inline 01354 const Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t) 01355 { 01356 Mat<Num_T> r(m.no_rows, m.no_cols); 01357 01358 for (int i=0; i<r.datasize; i++) 01359 r.data[i] = m.data[i] * t; 01360 01361 return r; 01362 } 01363 01364 template<class Num_T> inline 01365 const Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m) 01366 { 01367 Mat<Num_T> r(m.no_rows, m.no_cols); 01368 01369 for (int i=0; i<r.datasize; i++) 01370 r.data[i] = m.data[i] * t; 01371 01372 return r; 01373 } 01374 01375 template<class Num_T> inline 01376 const Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01377 { 01378 it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat<Num_T>::elem_mult: wrong sizes"); 01379 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01380 01381 for (int i=0; i<r.datasize; i++) 01382 r.data[i] = m1.data[i] * m2.data[i]; 01383 01384 return r; 01385 } 01386 01387 template<class Num_T> inline 01388 Mat<Num_T>& Mat<Num_T>::operator/=(Num_T t) 01389 { 01390 for (int i=0; i<datasize; i++) 01391 data[i] /= t; 01392 return *this; 01393 } 01394 01395 template<class Num_T> inline 01396 const Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t) 01397 { 01398 Mat<Num_T> r(m.no_rows, m.no_cols); 01399 01400 for (int i=0; i<r.datasize; i++) 01401 r.data[i] = m.data[i] / t; 01402 01403 return r; 01404 } 01405 01406 template<class Num_T> inline 01407 Mat<Num_T>& Mat<Num_T>::operator/=(const Mat<Num_T> &m) 01408 { 01409 it_assert1(m.no_rows==no_rows && m.no_cols==no_cols, "Mat<Num_T>::operator/=: wrong sizes"); 01410 01411 for (int i=0; i<datasize; i++) 01412 data[i] /= m.data[i]; 01413 return *this; 01414 } 01415 01416 template<class Num_T> inline 01417 const Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2) 01418 { 01419 it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat<Num_T>::elem_div: worng sizes"); 01420 Mat<Num_T> r(m1.no_rows, m1.no_cols); 01421 01422 for (int i=0; i<r.datasize; i++) 01423 r.data[i] = m1.data[i] / m2.data[i]; 01424 01425 return r; 01426 } 01427 01428 template<class Num_T> 01429 bool Mat<Num_T>::operator==(const Mat<Num_T> &m) const 01430 { 01431 if (no_rows!=m.no_rows || no_cols != m.no_cols) return false; 01432 for (int i=0;i<datasize;i++) { 01433 if (data[i]!=m.data[i]) return false; 01434 } 01435 return true; 01436 } 01437 01438 template<class Num_T> 01439 bool Mat<Num_T>::operator!=(const Mat<Num_T> &m) const 01440 { 01441 if (no_rows != m.no_rows || no_cols != m.no_cols) return true; 01442 for (int i=0;i<datasize;i++) { 01443 if (data[i]!=m.data[i]) return true; 01444 } 01445 return false; 01446 } 01447 01448 template <class Num_T> 01449 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m) 01450 { 01451 int i; 01452 01453 switch (m.rows()) { 01454 case 0 : 01455 os << "[]"; 01456 break; 01457 case 1 : 01458 os << '[' << m.get_row(0) << ']'; 01459 break; 01460 default: 01461 os << '[' << m.get_row(0) << std::endl; 01462 for (i=1; i<m.rows()-1; i++) 01463 os << ' ' << m.get_row(i) << std::endl; 01464 os << ' ' << m.get_row(m.rows()-1) << ']'; 01465 } 01466 01467 return os; 01468 } 01469 01470 template <class Num_T> 01471 std::istream &operator>>(std::istream &is, Mat<Num_T> &m) 01472 { 01473 std::ostringstream buffer; 01474 bool started = false; 01475 bool finished = false; 01476 bool brackets = false; 01477 bool within_double_brackets = false; 01478 char c; 01479 01480 while (!finished) { 01481 if (is.eof()) { 01482 finished = true; 01483 } else { 01484 c = is.get(); 01485 01486 if (is.eof() || (c == '\n')) { 01487 if (brackets) { 01488 // Right bracket missing 01489 is.setstate(std::ios_base::failbit); 01490 finished = true; 01491 } else if (!((c == '\n') && !started)) { 01492 finished = true; 01493 } 01494 } else if ((c == ' ') || (c == '\t')) { 01495 if (started) { 01496 buffer << ' '; 01497 } 01498 } else if (c == '[') { 01499 if ((started && !brackets) || within_double_brackets) { 01500 // Unexpected left bracket 01501 is.setstate(std::ios_base::failbit); 01502 finished = true; 01503 } else if (!started) { 01504 started = true; 01505 brackets = true; 01506 } else { 01507 within_double_brackets = true; 01508 } 01509 } else if (c == ']') { 01510 if (!started || !brackets) { 01511 // Unexpected right bracket 01512 is.setstate(std::ios_base::failbit); 01513 finished = true; 01514 } else if (within_double_brackets) { 01515 within_double_brackets = false; 01516 buffer << ';'; 01517 } else { 01518 finished = true; 01519 } 01520 while (!is.eof() && (((c = is.peek()) == ' ') || (c == '\t'))) { 01521 is.get(); 01522 } 01523 if (!is.eof() && (c == '\n')) { 01524 is.get(); 01525 } 01526 } else { 01527 started = true; 01528 buffer << c; 01529 } 01530 } 01531 } 01532 01533 if (!started) { 01534 m.set_size(0, false); 01535 } else { 01536 m.set(buffer.str()); 01537 } 01538 01539 return is; 01540 } 01541 01542 #ifndef _MSC_VER 01543 01544 //--------------------------------------------------------------------- 01545 // Instantiations 01546 //--------------------------------------------------------------------- 01547 01548 //------- class instantiations -------- 01549 01551 extern template class Mat<double>; 01553 extern template class Mat<std::complex<double> >; 01555 extern template class Mat<int>; 01557 extern template class Mat<short int>; 01559 extern template class Mat<bin>; 01560 01561 01562 //-------------------- Operator instantiations -------------------- 01563 01564 //-------- Addition operators --------------- 01565 01567 extern template const mat operator+(const mat &m1, const mat &m2); 01569 extern template const cmat operator+(const cmat &m1, const cmat &m2); 01571 extern template const imat operator+(const imat &m1, const imat &m2); 01573 extern template const smat operator+(const smat &m1, const smat &m2); 01575 extern template const bmat operator+(const bmat &m1, const bmat &m2); 01576 01578 extern template const mat operator+(const mat &m, double t); 01580 extern template const cmat operator+(const cmat &m, std::complex<double> t); 01582 extern template const imat operator+(const imat &m, int t); 01584 extern template const smat operator+(const smat &m, short t); 01586 extern template const bmat operator+(const bmat &m, bin t); 01587 01589 extern template const mat operator+(double t, const mat &m); 01591 extern template const cmat operator+(std::complex<double> t, const cmat &m); 01593 extern template const imat operator+(int t, const imat &m); 01595 extern template const smat operator+(short t, const smat &m); 01597 extern template const bmat operator+(bin t, const bmat &m); 01598 01599 //-------- Subraction operators --------------- 01600 01602 extern template const mat operator-(const mat &m1, const mat &m2); 01604 extern template const cmat operator-(const cmat &m1, const cmat &m2); 01606 extern template const imat operator-(const imat &m1, const imat &m2); 01608 extern template const smat operator-(const smat &m1, const smat &m2); 01610 extern template const bmat operator-(const bmat &m1, const bmat &m2); 01611 01613 extern template const mat operator-(const mat &m, double t); 01615 extern template const cmat operator-(const cmat &m, std::complex<double> t); 01617 extern template const imat operator-(const imat &m, int t); 01619 extern template const smat operator-(const smat &m, short t); 01621 extern template const bmat operator-(const bmat &m, bin t); 01622 01624 extern template const mat operator-(double t, const mat &m); 01626 extern template const cmat operator-(std::complex<double> t, const cmat &m); 01628 extern template const imat operator-(int t, const imat &m); 01630 extern template const smat operator-(short t, const smat &m); 01632 extern template const bmat operator-(bin t, const bmat &m); 01633 01634 //--------- Unary minus --------------- 01635 01637 extern template const mat operator-(const mat &m); 01639 extern template const cmat operator-(const cmat &m); 01641 extern template const imat operator-(const imat &m); 01643 extern template const smat operator-(const smat &m); 01645 extern template const bmat operator-(const bmat &m); 01646 01647 //-------- Multiplication operators --------------- 01648 01649 #if !defined(HAVE_CBLAS) 01651 extern template const mat operator*(const mat &m1, const mat &m2); 01653 extern template const cmat operator*(const cmat &m1, const cmat &m2); 01654 #endif 01656 extern template const imat operator*(const imat &m1, const imat &m2); 01658 extern template const smat operator*(const smat &m1, const smat &m2); 01660 extern template const bmat operator*(const bmat &m1, const bmat &m2); 01661 01662 #if !defined(HAVE_CBLAS) 01664 extern template const vec operator*(const mat &m, const vec &v); 01666 extern template const cvec operator*(const cmat &m, const cvec &v); 01667 #endif 01669 extern template const ivec operator*(const imat &m, const ivec &v); 01671 extern template const svec operator*(const smat &m, const svec &v); 01673 extern template const bvec operator*(const bmat &m, const bvec &v); 01674 01676 extern template const mat operator*(const vec &v, const mat &m); 01678 extern template const cmat operator*(const cvec &v, const cmat &m); 01680 extern template const imat operator*(const ivec &v, const imat &m); 01682 extern template const smat operator*(const svec &v, const smat &m); 01684 extern template const bmat operator*(const bvec &v, const bmat &m); 01685 01687 extern template const mat operator*(const mat &m, double t); 01689 extern template const cmat operator*(const cmat &m, std::complex<double> t); 01691 extern template const imat operator*(const imat &m, int t); 01693 extern template const smat operator*(const smat &m, short t); 01695 extern template const bmat operator*(const bmat &m, bin t); 01696 01698 extern template const mat operator*(double t, const mat &m); 01700 extern template const cmat operator*(std::complex<double> t, const cmat &m); 01702 extern template const imat operator*(int t, const imat &m); 01704 extern template const smat operator*(short t, const smat &m); 01706 extern template const bmat operator*(bin t, const bmat &m); 01707 01708 // ------------ Elementwise multiplication ----------- 01709 01711 extern template const mat elem_mult(const mat &m1, const mat &m2); 01713 extern template const cmat elem_mult(const cmat &m1, const cmat &m2); 01715 extern template const imat elem_mult(const imat &m1, const imat &m2); 01716 // Extern Template Const instantiation of elem_mult 01717 //extern template const llmat elem_mult(const llmat &m1, const llmat &m2); 01719 extern template const smat elem_mult(const smat &m1, const smat &m2); 01721 extern template const bmat elem_mult(const bmat &m1, const bmat &m2); 01722 01723 // ------------ Division operator ----------- 01724 01726 extern template const mat operator/(const mat &m, double t); 01728 extern template const cmat operator/(const cmat &m, std::complex<double> t); 01730 extern template const imat operator/(const imat &m, int t); 01732 extern template const smat operator/(const smat &m, short t); 01734 extern template const bmat operator/(const bmat &m, bin t); 01735 01736 // ------------ Elementwise division ----------- 01737 01739 extern template const mat elem_div(const mat &m1, const mat &m2); 01741 extern template const cmat elem_div(const cmat &m1, const cmat &m2); 01743 extern template const imat elem_div(const imat &m1, const imat &m2); 01745 extern template const smat elem_div(const smat &m1, const smat &m2); 01747 extern template const bmat elem_div(const bmat &m1, const bmat &m2); 01748 01749 // ------------- Concatenations ----------------- 01750 01752 extern template const mat concat_horizontal(const mat &m1, const mat &m2); 01754 extern template const cmat concat_horizontal(const cmat &m1, const cmat &m2); 01756 extern template const imat concat_horizontal(const imat &m1, const imat &m2); 01758 extern template const smat concat_horizontal(const smat &m1, const smat &m2); 01760 extern template const bmat concat_horizontal(const bmat &m1, const bmat &m2); 01761 01763 extern template const mat concat_vertical(const mat &m1, const mat &m2); 01765 extern template const cmat concat_vertical(const cmat &m1, const cmat &m2); 01767 extern template const imat concat_vertical(const imat &m1, const imat &m2); 01769 extern template const smat concat_vertical(const smat &m1, const smat &m2); 01771 extern template const bmat concat_vertical(const bmat &m1, const bmat &m2); 01772 01773 //----------- Output stream -------------- 01774 01776 extern template std::ostream &operator<<(std::ostream &os, const mat &m); 01778 extern template std::ostream &operator<<(std::ostream &os, const cmat &m); 01780 extern template std::ostream &operator<<(std::ostream &os, const imat &m); 01782 extern template std::ostream &operator<<(std::ostream &os, const smat &m); 01784 extern template std::ostream &operator<<(std::ostream &os, const bmat &m); 01785 01786 //----------- Input stream ------------- 01787 01789 extern template std::istream &operator>>(std::istream &is, mat &m); 01791 extern template std::istream &operator>>(std::istream &is, cmat &m); 01793 extern template std::istream &operator>>(std::istream &is, imat &m); 01795 extern template std::istream &operator>>(std::istream &is, smat &m); 01797 extern template std::istream &operator>>(std::istream &is, bmat &m); 01798 01799 #endif // #ifndef _MSC_VER 01800 01801 } // namespace itpp 01802 01803 #endif // #ifndef MAT_H
Generated on Thu Apr 19 14:43:43 2007 for IT++ by Doxygen 1.5.1