IT++ Logo Newcom Logo

mat.h

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

Generated on Sat Aug 25 23:40:03 2007 for IT++ by Doxygen 1.5.2