IT++ Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mat.h
Go to the documentation of this file.
1 
29 #ifndef MAT_H
30 #define MAT_H
31 
32 #include <itpp/base/itassert.h>
33 #include <itpp/base/math/misc.h>
34 #include <itpp/base/factory.h>
35 
36 
37 namespace itpp
38 {
39 
40 // Declaration of Vec
41 template<class Num_T> class Vec;
42 // Declaration of Mat
43 template<class Num_T> class Mat;
44 // Declaration of bin
45 class bin;
46 
48 template<class Num_T>
49 Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
51 template<class Num_T>
52 Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
53 
55 template<class Num_T>
56 Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
58 template<class Num_T>
59 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t);
61 template<class Num_T>
62 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m);
63 
65 template<class Num_T>
66 Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
68 template<class Num_T>
69 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t);
71 template<class Num_T>
72 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m);
74 template<class Num_T>
75 Mat<Num_T> operator-(const Mat<Num_T> &m);
76 
78 template<class Num_T>
79 Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
81 template<class Num_T>
82 Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v);
84 template<class Num_T>
85 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t);
87 template<class Num_T>
88 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m);
89 
91 template<class Num_T>
92 Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
94 template<class Num_T>
95 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
96  Mat<Num_T> &out);
98 template<class Num_T>
99 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
100  const Mat<Num_T> &m3, Mat<Num_T> &out);
102 template<class Num_T>
103 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
104  const Mat<Num_T> &m3, const Mat<Num_T> &m4,
105  Mat<Num_T> &out);
107 template<class Num_T>
108 void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2);
110 template<class Num_T>
111 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
112 
114 template<class Num_T>
115 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t);
117 template<class Num_T>
118 Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m);
119 
121 template<class Num_T>
122 Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
124 template<class Num_T>
125 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
126  Mat<Num_T> &out);
128 template<class Num_T>
129 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
130 
131 // -------------------------------------------------------------------------------------
132 // Declaration of Mat
133 // -------------------------------------------------------------------------------------
134 
200 template<class Num_T>
201 class Mat
202 {
203 public:
205  typedef Num_T value_type;
206 
208  explicit Mat(const Factory &f = DEFAULT_FACTORY);
210  Mat(int rows, int cols, const Factory &f = DEFAULT_FACTORY);
212  Mat(const Mat<Num_T> &m);
214  Mat(const Mat<Num_T> &m, const Factory &f);
216  Mat(const Vec<Num_T> &v, const Factory &f = DEFAULT_FACTORY);
218  Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY);
220  Mat(const char *str, const Factory &f = DEFAULT_FACTORY);
228  Mat(const Num_T *c_array, int rows, int cols, bool row_major = true,
229  const Factory &f = DEFAULT_FACTORY);
230 
232  ~Mat();
233 
235  int cols() const { return no_cols; }
237  int rows() const { return no_rows; }
239  int size() const { return datasize; }
241  void set_size(int rows, int cols, bool copy = false);
243  void zeros();
245  void clear() { zeros(); }
247  void ones();
249  void set(const std::string &str);
251  void set(const char *str);
252 
254  const Num_T &operator()(int r, int c) const;
256  Num_T &operator()(int r, int c);
258  const Num_T &operator()(int i) const;
260  Num_T &operator()(int i);
262  const Num_T &get(int r, int c) const;
264  const Num_T &get(int i) const;
266  void set(int r, int c, Num_T t);
267 
273  Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const;
279  Mat<Num_T> get(int r1, int r2, int c1, int c2) const;
280 
282  Vec<Num_T> get_row(int r) const;
284  Mat<Num_T> get_rows(int r1, int r2) const;
286  Mat<Num_T> get_rows(const Vec<int> &indexlist) const;
288  Vec<Num_T> get_col(int c) const;
290  Mat<Num_T> get_cols(int c1, int c2) const;
292  Mat<Num_T> get_cols(const Vec<int> &indexlist) const;
294  void set_row(int r, const Vec<Num_T> &v);
296  void set_col(int c, const Vec<Num_T> &v);
298  void set_rows(int r, const Mat<Num_T> &m);
300  void set_cols(int c, const Mat<Num_T> &m);
302  void copy_row(int to, int from);
304  void copy_col(int to, int from);
306  void swap_rows(int r1, int r2);
308  void swap_cols(int c1, int c2);
309 
311  void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m);
313  void set_submatrix(int r, int c, const Mat<Num_T> &m);
315  void set_submatrix(int r1, int r2, int c1, int c2, Num_T t);
316 
318  void del_row(int r);
320  void del_rows(int r1, int r2);
322  void del_col(int c);
324  void del_cols(int c1, int c2);
326  void ins_row(int r, const Vec<Num_T> &v);
328  void ins_col(int c, const Vec<Num_T> &v);
330  void append_row(const Vec<Num_T> &v);
332  void append_col(const Vec<Num_T> &v);
333 
335  Mat<Num_T> transpose() const;
337  Mat<Num_T> T() const { return this->transpose(); }
341  Mat<Num_T> H() const { return this->hermitian_transpose(); }
342 
344  friend Mat<Num_T> concat_horizontal<>(const Mat<Num_T> &m1,
345  const Mat<Num_T> &m2);
347  friend Mat<Num_T> concat_vertical<>(const Mat<Num_T> &m1,
348  const Mat<Num_T> &m2);
349 
351  Mat<Num_T>& operator=(Num_T t);
353  Mat<Num_T>& operator=(const Mat<Num_T> &m);
355  Mat<Num_T>& operator=(const Vec<Num_T> &v);
357  Mat<Num_T>& operator=(const std::string &str);
359  Mat<Num_T>& operator=(const char *str);
360 
362  Mat<Num_T>& operator+=(const Mat<Num_T> &m);
364  Mat<Num_T>& operator+=(Num_T t);
366  friend Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
368  friend Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t);
370  friend Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m);
371 
373  Mat<Num_T>& operator-=(const Mat<Num_T> &m);
375  Mat<Num_T>& operator-=(Num_T t);
377  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
379  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t);
381  friend Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m);
383  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m);
384 
386  Mat<Num_T>& operator*=(const Mat<Num_T> &m);
388  Mat<Num_T>& operator*=(Num_T t);
390  friend Mat<Num_T> operator*<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
392  friend Vec<Num_T> operator*<>(const Mat<Num_T> &m, const Vec<Num_T> &v);
394  friend Mat<Num_T> operator*<>(const Mat<Num_T> &m, Num_T t);
396  friend Mat<Num_T> operator*<>(Num_T t, const Mat<Num_T> &m);
397 
399  friend Mat<Num_T> elem_mult<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
401  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
402  Mat<Num_T> &out);
404  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
405  const Mat<Num_T> &m3, Mat<Num_T> &out);
407  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
408  const Mat<Num_T> &m3, const Mat<Num_T> &m4,
409  Mat<Num_T> &out);
411  friend void elem_mult_inplace<>(const Mat<Num_T> &m1, Mat<Num_T> &m2);
413  friend Num_T elem_mult_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
414 
416  Mat<Num_T>& operator/=(Num_T t);
418  Mat<Num_T>& operator/=(const Mat<Num_T> &m);
419 
421  friend Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t);
423  friend Mat<Num_T> operator/<>(Num_T t, const Mat<Num_T> &m);
424 
426  friend Mat<Num_T> elem_div<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
428  friend void elem_div_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
429  Mat<Num_T> &out);
431  friend Num_T elem_div_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
432 
434  bool operator==(const Mat<Num_T> &m) const;
436  bool operator!=(const Mat<Num_T> &m) const;
437 
439  Num_T &_elem(int r, int c) { return data[r+c*no_rows]; }
441  const Num_T &_elem(int r, int c) const { return data[r+c*no_rows]; }
443  Num_T &_elem(int i) { return data[i]; }
445  const Num_T &_elem(int i) const { return data[i]; }
446 
448  Num_T *_data() { return data; }
450  const Num_T *_data() const { return data; }
452  int _datasize() const { return datasize; }
453 
454 protected:
456  void alloc(int rows, int cols);
458  void free();
459 
464 
465  Num_T *data;
467  const Factory &factory;
468 
469 private:
471  bool in_range(int r, int c) const {
472  return ((r >= 0) && (r < no_rows) && (c >= 0) && (c < no_cols));
473  }
475  bool row_in_range(int r) const { return ((r >= 0) && (r < no_rows)); }
477  bool col_in_range(int c) const { return ((c >= 0) && (c < no_cols)); }
479  bool in_range(int i) const { return ((i >= 0) && (i < datasize)); }
480 };
481 
482 // -------------------------------------------------------------------------------------
483 // Type definitions of mat, cmat, imat, smat, and bmat
484 // -------------------------------------------------------------------------------------
485 
490 typedef Mat<double> mat;
491 
497 
502 typedef Mat<int> imat;
503 
509 
516 typedef Mat<bin> bmat;
517 
518 } //namespace itpp
519 
520 
521 #include <itpp/base/vec.h>
522 
523 namespace itpp
524 {
525 
526 // ----------------------------------------------------------------------
527 // Declaration of input and output streams for Mat
528 // ----------------------------------------------------------------------
529 
534 template <class Num_T>
535 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m);
536 
548 template <class Num_T>
549 std::istream &operator>>(std::istream &is, Mat<Num_T> &m);
550 
551 // ----------------------------------------------------------------------
552 // Implementation of templated Mat members and friends
553 // ----------------------------------------------------------------------
554 
555 template<class Num_T> inline
556 void Mat<Num_T>::alloc(int rows, int cols)
557 {
558  if ((rows > 0) && (cols > 0)) {
559  datasize = rows * cols;
560  no_rows = rows;
561  no_cols = cols;
562  create_elements(data, datasize, factory);
563  }
564  else {
565  data = 0;
566  datasize = 0;
567  no_rows = 0;
568  no_cols = 0;
569  }
570 }
571 
572 template<class Num_T> inline
574 {
575  destroy_elements(data, datasize);
576  datasize = 0;
577  no_rows = 0;
578  no_cols = 0;
579 }
580 
581 
582 template<class Num_T> inline
584  datasize(0), no_rows(0), no_cols(0), data(0), factory(f) {}
585 
586 template<class Num_T> inline
587 Mat<Num_T>::Mat(int rows, int cols, const Factory &f) :
588  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
589 {
590  it_assert_debug((rows >= 0) && (cols >= 0), "Mat<>::Mat(): Wrong size");
591  alloc(rows, cols);
592 }
593 
594 template<class Num_T> inline
596  datasize(0), no_rows(0), no_cols(0), data(0), factory(m.factory)
597 {
598  alloc(m.no_rows, m.no_cols);
599  copy_vector(m.datasize, m.data, data);
600 }
601 
602 template<class Num_T> inline
603 Mat<Num_T>::Mat(const Mat<Num_T> &m, const Factory &f) :
604  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
605 {
606  alloc(m.no_rows, m.no_cols);
607  copy_vector(m.datasize, m.data, data);
608 }
609 
610 template<class Num_T> inline
611 Mat<Num_T>::Mat(const Vec<Num_T> &v, const Factory &f) :
612  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
613 {
614  int size = v.size();
615  alloc(size, 1);
616  copy_vector(size, v._data(), data);
617 }
618 
619 template<class Num_T> inline
620 Mat<Num_T>::Mat(const std::string &str, const Factory &f) :
621  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
622 {
623  set(str);
624 }
625 
626 template<class Num_T> inline
627 Mat<Num_T>::Mat(const char *str, const Factory &f) :
628  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
629 {
630  set(std::string(str));
631 }
632 
633 template<class Num_T>
634 Mat<Num_T>::Mat(const Num_T *c_array, int rows, int cols, bool row_major,
635  const Factory &f):
636  datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
637 {
638  alloc(rows, cols);
639  if (!row_major)
640  copy_vector(datasize, c_array, data);
641  else
642  for (int i = 0; i < rows; i++)
643  for (int j = 0; j < cols; j++)
644  data[i+j*no_rows] = c_array[i*no_cols+j];
645 }
646 
647 template<class Num_T> inline
649 {
650  free();
651 }
652 
653 
654 template<class Num_T>
655 void Mat<Num_T>::set_size(int rows, int cols, bool copy)
656 {
657  it_assert_debug((rows >= 0) && (cols >= 0),
658  "Mat<>::set_size(): Wrong size");
659  // check if we have to resize the current matrix
660  if ((no_rows == rows) && (no_cols == cols))
661  return;
662  // check if one of dimensions is zero
663  if ((rows == 0) || (cols == 0)) {
664  free();
665  return;
666  }
667  // conditionally copy previous matrix content
668  if (copy) {
669  // create a temporary pointer to the allocated data
670  Num_T* tmp = data;
671  // store the current number of elements and number of rows
672  int old_datasize = datasize;
673  int old_rows = no_rows;
674  // check the boundaries of the copied data
675  int min_r = (no_rows < rows) ? no_rows : rows;
676  int min_c = (no_cols < cols) ? no_cols : cols;
677  // allocate new memory
678  alloc(rows, cols);
679  // copy the previous data into the allocated memory
680  for (int i = 0; i < min_c; ++i) {
681  copy_vector(min_r, &tmp[i*old_rows], &data[i*no_rows]);
682  }
683  // fill-in the rest of matrix with zeros
684  for (int i = min_r; i < rows; ++i)
685  for (int j = 0; j < cols; ++j)
686  data[i+j*rows] = Num_T(0);
687  for (int j = min_c; j < cols; ++j)
688  for (int i = 0; i < min_r; ++i)
689  data[i+j*rows] = Num_T(0);
690  // delete old elements
691  destroy_elements(tmp, old_datasize);
692  }
693  // if possible, reuse the allocated memory
694  else if (datasize == rows * cols) {
695  no_rows = rows;
696  no_cols = cols;
697  }
698  // finally release old memory and allocate a new one
699  else {
700  free();
701  alloc(rows, cols);
702  }
703 }
704 
705 template<class Num_T> inline
707 {
708  for (int i = 0; i < datasize; i++)
709  data[i] = Num_T(0);
710 }
711 
712 template<class Num_T> inline
714 {
715  for (int i = 0; i < datasize; i++)
716  data[i] = Num_T(1);
717 }
718 
719 template<class Num_T> inline
720 const Num_T& Mat<Num_T>::operator()(int r, int c) const
721 {
722  it_assert_debug(in_range(r, c),
723  "Mat<>::operator(): Indexing out of range");
724  return data[r+c*no_rows];
725 }
726 
727 template<class Num_T> inline
728 Num_T& Mat<Num_T>::operator()(int r, int c)
729 {
730  it_assert_debug(in_range(r, c),
731  "Mat<>::operator(): Indexing out of range");
732  return data[r+c*no_rows];
733 }
734 
735 template<class Num_T> inline
737 {
738  it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range");
739  return data[i];
740 }
741 
742 template<class Num_T> inline
743 const Num_T& Mat<Num_T>::operator()(int i) const
744 {
745  it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range");
746  return data[i];
747 }
748 
749 template<class Num_T> inline
750 const Num_T& Mat<Num_T>::get(int r, int c) const
751 {
752  return (*this)(r, c);
753 }
754 
755 template<class Num_T> inline
756 const Num_T& Mat<Num_T>::get(int i) const
757 {
758  return (*this)(i);
759 }
760 
761 template<class Num_T> inline
762 void Mat<Num_T>::set(int r, int c, Num_T t)
763 {
764  it_assert_debug(in_range(r, c), "Mat<>::set(): Indexing out of range");
765  data[r+c*no_rows] = t;
766 }
767 
768 
769 template<class Num_T>
770 void Mat<Num_T>::set(const std::string &str)
771 {
772  // actual row counter
773  int rows = 0;
774  // number of rows to allocate next time (8, 16, 32, 64, etc.)
775  int maxrows = 8;
776 
777  // clean the current matrix content
778  free();
779 
780  // variable to store the start of a current vector
781  std::string::size_type beg = 0;
782  std::string::size_type end = 0;
783  while (end != std::string::npos) {
784  // find next occurrence of a semicolon in string str
785  end = str.find(';', beg);
786  // parse first row into a vector v
787  Vec<Num_T> v(str.substr(beg, end - beg));
788  int v_size = v.size();
789 
790  // this check is necessary to parse the following two strings as the
791  // same matrix: "1 0 1; ; 1 1; " and "1 0 1; 0 0 0; 1 1 0"
792  if ((end != std::string::npos) || (v_size > 0)) {
793  // matrix empty -> insert v as a first row and allocate maxrows
794  if (rows == 0) {
795  set_size(maxrows, v_size, true);
796  set_row(rows++, v);
797  }
798  else {
799  // check if we need to resize the matrix
800  if ((rows == maxrows) || (v_size != no_cols)) {
801  // we need to add new rows
802  if (rows == maxrows) {
803  maxrows *= 2;
804  }
805  // check if we need to add new columns
806  if (v_size > no_cols) {
807  set_size(maxrows, v_size, true);
808  }
809  else {
810  set_size(maxrows, no_cols, true);
811  // set the size of the parsed vector to the number of columns
812  v.set_size(no_cols, true);
813  }
814  }
815  // set the parsed vector as the next row
816  set_row(rows++, v);
817  }
818  }
819  // update the starting position of the next vector in the parsed
820  // string
821  beg = end + 1;
822  } // if ((end != std::string::npos) || (v.size > 0))
823 
824  set_size(rows, no_cols, true);
825 }
826 
827 template<class Num_T> inline
828 void Mat<Num_T>::set(const char *str)
829 {
830  set(std::string(str));
831 }
832 
833 template<class Num_T> inline
834 Mat<Num_T> Mat<Num_T>::operator()(int r1, int r2, int c1, int c2) const
835 {
836  if (r1 == -1) r1 = no_rows - 1;
837  if (r2 == -1) r2 = no_rows - 1;
838  if (c1 == -1) c1 = no_cols - 1;
839  if (c2 == -1) c2 = no_cols - 1;
840 
841  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) &&
842  (c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
843  "Mat<>::operator()(r1, r2, c1, c2): Wrong indexing");
844 
845  Mat<Num_T> s(r2 - r1 + 1, c2 - c1 + 1);
846 
847  for (int i = 0;i < s.no_cols;i++)
848  copy_vector(s.no_rows, data + r1 + (c1 + i)*no_rows, s.data + i*s.no_rows);
849 
850  return s;
851 }
852 
853 template<class Num_T> inline
854 Mat<Num_T> Mat<Num_T>::get(int r1, int r2, int c1, int c2) const
855 {
856  return (*this)(r1, r2, c1, c2);
857 }
858 
859 template<class Num_T> inline
861 {
862  it_assert_debug(row_in_range(r), "Mat<>::get_row(): Index out of range");
863  Vec<Num_T> a(no_cols);
864 
865  copy_vector(no_cols, data + r, no_rows, a._data(), 1);
866  return a;
867 }
868 
869 template<class Num_T>
870 Mat<Num_T> Mat<Num_T>::get_rows(int r1, int r2) const
871 {
872  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows),
873  "Mat<>::get_rows(): Wrong indexing");
874  Mat<Num_T> m(r2 - r1 + 1, no_cols);
875 
876  for (int i = 0; i < m.rows(); i++)
877  copy_vector(no_cols, data + i + r1, no_rows, m.data + i, m.no_rows);
878 
879  return m;
880 }
881 
882 template<class Num_T>
883 Mat<Num_T> Mat<Num_T>::get_rows(const Vec<int> &indexlist) const
884 {
885  Mat<Num_T> m(indexlist.size(), no_cols);
886 
887  for (int i = 0;i < indexlist.size();i++) {
888  it_assert_debug(row_in_range(indexlist(i)),
889  "Mat<>::get_rows(indexlist): Indexing out of range");
890  copy_vector(no_cols, data + indexlist(i), no_rows, m.data + i, m.no_rows);
891  }
892 
893  return m;
894 }
895 
896 template<class Num_T> inline
898 {
899  it_assert_debug(col_in_range(c), "Mat<>::get_col(): Index out of range");
900  Vec<Num_T> a(no_rows);
901 
902  copy_vector(no_rows, data + c*no_rows, a._data());
903 
904  return a;
905 }
906 
907 template<class Num_T>
908 Mat<Num_T> Mat<Num_T>::get_cols(int c1, int c2) const
909 {
910  it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
911  "Mat<>::get_cols(): Wrong indexing");
912  Mat<Num_T> m(no_rows, c2 - c1 + 1);
913 
914  for (int i = 0; i < m.cols(); i++)
915  copy_vector(no_rows, data + (i + c1)*no_rows, m.data + i*m.no_rows);
916 
917  return m;
918 }
919 
920 template<class Num_T>
921 Mat<Num_T> Mat<Num_T>::get_cols(const Vec<int> &indexlist) const
922 {
923  Mat<Num_T> m(no_rows, indexlist.size());
924 
925  for (int i = 0; i < indexlist.size(); i++) {
926  it_assert_debug(col_in_range(indexlist(i)),
927  "Mat<>::get_cols(indexlist): Indexing out of range");
928  copy_vector(no_rows, data + indexlist(i)*no_rows, m.data + i*m.no_rows);
929  }
930 
931  return m;
932 }
933 
934 template<class Num_T> inline
935 void Mat<Num_T>::set_row(int r, const Vec<Num_T> &v)
936 {
937  it_assert_debug(row_in_range(r), "Mat<>::set_row(): Index out of range");
938  it_assert_debug(v.size() == no_cols,
939  "Mat<>::set_row(): Wrong size of input vector");
940  copy_vector(v.size(), v._data(), 1, data + r, no_rows);
941 }
942 
943 template<class Num_T> inline
944 void Mat<Num_T>::set_col(int c, const Vec<Num_T> &v)
945 {
946  it_assert_debug(col_in_range(c), "Mat<>::set_col(): Index out of range");
947  it_assert_debug(v.size() == no_rows,
948  "Mat<>::set_col(): Wrong size of input vector");
949  copy_vector(v.size(), v._data(), data + c*no_rows);
950 }
951 
952 
953 template<class Num_T>
954 void Mat<Num_T>::set_rows(int r, const Mat<Num_T> &m)
955 {
956  it_assert_debug(row_in_range(r), "Mat<>::set_rows(): Index out of range");
957  it_assert_debug(no_cols == m.cols(),
958  "Mat<>::set_rows(): Column sizes do not match");
959  it_assert_debug(m.rows() + r <= no_rows,
960  "Mat<>::set_rows(): Not enough rows");
961 
962  for (int i = 0; i < m.rows(); ++i) {
963  copy_vector(no_cols, m.data + i, m.no_rows, data + i + r, no_rows);
964  }
965 }
966 
967 template<class Num_T>
968 void Mat<Num_T>::set_cols(int c, const Mat<Num_T> &m)
969 {
970  it_assert_debug(col_in_range(c), "Mat<>::set_cols(): Index out of range");
971  it_assert_debug(no_rows == m.rows(),
972  "Mat<>::set_cols(): Row sizes do not match");
973  it_assert_debug(m.cols() + c <= no_cols,
974  "Mat<>::set_cols(): Not enough colums");
975 
976  for (int i = 0; i < m.cols(); ++i) {
977  copy_vector(no_rows, m.data + i*no_rows, data + (i + c)*no_rows);
978  }
979 }
980 
981 
982 template<class Num_T> inline
983 void Mat<Num_T>::copy_row(int to, int from)
984 {
985  it_assert_debug(row_in_range(to) && row_in_range(from),
986  "Mat<>::copy_row(): Indexing out of range");
987  if (from == to)
988  return;
989 
990  copy_vector(no_cols, data + from, no_rows, data + to, no_rows);
991 }
992 
993 template<class Num_T> inline
994 void Mat<Num_T>::copy_col(int to, int from)
995 {
996  it_assert_debug(col_in_range(to) && col_in_range(from),
997  "Mat<>::copy_col(): Indexing out of range");
998  if (from == to)
999  return;
1000 
1001  copy_vector(no_rows, data + from*no_rows, data + to*no_rows);
1002 }
1003 
1004 template<class Num_T> inline
1005 void Mat<Num_T>::swap_rows(int r1, int r2)
1006 {
1007  it_assert_debug(row_in_range(r1) && row_in_range(r2),
1008  "Mat<>::swap_rows(): Indexing out of range");
1009  if (r1 == r2)
1010  return;
1011 
1012  swap_vector(no_cols, data + r1, no_rows, data + r2, no_rows);
1013 }
1014 
1015 template<class Num_T> inline
1016 void Mat<Num_T>::swap_cols(int c1, int c2)
1017 {
1018  it_assert_debug(col_in_range(c1) && col_in_range(c2),
1019  "Mat<>::swap_cols(): Indexing out of range");
1020  if (c1 == c2)
1021  return;
1022 
1023  swap_vector(no_rows, data + c1*no_rows, data + c2*no_rows);
1024 }
1025 
1026 template<class Num_T>
1027 void Mat<Num_T>::set_submatrix(int r1, int, int c1, int, const Mat<Num_T> &m)
1028 {
1029  it_warning("Mat<>::set_submatrix(r1, r2, r3, r4, m): This function is "
1030  "deprecated and might be removed from future IT++ releases. "
1031  "Please use Mat<>::set_submatrix(r, c, m) function instead.");
1032  set_submatrix(r1, c1, m);
1033 }
1034 
1035 template<class Num_T> inline
1036 void Mat<Num_T>::set_submatrix(int r, int c, const Mat<Num_T> &m)
1037 {
1038  it_assert_debug((r >= 0) && (r + m.no_rows <= no_rows) &&
1039  (c >= 0) && (c + m.no_cols <= no_cols),
1040  "Mat<>::set_submatrix(): Indexing out of range "
1041  "or wrong input matrix");
1042  for (int i = 0; i < m.no_cols; i++)
1043  copy_vector(m.no_rows, m.data + i*m.no_rows, data + (c + i)*no_rows + r);
1044 }
1045 
1046 
1047 
1048 template<class Num_T> inline
1049 void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, Num_T t)
1050 {
1051  if (r1 == -1) r1 = no_rows - 1;
1052  if (r2 == -1) r2 = no_rows - 1;
1053  if (c1 == -1) c1 = no_cols - 1;
1054  if (c2 == -1) c2 = no_cols - 1;
1055  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) &&
1056  (c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
1057  "Mat<>::set_submatrix(): Wrong indexing");
1058  for (int i = c1; i <= c2; i++) {
1059  int pos = i * no_rows + r1;
1060  for (int j = r1; j <= r2; j++)
1061  data[pos++] = t;
1062  }
1063 }
1064 
1065 template<class Num_T>
1067 {
1068  it_assert_debug(row_in_range(r), "Mat<>::del_row(): Index out of range");
1069  Mat<Num_T> Temp(*this);
1070  set_size(no_rows - 1, no_cols, false);
1071  for (int i = 0 ; i < r ; i++) {
1072  copy_vector(no_cols, &Temp.data[i], no_rows + 1, &data[i], no_rows);
1073  }
1074  for (int i = r ; i < no_rows ; i++) {
1075  copy_vector(no_cols, &Temp.data[i+1], no_rows + 1, &data[i], no_rows);
1076  }
1077 
1078 }
1079 
1080 template<class Num_T>
1081 void Mat<Num_T>::del_rows(int r1, int r2)
1082 {
1083  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows),
1084  "Mat<>::del_rows(): Indexing out of range");
1085  Mat<Num_T> Temp(*this);
1086  int no_del_rows = r2 - r1 + 1;
1087  set_size(no_rows - no_del_rows, no_cols, false);
1088  for (int i = 0; i < r1 ; ++i) {
1089  copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows);
1090  }
1091  for (int i = r2 + 1; i < Temp.no_rows; ++i) {
1092  copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i-no_del_rows],
1093  no_rows);
1094  }
1095 }
1096 
1097 template<class Num_T>
1099 {
1100  it_assert_debug(col_in_range(c), "Mat<>::del_col(): Index out of range");
1101  Mat<Num_T> Temp(*this);
1102 
1103  set_size(no_rows, no_cols - 1, false);
1104  copy_vector(c*no_rows, Temp.data, data);
1105  copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]);
1106 }
1107 
1108 template<class Num_T>
1109 void Mat<Num_T>::del_cols(int c1, int c2)
1110 {
1111  it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
1112  "Mat<>::del_cols(): Indexing out of range");
1113  Mat<Num_T> Temp(*this);
1114  int n_deleted_cols = c2 - c1 + 1;
1115  set_size(no_rows, no_cols - n_deleted_cols, false);
1116  copy_vector(c1*no_rows, Temp.data, data);
1117  copy_vector((no_cols - c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]);
1118 }
1119 
1120 template<class Num_T>
1121 void Mat<Num_T>::ins_row(int r, const Vec<Num_T> &v)
1122 {
1123  it_assert_debug((r >= 0) && (r <= no_rows),
1124  "Mat<>::ins_row(): Index out of range");
1125  it_assert_debug((v.size() == no_cols) || (no_rows == 0),
1126  "Mat<>::ins_row(): Wrong size of the input vector");
1127 
1128  if (no_cols == 0) {
1129  no_cols = v.size();
1130  }
1131 
1132  Mat<Num_T> Temp(*this);
1133  set_size(no_rows + 1, no_cols, false);
1134 
1135  for (int i = 0 ; i < r ; i++) {
1136  copy_vector(no_cols, &Temp.data[i], no_rows - 1, &data[i], no_rows);
1137  }
1138  copy_vector(no_cols, v._data(), 1, &data[r], no_rows);
1139  for (int i = r + 1 ; i < no_rows ; i++) {
1140  copy_vector(no_cols, &Temp.data[i-1], no_rows - 1, &data[i], no_rows);
1141  }
1142 }
1143 
1144 template<class Num_T>
1145 void Mat<Num_T>::ins_col(int c, const Vec<Num_T> &v)
1146 {
1147  it_assert_debug((c >= 0) && (c <= no_cols),
1148  "Mat<>::ins_col(): Index out of range");
1149  it_assert_debug((v.size() == no_rows) || (no_cols == 0),
1150  "Mat<>::ins_col(): Wrong size of the input vector");
1151 
1152  if (no_rows == 0) {
1153  no_rows = v.size();
1154  }
1155 
1156  Mat<Num_T> Temp(*this);
1157  set_size(no_rows, no_cols + 1, false);
1158 
1159  copy_vector(c*no_rows, Temp.data, data);
1160  copy_vector(no_rows, v._data(), &data[c*no_rows]);
1161  copy_vector((no_cols - c - 1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]);
1162 }
1163 
1164 template<class Num_T> inline
1166 {
1167  ins_row(no_rows, v);
1168 }
1169 
1170 template<class Num_T> inline
1172 {
1173  ins_col(no_cols, v);
1174 }
1175 
1176 template<class Num_T>
1178 {
1179  Mat<Num_T> temp(no_cols, no_rows);
1180  for (int i = 0; i < no_rows; ++i) {
1181  copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1);
1182  }
1183  return temp;
1184 }
1185 
1187 template<>
1188 cmat cmat::hermitian_transpose() const;
1190 
1191 template<class Num_T>
1193 {
1194  Mat<Num_T> temp(no_cols, no_rows);
1195  for (int i = 0; i < no_rows; ++i) {
1196  copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1);
1197  }
1198  return temp;
1199 }
1200 
1201 template<class Num_T>
1203 {
1204  // if one of the input matrix is empty just copy the other one as a result
1205  if (m1.no_cols == 0)
1206  return m2;
1207  if (m2.no_cols == 0)
1208  return m1;
1209  it_assert_debug(m1.no_rows == m2.no_rows,
1210  "Mat<>::concat_horizontal(): Wrong sizes");
1211  int no_rows = m1.no_rows;
1212  Mat<Num_T> temp(no_rows, m1.no_cols + m2.no_cols);
1213  for (int i = 0; i < m1.no_cols; ++i) {
1214  copy_vector(no_rows, &m1.data[i * no_rows], &temp.data[i * no_rows]);
1215  }
1216  for (int i = 0; i < m2.no_cols; ++i) {
1217  copy_vector(no_rows, &m2.data[i * no_rows], &temp.data[(m1.no_cols + i)
1218  * no_rows]);
1219  }
1220  return temp;
1221 }
1222 
1223 template<class Num_T>
1225 {
1226  // if one of the input matrix is empty just copy the other one as a result
1227  if (m1.no_rows == 0)
1228  return m2;
1229  if (m2.no_rows == 0)
1230  return m1;
1231  it_assert_debug(m1.no_cols == m2.no_cols,
1232  "Mat<>::concat_vertical(): Wrong sizes");
1233  int no_cols = m1.no_cols;
1234  Mat<Num_T> temp(m1.no_rows + m2.no_rows, no_cols);
1235  for (int i = 0; i < no_cols; ++i) {
1236  copy_vector(m1.no_rows, &m1.data[i * m1.no_rows],
1237  &temp.data[i * temp.no_rows]);
1238  copy_vector(m2.no_rows, &m2.data[i * m2.no_rows],
1239  &temp.data[i * temp.no_rows + m1.no_rows]);
1240  }
1241  return temp;
1242 }
1243 
1244 template<class Num_T> inline
1246 {
1247  for (int i = 0; i < datasize; i++)
1248  data[i] = t;
1249  return *this;
1250 }
1251 
1252 template<class Num_T> inline
1254 {
1255  if (this != &m) {
1256  set_size(m.no_rows, m.no_cols, false);
1257  if (m.datasize != 0)
1258  copy_vector(m.datasize, m.data, data);
1259  }
1260  return *this;
1261 }
1262 
1263 template<class Num_T> inline
1265 {
1266  it_assert_debug(((no_rows == 1) && (no_cols == v.size()))
1267  || ((no_cols == 1) && (no_rows == v.size())),
1268  "Mat<>::operator=(): Wrong size of the input vector");
1269  set_size(v.size(), 1, false);
1270  copy_vector(v.size(), v._data(), data);
1271  return *this;
1272 }
1273 
1274 template<class Num_T> inline
1275 Mat<Num_T>& Mat<Num_T>::operator=(const std::string &str)
1276 {
1277  set(str);
1278  return *this;
1279 }
1280 
1281 template<class Num_T> inline
1283 {
1284  set(std::string(str));
1285  return *this;
1286 }
1287 
1288 //-------------------- Templated friend functions --------------------------
1289 
1290 template<class Num_T>
1292 {
1293  if (datasize == 0)
1294  operator=(m);
1295  else {
1296  int i, j, m_pos = 0, pos = 0;
1297  it_assert_debug(m.no_rows == no_rows && m.no_cols == no_cols, "Mat<Num_T>::operator+=: wrong sizes");
1298  for (i = 0; i < no_cols; i++) {
1299  for (j = 0; j < no_rows; j++)
1300  data[pos+j] += m.data[m_pos+j];
1301  pos += no_rows;
1302  m_pos += m.no_rows;
1303  }
1304  }
1305  return *this;
1306 }
1307 
1308 template<class Num_T> inline
1310 {
1311  for (int i = 0; i < datasize; i++)
1312  data[i] += t;
1313  return *this;
1314 }
1315 
1316 template<class Num_T>
1318 {
1319  Mat<Num_T> r(m1.no_rows, m1.no_cols);
1320  int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0;
1321 
1322  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1323  "Mat<>::operator+(): Wrong sizes");
1324 
1325  for (i = 0; i < r.no_cols; i++) {
1326  for (j = 0; j < r.no_rows; j++)
1327  r.data[r_pos+j] = m1.data[m1_pos+j] + m2.data[m2_pos+j];
1328  // next column
1329  m1_pos += m1.no_rows;
1330  m2_pos += m2.no_rows;
1331  r_pos += r.no_rows;
1332  }
1333 
1334  return r;
1335 }
1336 
1337 
1338 template<class Num_T>
1339 Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t)
1340 {
1341  Mat<Num_T> r(m.no_rows, m.no_cols);
1342 
1343  for (int i = 0; i < r.datasize; i++)
1344  r.data[i] = m.data[i] + t;
1345 
1346  return r;
1347 }
1348 
1349 template<class Num_T>
1350 Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m)
1351 {
1352  Mat<Num_T> r(m.no_rows, m.no_cols);
1353 
1354  for (int i = 0; i < r.datasize; i++)
1355  r.data[i] = t + m.data[i];
1356 
1357  return r;
1358 }
1359 
1360 template<class Num_T>
1362 {
1363  int i, j, m_pos = 0, pos = 0;
1364 
1365  if (datasize == 0) {
1366  set_size(m.no_rows, m.no_cols, false);
1367  for (i = 0; i < no_cols; i++) {
1368  for (j = 0; j < no_rows; j++)
1369  data[pos+j] = -m.data[m_pos+j];
1370  // next column
1371  m_pos += m.no_rows;
1372  pos += no_rows;
1373  }
1374  }
1375  else {
1376  it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols),
1377  "Mat<>::operator-=(): Wrong sizes");
1378  for (i = 0; i < no_cols; i++) {
1379  for (j = 0; j < no_rows; j++)
1380  data[pos+j] -= m.data[m_pos+j];
1381  // next column
1382  m_pos += m.no_rows;
1383  pos += no_rows;
1384  }
1385  }
1386  return *this;
1387 }
1388 
1389 template<class Num_T>
1391 {
1392  Mat<Num_T> r(m1.no_rows, m1.no_cols);
1393  int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0;
1394  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1395  "Mat<>::operator-(): Wrong sizes");
1396 
1397  for (i = 0; i < r.no_cols; i++) {
1398  for (j = 0; j < r.no_rows; j++)
1399  r.data[r_pos+j] = m1.data[m1_pos+j] - m2.data[m2_pos+j];
1400  // next column
1401  m1_pos += m1.no_rows;
1402  m2_pos += m2.no_rows;
1403  r_pos += r.no_rows;
1404  }
1405 
1406  return r;
1407 }
1408 
1409 template<class Num_T> inline
1411 {
1412  for (int i = 0; i < datasize; i++)
1413  data[i] -= t;
1414  return *this;
1415 }
1416 
1417 template<class Num_T>
1418 Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t)
1419 {
1420  Mat<Num_T> r(m.no_rows, m.no_cols);
1421  int i, j, m_pos = 0, r_pos = 0;
1422 
1423  for (i = 0; i < r.no_cols; i++) {
1424  for (j = 0; j < r.no_rows; j++)
1425  r.data[r_pos+j] = m.data[m_pos+j] - t;
1426  // next column
1427  m_pos += m.no_rows;
1428  r_pos += r.no_rows;
1429  }
1430 
1431  return r;
1432 }
1433 
1434 template<class Num_T>
1435 Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m)
1436 {
1437  Mat<Num_T> r(m.no_rows, m.no_cols);
1438  int i, j, m_pos = 0, r_pos = 0;
1439 
1440  for (i = 0; i < r.no_cols; i++) {
1441  for (j = 0; j < r.no_rows; j++)
1442  r.data[r_pos+j] = t - m.data[m_pos+j];
1443  // next column
1444  m_pos += m.no_rows;
1445  r_pos += r.no_rows;
1446  }
1447 
1448  return r;
1449 }
1450 
1451 template<class Num_T>
1453 {
1454  Mat<Num_T> r(m.no_rows, m.no_cols);
1455  int i, j, m_pos = 0, r_pos = 0;
1456 
1457  for (i = 0; i < r.no_cols; i++) {
1458  for (j = 0; j < r.no_rows; j++)
1459  r.data[r_pos+j] = -m.data[m_pos+j];
1460  // next column
1461  m_pos += m.no_rows;
1462  r_pos += r.no_rows;
1463  }
1464 
1465  return r;
1466 }
1467 
1468 
1469 template<> mat& mat::operator*=(const mat &m);
1470 template<> cmat& cmat::operator*=(const cmat &m);
1471 
1472 template<class Num_T>
1474 {
1475  it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
1476  Mat<Num_T> r(no_rows, m.no_cols);
1477 
1478  Num_T tmp;
1479 
1480  int i, j, k, r_pos = 0, pos = 0, m_pos = 0;
1481 
1482  for (i = 0; i < r.no_cols; i++) {
1483  for (j = 0; j < r.no_rows; j++) {
1484  tmp = Num_T(0);
1485  pos = 0;
1486  for (k = 0; k < no_cols; k++) {
1487  tmp += data[pos+j] * m.data[m_pos+k];
1488  pos += no_rows;
1489  }
1490  r.data[r_pos+j] = tmp;
1491  }
1492  r_pos += r.no_rows;
1493  m_pos += m.no_rows;
1494  }
1495  operator=(r); // time consuming
1496  return *this;
1497 }
1498 
1499 template<class Num_T> inline
1501 {
1502  scal_vector(datasize, t, data);
1503  return *this;
1504 }
1505 
1506 
1507 template<> mat operator*(const mat &m1, const mat &m2);
1508 template<> cmat operator*(const cmat &m1, const cmat &m2);
1509 
1510 template<class Num_T>
1512 {
1513  it_assert_debug(m1.no_cols == m2.no_rows,
1514  "Mat<>::operator*(): Wrong sizes");
1515  Mat<Num_T> r(m1.no_rows, m2.no_cols);
1516 
1517  Num_T tmp;
1518  int i, j, k;
1519  Num_T *tr = r.data, *t1, *t2 = m2.data;
1520 
1521  for (i = 0; i < r.no_cols; i++) {
1522  for (j = 0; j < r.no_rows; j++) {
1523  tmp = Num_T(0);
1524  t1 = m1.data + j;
1525  for (k = m1.no_cols; k > 0; k--) {
1526  tmp += *(t1) * *(t2++);
1527  t1 += m1.no_rows;
1528  }
1529  *(tr++) = tmp;
1530  t2 -= m2.no_rows;
1531  }
1532  t2 += m2.no_rows;
1533  }
1534 
1535  return r;
1536 }
1537 
1538 
1539 template<> vec operator*(const mat &m, const vec &v);
1540 template<> cvec operator*(const cmat &m, const cvec &v);
1541 
1542 template<class Num_T>
1544 {
1545  it_assert_debug(m.no_cols == v.size(),
1546  "Mat<>::operator*(): Wrong sizes");
1547  Vec<Num_T> r(m.no_rows);
1548  int i, k, m_pos;
1549 
1550  for (i = 0; i < m.no_rows; i++) {
1551  r(i) = Num_T(0);
1552  m_pos = 0;
1553  for (k = 0; k < m.no_cols; k++) {
1554  r(i) += m.data[m_pos+i] * v(k);
1555  m_pos += m.no_rows;
1556  }
1557  }
1558 
1559  return r;
1560 }
1561 
1562 template<class Num_T>
1563 Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t)
1564 {
1565  Mat<Num_T> r(m.no_rows, m.no_cols);
1566 
1567  for (int i = 0; i < r.datasize; i++)
1568  r.data[i] = m.data[i] * t;
1569 
1570  return r;
1571 }
1572 
1573 template<class Num_T> inline
1574 Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m)
1575 {
1576  return operator*(m, t);
1577 }
1578 
1579 template<class Num_T> inline
1581 {
1582  Mat<Num_T> out;
1583  elem_mult_out(m1, m2, out);
1584  return out;
1585 }
1586 
1587 template<class Num_T>
1588 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
1589  Mat<Num_T> &out)
1590 {
1591  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1592  "Mat<>::elem_mult_out(): Wrong sizes");
1593  out.set_size(m1.no_rows, m1.no_cols);
1594  for (int i = 0; i < out.datasize; i++)
1595  out.data[i] = m1.data[i] * m2.data[i];
1596 }
1597 
1598 template<class Num_T>
1599 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
1600  const Mat<Num_T> &m3, Mat<Num_T> &out)
1601 {
1602  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows)
1603  && (m1.no_cols == m2.no_cols) && (m1.no_cols == m3.no_cols),
1604  "Mat<>::elem_mult_out(): Wrong sizes");
1605  out.set_size(m1.no_rows, m1.no_cols);
1606  for (int i = 0; i < out.datasize; i++)
1607  out.data[i] = m1.data[i] * m2.data[i] * m3.data[i];
1608 }
1609 
1610 template<class Num_T>
1611 void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
1612  const Mat<Num_T> &m3, const Mat<Num_T> &m4,
1613  Mat<Num_T> &out)
1614 {
1615  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows)
1616  && (m1.no_rows == m4.no_rows) && (m1.no_cols == m2.no_cols)
1617  && (m1.no_cols == m3.no_cols) && (m1.no_cols == m4.no_cols),
1618  "Mat<>::elem_mult_out(): Wrong sizes");
1619  out.set_size(m1.no_rows, m1.no_cols);
1620  for (int i = 0; i < out.datasize; i++)
1621  out.data[i] = m1.data[i] * m2.data[i] * m3.data[i] * m4.data[i];
1622 }
1623 
1624 template<class Num_T>
1625 #ifndef _MSC_VER
1626 inline
1627 #endif
1629 {
1630  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1631  "Mat<>::elem_mult_inplace(): Wrong sizes");
1632  for (int i = 0; i < m2.datasize; i++)
1633  m2.data[i] *= m1.data[i];
1634 }
1635 
1636 template<class Num_T> inline
1637 Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
1638 {
1639  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1640  "Mat<>::elem_mult_sum(): Wrong sizes");
1641  Num_T acc = 0;
1642 
1643  for (int i = 0; i < m1.datasize; i++)
1644  acc += m1.data[i] * m2.data[i];
1645 
1646  return acc;
1647 }
1648 
1649 template<class Num_T> inline
1651 {
1652  for (int i = 0; i < datasize; i++)
1653  data[i] /= t;
1654  return *this;
1655 }
1656 
1657 template<class Num_T> inline
1659 {
1660  it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols),
1661  "Mat<>::operator/=(): Wrong sizes");
1662  for (int i = 0; i < datasize; i++)
1663  data[i] /= m.data[i];
1664  return *this;
1665 }
1666 
1667 template<class Num_T>
1668 Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t)
1669 {
1670  Mat<Num_T> r(m.no_rows, m.no_cols);
1671  for (int i = 0; i < r.datasize; ++i)
1672  r.data[i] = m.data[i] / t;
1673  return r;
1674 }
1675 
1676 template<class Num_T>
1677 Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m)
1678 {
1679  Mat<Num_T> r(m.no_rows, m.no_cols);
1680  for (int i = 0; i < r.datasize; ++i)
1681  r.data[i] = t / m.data[i];
1682  return r;
1683 }
1684 
1685 template<class Num_T> inline
1687 {
1688  Mat<Num_T> out;
1689  elem_div_out(m1, m2, out);
1690  return out;
1691 }
1692 
1693 template<class Num_T>
1694 void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
1695  Mat<Num_T> &out)
1696 {
1697  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1698  "Mat<>::elem_div_out(): Wrong sizes");
1699 
1700  if ((out.no_rows != m1.no_rows) || (out.no_cols != m1.no_cols))
1701  out.set_size(m1.no_rows, m1.no_cols);
1702 
1703  for (int i = 0; i < out.datasize; i++)
1704  out.data[i] = m1.data[i] / m2.data[i];
1705 }
1706 
1707 template<class Num_T> inline
1708 Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
1709 {
1710  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
1711  "Mat<>::elem_div_sum(): Wrong sizes");
1712  Num_T acc = 0;
1713 
1714  for (int i = 0; i < m1.datasize; i++)
1715  acc += m1.data[i] / m2.data[i];
1716 
1717  return acc;
1718 }
1719 
1720 template<class Num_T>
1722 {
1723  if (no_rows != m.no_rows || no_cols != m.no_cols) return false;
1724  for (int i = 0;i < datasize;i++) {
1725  if (data[i] != m.data[i]) return false;
1726  }
1727  return true;
1728 }
1729 
1730 template<class Num_T>
1732 {
1733  if (no_rows != m.no_rows || no_cols != m.no_cols) return true;
1734  for (int i = 0;i < datasize;i++) {
1735  if (data[i] != m.data[i]) return true;
1736  }
1737  return false;
1738 }
1739 
1740 template <class Num_T>
1741 std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m)
1742 {
1743  int i;
1744 
1745  switch (m.rows()) {
1746  case 0 :
1747  os << "[]";
1748  break;
1749  case 1 :
1750  os << '[' << m.get_row(0) << ']';
1751  break;
1752  default:
1753  os << '[' << m.get_row(0) << std::endl;
1754  for (i = 1; i < m.rows() - 1; i++)
1755  os << ' ' << m.get_row(i) << std::endl;
1756  os << ' ' << m.get_row(m.rows() - 1) << ']';
1757  }
1758 
1759  return os;
1760 }
1761 
1762 template <class Num_T>
1763 std::istream &operator>>(std::istream &is, Mat<Num_T> &m)
1764 {
1765  std::ostringstream buffer;
1766  bool started = false;
1767  bool finished = false;
1768  bool brackets = false;
1769  bool within_double_brackets = false;
1770  char c;
1771 
1772  while (!finished) {
1773  if (is.eof()) {
1774  finished = true;
1775  }
1776  else {
1777  is.get(c);
1778 
1779  if (is.eof() || (c == '\n')) {
1780  if (brackets) {
1781  // Right bracket missing
1782  is.setstate(std::ios_base::failbit);
1783  finished = true;
1784  }
1785  else if (!((c == '\n') && !started)) {
1786  finished = true;
1787  }
1788  }
1789  else if ((c == ' ') || (c == '\t')) {
1790  if (started) {
1791  buffer << ' ';
1792  }
1793  }
1794  else if (c == '[') {
1795  if ((started && !brackets) || within_double_brackets) {
1796  // Unexpected left bracket
1797  is.setstate(std::ios_base::failbit);
1798  finished = true;
1799  }
1800  else if (!started) {
1801  started = true;
1802  brackets = true;
1803  }
1804  else {
1805  within_double_brackets = true;
1806  }
1807  }
1808  else if (c == ']') {
1809  if (!started || !brackets) {
1810  // Unexpected right bracket
1811  is.setstate(std::ios_base::failbit);
1812  finished = true;
1813  }
1814  else if (within_double_brackets) {
1815  within_double_brackets = false;
1816  buffer << ';';
1817  }
1818  else {
1819  finished = true;
1820  }
1821  while (!is.eof() && (((c = static_cast<char>(is.peek())) == ' ')
1822  || (c == '\t'))) {
1823  is.get();
1824  }
1825  if (!is.eof() && (c == '\n')) {
1826  is.get();
1827  }
1828  }
1829  else {
1830  started = true;
1831  buffer << c;
1832  }
1833  }
1834  }
1835 
1836  if (!started) {
1837  m.set_size(0, false);
1838  }
1839  else {
1840  m.set(buffer.str());
1841  }
1842 
1843  return is;
1844 }
1845 
1847 
1848 // ---------------------------------------------------------------------
1849 // Instantiations
1850 // ---------------------------------------------------------------------
1851 
1852 #ifndef _MSC_VER
1853 
1854 // class instantiations
1855 
1856 extern template class Mat<double>;
1857 extern template class Mat<std::complex<double> >;
1858 extern template class Mat<int>;
1859 extern template class Mat<short int>;
1860 extern template class Mat<bin>;
1861 
1862 // addition operators
1863 
1864 extern template mat operator+(const mat &m1, const mat &m2);
1865 extern template cmat operator+(const cmat &m1, const cmat &m2);
1866 extern template imat operator+(const imat &m1, const imat &m2);
1867 extern template smat operator+(const smat &m1, const smat &m2);
1868 extern template bmat operator+(const bmat &m1, const bmat &m2);
1869 
1870 extern template mat operator+(const mat &m, double t);
1871 extern template cmat operator+(const cmat &m, std::complex<double> t);
1872 extern template imat operator+(const imat &m, int t);
1873 extern template smat operator+(const smat &m, short t);
1874 extern template bmat operator+(const bmat &m, bin t);
1875 
1876 extern template mat operator+(double t, const mat &m);
1877 extern template cmat operator+(std::complex<double> t, const cmat &m);
1878 extern template imat operator+(int t, const imat &m);
1879 extern template smat operator+(short t, const smat &m);
1880 extern template bmat operator+(bin t, const bmat &m);
1881 
1882 // subtraction operators
1883 
1884 extern template mat operator-(const mat &m1, const mat &m2);
1885 extern template cmat operator-(const cmat &m1, const cmat &m2);
1886 extern template imat operator-(const imat &m1, const imat &m2);
1887 extern template smat operator-(const smat &m1, const smat &m2);
1888 extern template bmat operator-(const bmat &m1, const bmat &m2);
1889 
1890 extern template mat operator-(const mat &m, double t);
1891 extern template cmat operator-(const cmat &m, std::complex<double> t);
1892 extern template imat operator-(const imat &m, int t);
1893 extern template smat operator-(const smat &m, short t);
1894 extern template bmat operator-(const bmat &m, bin t);
1895 
1896 extern template mat operator-(double t, const mat &m);
1897 extern template cmat operator-(std::complex<double> t, const cmat &m);
1898 extern template imat operator-(int t, const imat &m);
1899 extern template smat operator-(short t, const smat &m);
1900 extern template bmat operator-(bin t, const bmat &m);
1901 
1902 // unary minus
1903 
1904 extern template mat operator-(const mat &m);
1905 extern template cmat operator-(const cmat &m);
1906 extern template imat operator-(const imat &m);
1907 extern template smat operator-(const smat &m);
1908 extern template bmat operator-(const bmat &m);
1909 
1910 // multiplication operators
1911 
1912 extern template imat operator*(const imat &m1, const imat &m2);
1913 extern template smat operator*(const smat &m1, const smat &m2);
1914 extern template bmat operator*(const bmat &m1, const bmat &m2);
1915 
1916 extern template ivec operator*(const imat &m, const ivec &v);
1917 extern template svec operator*(const smat &m, const svec &v);
1918 extern template bvec operator*(const bmat &m, const bvec &v);
1919 
1920 extern template mat operator*(const mat &m, double t);
1921 extern template cmat operator*(const cmat &m, std::complex<double> t);
1922 extern template imat operator*(const imat &m, int t);
1923 extern template smat operator*(const smat &m, short t);
1924 extern template bmat operator*(const bmat &m, bin t);
1925 
1926 extern template mat operator*(double t, const mat &m);
1927 extern template cmat operator*(std::complex<double> t, const cmat &m);
1928 extern template imat operator*(int t, const imat &m);
1929 extern template smat operator*(short t, const smat &m);
1930 extern template bmat operator*(bin t, const bmat &m);
1931 
1932 // element-wise multiplication
1933 
1934 extern template mat elem_mult(const mat &m1, const mat &m2);
1935 extern template cmat elem_mult(const cmat &m1, const cmat &m2);
1936 extern template imat elem_mult(const imat &m1, const imat &m2);
1937 extern template smat elem_mult(const smat &m1, const smat &m2);
1938 extern template bmat elem_mult(const bmat &m1, const bmat &m2);
1939 
1940 extern template void elem_mult_out(const mat &m1, const mat &m2, mat &out);
1941 extern template void elem_mult_out(const cmat &m1, const cmat &m2,
1942  cmat &out);
1943 extern template void elem_mult_out(const imat &m1, const imat &m2,
1944  imat &out);
1945 extern template void elem_mult_out(const smat &m1, const smat &m2,
1946  smat &out);
1947 extern template void elem_mult_out(const bmat &m1, const bmat &m2,
1948  bmat &out);
1949 
1950 extern template void elem_mult_out(const mat &m1, const mat &m2,
1951  const mat &m3, mat &out);
1952 extern template void elem_mult_out(const cmat &m1, const cmat &m2,
1953  const cmat &m3, cmat &out);
1954 extern template void elem_mult_out(const imat &m1, const imat &m2,
1955  const imat &m3, imat &out);
1956 extern template void elem_mult_out(const smat &m1, const smat &m2,
1957  const smat &m3, smat &out);
1958 extern template void elem_mult_out(const bmat &m1, const bmat &m2,
1959  const bmat &m3, bmat &out);
1960 
1961 extern template void elem_mult_out(const mat &m1, const mat &m2,
1962  const mat &m3, const mat &m4, mat &out);
1963 extern template void elem_mult_out(const cmat &m1, const cmat &m2,
1964  const cmat &m3, const cmat &m4,
1965  cmat &out);
1966 extern template void elem_mult_out(const imat &m1, const imat &m2,
1967  const imat &m3, const imat &m4,
1968  imat &out);
1969 extern template void elem_mult_out(const smat &m1, const smat &m2,
1970  const smat &m3, const smat &m4,
1971  smat &out);
1972 extern template void elem_mult_out(const bmat &m1, const bmat &m2,
1973  const bmat &m3, const bmat &m4,
1974  bmat &out);
1975 
1976 extern template void elem_mult_inplace(const mat &m1, mat &m2);
1977 extern template void elem_mult_inplace(const cmat &m1, cmat &m2);
1978 extern template void elem_mult_inplace(const imat &m1, imat &m2);
1979 extern template void elem_mult_inplace(const smat &m1, smat &m2);
1980 extern template void elem_mult_inplace(const bmat &m1, bmat &m2);
1981 
1982 extern template double elem_mult_sum(const mat &m1, const mat &m2);
1983 extern template std::complex<double> elem_mult_sum(const cmat &m1,
1984  const cmat &m2);
1985 extern template int elem_mult_sum(const imat &m1, const imat &m2);
1986 extern template short elem_mult_sum(const smat &m1, const smat &m2);
1987 extern template bin elem_mult_sum(const bmat &m1, const bmat &m2);
1988 
1989 // division operator
1990 
1991 extern template mat operator/(double t, const mat &m);
1992 extern template cmat operator/(std::complex<double> t, const cmat &m);
1993 extern template imat operator/(int t, const imat &m);
1994 extern template smat operator/(short t, const smat &m);
1995 extern template bmat operator/(bin t, const bmat &m);
1996 
1997 extern template mat operator/(const mat &m, double t);
1998 extern template cmat operator/(const cmat &m, std::complex<double> t);
1999 extern template imat operator/(const imat &m, int t);
2000 extern template smat operator/(const smat &m, short t);
2001 extern template bmat operator/(const bmat &m, bin t);
2002 
2003 // element-wise division
2004 
2005 extern template mat elem_div(const mat &m1, const mat &m2);
2006 extern template cmat elem_div(const cmat &m1, const cmat &m2);
2007 extern template imat elem_div(const imat &m1, const imat &m2);
2008 extern template smat elem_div(const smat &m1, const smat &m2);
2009 extern template bmat elem_div(const bmat &m1, const bmat &m2);
2010 
2011 extern template void elem_div_out(const mat &m1, const mat &m2, mat &out);
2012 extern template void elem_div_out(const cmat &m1, const cmat &m2, cmat &out);
2013 extern template void elem_div_out(const imat &m1, const imat &m2, imat &out);
2014 extern template void elem_div_out(const smat &m1, const smat &m2, smat &out);
2015 extern template void elem_div_out(const bmat &m1, const bmat &m2, bmat &out);
2016 
2017 extern template double elem_div_sum(const mat &m1, const mat &m2);
2018 extern template std::complex<double> elem_div_sum(const cmat &m1,
2019  const cmat &m2);
2020 extern template int elem_div_sum(const imat &m1, const imat &m2);
2021 extern template short elem_div_sum(const smat &m1, const smat &m2);
2022 extern template bin elem_div_sum(const bmat &m1, const bmat &m2);
2023 
2024 // concatenation
2025 
2026 extern template mat concat_horizontal(const mat &m1, const mat &m2);
2027 extern template cmat concat_horizontal(const cmat &m1, const cmat &m2);
2028 extern template imat concat_horizontal(const imat &m1, const imat &m2);
2029 extern template smat concat_horizontal(const smat &m1, const smat &m2);
2030 extern template bmat concat_horizontal(const bmat &m1, const bmat &m2);
2031 
2032 extern template mat concat_vertical(const mat &m1, const mat &m2);
2033 extern template cmat concat_vertical(const cmat &m1, const cmat &m2);
2034 extern template imat concat_vertical(const imat &m1, const imat &m2);
2035 extern template smat concat_vertical(const smat &m1, const smat &m2);
2036 extern template bmat concat_vertical(const bmat &m1, const bmat &m2);
2037 
2038 // I/O streams
2039 
2040 extern template std::ostream &operator<<(std::ostream &os, const mat &m);
2041 extern template std::ostream &operator<<(std::ostream &os, const cmat &m);
2042 extern template std::ostream &operator<<(std::ostream &os, const imat &m);
2043 extern template std::ostream &operator<<(std::ostream &os, const smat &m);
2044 extern template std::ostream &operator<<(std::ostream &os, const bmat &m);
2045 
2046 extern template std::istream &operator>>(std::istream &is, mat &m);
2047 extern template std::istream &operator>>(std::istream &is, cmat &m);
2048 extern template std::istream &operator>>(std::istream &is, imat &m);
2049 extern template std::istream &operator>>(std::istream &is, smat &m);
2050 extern template std::istream &operator>>(std::istream &is, bmat &m);
2051 
2052 #endif // _MSC_VER
2053 
2055 
2056 } // namespace itpp
2057 
2058 #endif // #ifndef MAT_H
SourceForge Logo

Generated on Fri Mar 21 2014 17:14:12 for IT++ by Doxygen 1.8.1.2