mixpp: mat.h Source File

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