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

Generated on Wed Jul 27 2011 16:27:04 for IT++ by Doxygen 1.7.4