mat.h

Go to the documentation of this file.

#ifndef MAT_H
#define MAT_H

#include <itpp/base/itassert.h>
#include <itpp/base/math/misc.h>
#include <itpp/base/factory.h>
#include <itpp/itexports.h>

namespace itpp
{

// Declaration of Vec
template<class Num_T> class Vec;
// Declaration of Mat
template<class Num_T> class Mat;
// Declaration of bin
class bin;

template<class Num_T>
Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

template<class Num_T>
Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t);
template<class Num_T>
Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m);

template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t);
template<class Num_T>
Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m);
template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m);

template<class Num_T>
Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v);
template<class Num_T>
Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t);
template<class Num_T>
Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m);

template<class Num_T>
Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   Mat<Num_T> &out);
template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   const Mat<Num_T> &m3, Mat<Num_T> &out);
template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   const Mat<Num_T> &m3, const Mat<Num_T> &m4,
                   Mat<Num_T> &out);
template<class Num_T>
void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2);
template<class Num_T>
Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

template<class Num_T>
Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t);
template<class Num_T>
Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m);

template<class Num_T>
Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
template<class Num_T>
void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                  Mat<Num_T> &out);
template<class Num_T>
Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

// -------------------------------------------------------------------------------------
// Declaration of Mat
// -------------------------------------------------------------------------------------

template<class Num_T>
class Mat
{
public:
  typedef Num_T value_type;

  explicit Mat(const Factory &f = DEFAULT_FACTORY);
  Mat(int rows, int cols, const Factory &f = DEFAULT_FACTORY);
  Mat(const Mat<Num_T> &m);
  Mat(const Mat<Num_T> &m, const Factory &f);
  Mat(const Vec<Num_T> &v, const Factory &f = DEFAULT_FACTORY);
  Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY);
  Mat(const char *str, const Factory &f = DEFAULT_FACTORY);
  Mat(const Num_T *c_array, int rows, int cols, bool row_major = true,
      const Factory &f = DEFAULT_FACTORY);

  ~Mat();

  int cols() const { return no_cols; }
  int rows() const { return no_rows; }
  int size() const { return datasize; }
  void set_size(int rows, int cols, bool copy = false);
  void zeros();
  void clear() { zeros(); }
  void ones();
  void set(const std::string &str);
  void set(const char *str);

  const Num_T &operator()(int r, int c) const;
  Num_T &operator()(int r, int c);
  const Num_T &operator()(int i) const;
  Num_T &operator()(int i);
  const Num_T &get(int r, int c) const;
  const Num_T &get(int i) const;
  void set(int r, int c, Num_T t);

  Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const;
  Mat<Num_T> get(int r1, int r2, int c1, int c2) const;

  Vec<Num_T> get_row(int r) const;
  Mat<Num_T> get_rows(int r1, int r2) const;
  Mat<Num_T> get_rows(const Vec<int> &indexlist) const;
  Vec<Num_T> get_col(int c) const;
  Mat<Num_T> get_cols(int c1, int c2) const;
  Mat<Num_T> get_cols(const Vec<int> &indexlist) const;
  void set_row(int r, const Vec<Num_T> &v);
  void set_col(int c, const Vec<Num_T> &v);
  void set_rows(int r, const Mat<Num_T> &m);
  void set_cols(int c, const Mat<Num_T> &m);
  void copy_row(int to, int from);
  void copy_col(int to, int from);
  void swap_rows(int r1, int r2);
  void swap_cols(int c1, int c2);

  void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m);
  void set_submatrix(int r, int c, const Mat<Num_T> &m);
  void set_submatrix(int r1, int r2, int c1, int c2, Num_T t);

  void del_row(int r);
  void del_rows(int r1, int r2);
  void del_col(int c);
  void del_cols(int c1, int c2);
  void ins_row(int r, const Vec<Num_T> &v);
  void ins_col(int c, const Vec<Num_T> &v);
  void append_row(const Vec<Num_T> &v);
  void append_col(const Vec<Num_T> &v);

  Mat<Num_T> transpose() const;
  Mat<Num_T> T() const { return this->transpose(); }
  Mat<Num_T> hermitian_transpose() const;
  Mat<Num_T> H() const { return this->hermitian_transpose(); }

  friend Mat<Num_T> concat_horizontal<>(const Mat<Num_T> &m1,
                                        const Mat<Num_T> &m2);
  friend Mat<Num_T> concat_vertical<>(const Mat<Num_T> &m1,
                                      const Mat<Num_T> &m2);

  Mat<Num_T>& operator=(Num_T t);
  Mat<Num_T>& operator=(const Mat<Num_T> &m);
  Mat<Num_T>& operator=(const Vec<Num_T> &v);
  Mat<Num_T>& operator=(const std::string &str);
  Mat<Num_T>& operator=(const char *str);

  Mat<Num_T>& operator+=(const Mat<Num_T> &m);
  Mat<Num_T>& operator+=(Num_T t);
  friend Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  friend Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t);
  friend Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m);

  Mat<Num_T>& operator-=(const Mat<Num_T> &m);
  Mat<Num_T>& operator-=(Num_T t);
  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t);
  friend Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m);
  friend Mat<Num_T> operator-<>(const Mat<Num_T> &m);

  Mat<Num_T>& operator*=(const Mat<Num_T> &m);
  Mat<Num_T>& operator*=(Num_T t);

  friend Mat<Num_T> elem_mult<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                              Mat<Num_T> &out);
  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                              const Mat<Num_T> &m3, Mat<Num_T> &out);
  friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                              const Mat<Num_T> &m3, const Mat<Num_T> &m4,
                              Mat<Num_T> &out);
  friend void elem_mult_inplace<>(const Mat<Num_T> &m1, Mat<Num_T> &m2);
  friend Num_T elem_mult_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

  Mat<Num_T>& operator/=(Num_T t);
  Mat<Num_T>& operator/=(const Mat<Num_T> &m);

  friend Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t);
  friend Mat<Num_T> operator/<>(Num_T t, const Mat<Num_T> &m);

  friend Mat<Num_T> elem_div<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  friend void elem_div_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                             Mat<Num_T> &out);
  friend Num_T elem_div_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

  bool operator==(const Mat<Num_T> &m) const;
  bool operator!=(const Mat<Num_T> &m) const;

  Num_T &_elem(int r, int c) { return data[r+c*no_rows]; }
  const Num_T &_elem(int r, int c) const { return data[r+c*no_rows]; }
  Num_T &_elem(int i) { return data[i]; }
  const Num_T &_elem(int i) const { return data[i]; }

  Num_T *_data() { return data; }
  const Num_T *_data() const { return data; }
  int _datasize() const { return datasize; }

protected:
  void alloc(int rows, int cols);
  void free();

  int datasize, no_rows, no_cols;

  Num_T *data;
  const Factory &factory;

private:
  bool in_range(int r, int c) const {
    return ((r >= 0) && (r < no_rows) && (c >= 0) && (c < no_cols));
  }
  bool row_in_range(int r) const { return ((r >= 0) && (r < no_rows)); }
  bool col_in_range(int c) const { return ((c >= 0) && (c < no_cols)); }
  bool in_range(int i) const { return ((i >= 0) && (i < datasize)); }
};

// -------------------------------------------------------------------------------------
// Type definitions of mat, cmat, imat, smat, and bmat
// -------------------------------------------------------------------------------------

typedef Mat<double> mat;

typedef Mat<std::complex<double> > cmat;

typedef Mat<int> imat;

typedef Mat<short int> smat;

typedef Mat<bin> bmat;

} //namespace itpp


#include <itpp/base/vec.h>

namespace itpp
{

// ----------------------------------------------------------------------
// Declaration of input and output streams for Mat
// ----------------------------------------------------------------------

template <class Num_T>
std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m);

template <class Num_T>
std::istream &operator>>(std::istream &is, Mat<Num_T> &m);

// ----------------------------------------------------------------------
// Implementation of templated Mat members and friends
// ----------------------------------------------------------------------

template<class Num_T> inline
void Mat<Num_T>::alloc(int rows, int cols)
{
  if ((rows > 0) && (cols > 0)) {
    datasize = rows * cols;
    no_rows = rows;
    no_cols = cols;
    create_elements(data, datasize, factory);
  }
  else {
    data = 0;
    datasize = 0;
    no_rows = 0;
    no_cols = 0;
  }
}

template<class Num_T> inline
void Mat<Num_T>::free()
{
  destroy_elements(data, datasize);
  datasize = 0;
  no_rows = 0;
  no_cols = 0;
}


template<class Num_T> inline
Mat<Num_T>::Mat(const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f) {}

template<class Num_T> inline
Mat<Num_T>::Mat(int rows, int cols, const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  it_assert_debug((rows >= 0) && (cols >= 0), "Mat<>::Mat(): Wrong size");
  alloc(rows, cols);
}

template<class Num_T> inline
Mat<Num_T>::Mat(const Mat<Num_T> &m) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(m.factory)
{
  alloc(m.no_rows, m.no_cols);
  copy_vector(m.datasize, m.data, data);
}

template<class Num_T> inline
Mat<Num_T>::Mat(const Mat<Num_T> &m, const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  alloc(m.no_rows, m.no_cols);
  copy_vector(m.datasize, m.data, data);
}

template<class Num_T> inline
Mat<Num_T>::Mat(const Vec<Num_T> &v, const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  int size = v.size();
  alloc(size, 1);
  copy_vector(size, v._data(), data);
}

template<class Num_T> inline
Mat<Num_T>::Mat(const std::string &str, const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  set(str);
}

template<class Num_T> inline
Mat<Num_T>::Mat(const char *str, const Factory &f) :
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  set(std::string(str));
}

template<class Num_T>
Mat<Num_T>::Mat(const Num_T *c_array, int rows, int cols, bool row_major,
                const Factory &f):
    datasize(0), no_rows(0), no_cols(0), data(0), factory(f)
{
  alloc(rows, cols);
  if (!row_major)
    copy_vector(datasize, c_array, data);
  else
    for (int i = 0; i < rows; i++)
      for (int j = 0; j < cols; j++)
        data[i+j*no_rows] = c_array[i*no_cols+j];
}

template<class Num_T> inline
Mat<Num_T>::~Mat()
{
  free();
}


template<class Num_T>
void Mat<Num_T>::set_size(int rows, int cols, bool copy)
{
  it_assert_debug((rows >= 0) && (cols >= 0),
                  "Mat<>::set_size(): Wrong size");
  // check if we have to resize the current matrix
  if ((no_rows == rows) && (no_cols == cols))
    return;
  // check if one of dimensions is zero
  if ((rows == 0) || (cols == 0)) {
    free();
    return;
  }
  // conditionally copy previous matrix content
  if (copy) {
    // create a temporary pointer to the allocated data
    Num_T* tmp = data;
    // store the current number of elements and number of rows
    int old_datasize = datasize;
    int old_rows = no_rows;
    // check the boundaries of the copied data
    int min_r = (no_rows < rows) ? no_rows : rows;
    int min_c = (no_cols < cols) ? no_cols : cols;
    // allocate new memory
    alloc(rows, cols);
    // copy the previous data into the allocated memory
    for (int i = 0; i < min_c; ++i) {
      copy_vector(min_r, &tmp[i*old_rows], &data[i*no_rows]);
    }
    // fill-in the rest of matrix with zeros
    for (int i = min_r; i < rows; ++i)
      for (int j = 0; j < cols; ++j)
        data[i+j*rows] = Num_T(0);
    for (int j = min_c; j < cols; ++j)
      for (int i = 0; i < min_r; ++i)
        data[i+j*rows] = Num_T(0);
    // delete old elements
    destroy_elements(tmp, old_datasize);
  }
  // if possible, reuse the allocated memory
  else if (datasize == rows * cols) {
    no_rows = rows;
    no_cols = cols;
  }
  // finally release old memory and allocate a new one
  else {
    free();
    alloc(rows, cols);
  }
}

template<class Num_T> inline
void Mat<Num_T>::zeros()
{
  for (int i = 0; i < datasize; i++)
    data[i] = Num_T(0);
}

template<class Num_T> inline
void Mat<Num_T>::ones()
{
  for (int i = 0; i < datasize; i++)
    data[i] = Num_T(1);
}

template<class Num_T> inline
const Num_T& Mat<Num_T>::operator()(int r, int c) const
{
  it_assert_debug(in_range(r, c),
                  "Mat<>::operator(): Indexing out of range");
  return data[r+c*no_rows];
}

template<class Num_T> inline
Num_T& Mat<Num_T>::operator()(int r, int c)
{
  it_assert_debug(in_range(r, c),
                  "Mat<>::operator(): Indexing out of range");
  return data[r+c*no_rows];
}

template<class Num_T> inline
Num_T& Mat<Num_T>::operator()(int i)
{
  it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range");
  return data[i];
}

template<class Num_T> inline
const Num_T& Mat<Num_T>::operator()(int i) const
{
  it_assert_debug(in_range(i), "Mat<>::operator(): Index out of range");
  return data[i];
}

template<class Num_T> inline
const Num_T& Mat<Num_T>::get(int r, int c) const
{
  return (*this)(r, c);
}

template<class Num_T> inline
const Num_T& Mat<Num_T>::get(int i) const
{
  return (*this)(i);
}

template<class Num_T> inline
void Mat<Num_T>::set(int r, int c, Num_T t)
{
  it_assert_debug(in_range(r, c), "Mat<>::set(): Indexing out of range");
  data[r+c*no_rows] = t;
}


template<class Num_T>
void Mat<Num_T>::set(const std::string &str)
{
  // actual row counter
  int rows = 0;
  // number of rows to allocate next time (8, 16, 32, 64, etc.)
  int maxrows = 8;

  // clean the current matrix content
  free();

  // variable to store the start of a current vector
  std::string::size_type beg = 0;
  std::string::size_type end = 0;
  while (end != std::string::npos) {
    // find next occurrence of a semicolon in string str
    end = str.find(';', beg);
    // parse first row into a vector v
    Vec<Num_T> v(str.substr(beg, end - beg));
    int v_size = v.size();

    // this check is necessary to parse the following two strings as the
    // same matrix: "1 0 1; ; 1 1; " and "1 0 1; 0 0 0; 1 1 0"
    if ((end != std::string::npos) || (v_size > 0)) {
      // matrix empty -> insert v as a first row and allocate maxrows
      if (rows == 0) {
        set_size(maxrows, v_size, true);
        set_row(rows++, v);
      }
      else {
        // check if we need to resize the matrix
        if ((rows == maxrows) || (v_size != no_cols)) {
          // we need to add new rows
          if (rows == maxrows) {
            maxrows *= 2;
          }
          // check if we need to add new columns
          if (v_size > no_cols) {
            set_size(maxrows, v_size, true);
          }
          else {
            set_size(maxrows, no_cols, true);
            // set the size of the parsed vector to the number of columns
            v.set_size(no_cols, true);
          }
        }
        // set the parsed vector as the next row
        set_row(rows++, v);
      }
    }
    // update the starting position of the next vector in the parsed
    // string
    beg = end + 1;
  } // if ((end != std::string::npos) || (v.size > 0))

  set_size(rows, no_cols, true);
}

template<class Num_T> inline
void Mat<Num_T>::set(const char *str)
{
  set(std::string(str));
}

template<class Num_T> inline
Mat<Num_T> Mat<Num_T>::operator()(int r1, int r2, int c1, int c2) const
{
  if (r1 == -1) r1 = no_rows - 1;
  if (r2 == -1) r2 = no_rows - 1;
  if (c1 == -1) c1 = no_cols - 1;
  if (c2 == -1) c2 = no_cols - 1;

  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) &&
                  (c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
                  "Mat<>::operator()(r1, r2, c1, c2): Wrong indexing");

  Mat<Num_T> s(r2 - r1 + 1, c2 - c1 + 1);

  for (int i = 0;i < s.no_cols;i++)
    copy_vector(s.no_rows, data + r1 + (c1 + i)*no_rows, s.data + i*s.no_rows);

  return s;
}

template<class Num_T> inline
Mat<Num_T> Mat<Num_T>::get(int r1, int r2, int c1, int c2) const
{
  return (*this)(r1, r2, c1, c2);
}

template<class Num_T> inline
Vec<Num_T> Mat<Num_T>::get_row(int r) const
{
  it_assert_debug(row_in_range(r), "Mat<>::get_row(): Index out of range");
  Vec<Num_T> a(no_cols);

  copy_vector(no_cols, data + r, no_rows, a._data(), 1);
  return a;
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::get_rows(int r1, int r2) const
{
  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows),
                  "Mat<>::get_rows(): Wrong indexing");
  Mat<Num_T> m(r2 - r1 + 1, no_cols);

  for (int i = 0; i < m.rows(); i++)
    copy_vector(no_cols, data + i + r1, no_rows, m.data + i, m.no_rows);

  return m;
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::get_rows(const Vec<int> &indexlist) const
{
  Mat<Num_T> m(indexlist.size(), no_cols);

  for (int i = 0;i < indexlist.size();i++) {
    it_assert_debug(row_in_range(indexlist(i)),
                    "Mat<>::get_rows(indexlist): Indexing out of range");
    copy_vector(no_cols, data + indexlist(i), no_rows, m.data + i, m.no_rows);
  }

  return m;
}

template<class Num_T> inline
Vec<Num_T> Mat<Num_T>::get_col(int c) const
{
  it_assert_debug(col_in_range(c), "Mat<>::get_col(): Index out of range");
  Vec<Num_T> a(no_rows);

  copy_vector(no_rows, data + c*no_rows, a._data());

  return a;
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::get_cols(int c1, int c2) const
{
  it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
                  "Mat<>::get_cols(): Wrong indexing");
  Mat<Num_T> m(no_rows, c2 - c1 + 1);

  for (int i = 0; i < m.cols(); i++)
    copy_vector(no_rows, data + (i + c1)*no_rows, m.data + i*m.no_rows);

  return m;
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::get_cols(const Vec<int> &indexlist) const
{
  Mat<Num_T> m(no_rows, indexlist.size());

  for (int i = 0; i < indexlist.size(); i++) {
    it_assert_debug(col_in_range(indexlist(i)),
                    "Mat<>::get_cols(indexlist): Indexing out of range");
    copy_vector(no_rows, data + indexlist(i)*no_rows, m.data + i*m.no_rows);
  }

  return m;
}

template<class Num_T> inline
void Mat<Num_T>::set_row(int r, const Vec<Num_T> &v)
{
  it_assert_debug(row_in_range(r), "Mat<>::set_row(): Index out of range");
  it_assert_debug(v.size() == no_cols,
                  "Mat<>::set_row(): Wrong size of input vector");
  copy_vector(v.size(), v._data(), 1, data + r, no_rows);
}

template<class Num_T> inline
void Mat<Num_T>::set_col(int c, const Vec<Num_T> &v)
{
  it_assert_debug(col_in_range(c), "Mat<>::set_col(): Index out of range");
  it_assert_debug(v.size() == no_rows,
                  "Mat<>::set_col(): Wrong size of input vector");
  copy_vector(v.size(), v._data(), data + c*no_rows);
}


template<class Num_T>
void Mat<Num_T>::set_rows(int r, const Mat<Num_T> &m)
{
  it_assert_debug(row_in_range(r), "Mat<>::set_rows(): Index out of range");
  it_assert_debug(no_cols == m.cols(),
                  "Mat<>::set_rows(): Column sizes do not match");
  it_assert_debug(m.rows() + r <= no_rows,
                  "Mat<>::set_rows(): Not enough rows");

  for (int i = 0; i < m.rows(); ++i) {
    copy_vector(no_cols, m.data + i, m.no_rows, data + i + r, no_rows);
  }
}

template<class Num_T>
void Mat<Num_T>::set_cols(int c, const Mat<Num_T> &m)
{
  it_assert_debug(col_in_range(c), "Mat<>::set_cols(): Index out of range");
  it_assert_debug(no_rows == m.rows(),
                  "Mat<>::set_cols(): Row sizes do not match");
  it_assert_debug(m.cols() + c <= no_cols,
                  "Mat<>::set_cols(): Not enough colums");

  for (int i = 0; i < m.cols(); ++i) {
    copy_vector(no_rows, m.data + i*no_rows, data + (i + c)*no_rows);
  }
}


template<class Num_T> inline
void Mat<Num_T>::copy_row(int to, int from)
{
  it_assert_debug(row_in_range(to) && row_in_range(from),
                  "Mat<>::copy_row(): Indexing out of range");
  if (from == to)
    return;

  copy_vector(no_cols, data + from, no_rows, data + to, no_rows);
}

template<class Num_T> inline
void Mat<Num_T>::copy_col(int to, int from)
{
  it_assert_debug(col_in_range(to) && col_in_range(from),
                  "Mat<>::copy_col(): Indexing out of range");
  if (from == to)
    return;

  copy_vector(no_rows, data + from*no_rows, data + to*no_rows);
}

template<class Num_T> inline
void Mat<Num_T>::swap_rows(int r1, int r2)
{
  it_assert_debug(row_in_range(r1) && row_in_range(r2),
                  "Mat<>::swap_rows(): Indexing out of range");
  if (r1 == r2)
    return;

  swap_vector(no_cols, data + r1, no_rows, data + r2, no_rows);
}

template<class Num_T> inline
void Mat<Num_T>::swap_cols(int c1, int c2)
{
  it_assert_debug(col_in_range(c1) && col_in_range(c2),
                  "Mat<>::swap_cols(): Indexing out of range");
  if (c1 == c2)
    return;

  swap_vector(no_rows, data + c1*no_rows, data + c2*no_rows);
}

template<class Num_T>
void Mat<Num_T>::set_submatrix(int r1, int, int c1, int, const Mat<Num_T> &m)
{
  it_warning("Mat<>::set_submatrix(r1, r2, r3, r4, m): This function is "
             "deprecated and might be removed from future IT++ releases. "
             "Please use Mat<>::set_submatrix(r, c, m) function instead.");
  set_submatrix(r1, c1, m);
}

template<class Num_T> inline
void Mat<Num_T>::set_submatrix(int r, int c, const Mat<Num_T> &m)
{
  it_assert_debug((r >= 0) && (r + m.no_rows <= no_rows) &&
                  (c >= 0) && (c + m.no_cols <= no_cols),
                  "Mat<>::set_submatrix(): Indexing out of range "
                  "or wrong input matrix");
  for (int i = 0; i < m.no_cols; i++)
    copy_vector(m.no_rows, m.data + i*m.no_rows, data + (c + i)*no_rows + r);
}



template<class Num_T> inline
void Mat<Num_T>::set_submatrix(int r1, int r2, int c1, int c2, Num_T t)
{
  if (r1 == -1) r1 = no_rows - 1;
  if (r2 == -1) r2 = no_rows - 1;
  if (c1 == -1) c1 = no_cols - 1;
  if (c2 == -1) c2 = no_cols - 1;
  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows) &&
                  (c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
                  "Mat<>::set_submatrix(): Wrong indexing");
  for (int i = c1; i <= c2; i++) {
    int pos = i * no_rows + r1;
    for (int j = r1; j <= r2; j++)
      data[pos++] = t;
  }
}

template<class Num_T>
void Mat<Num_T>::del_row(int r)
{
  it_assert_debug(row_in_range(r), "Mat<>::del_row(): Index out of range");
  Mat<Num_T> Temp(*this);
  set_size(no_rows - 1, no_cols, false);
  for (int i = 0 ; i < r ; i++) {
    copy_vector(no_cols, &Temp.data[i], no_rows + 1, &data[i], no_rows);
  }
  for (int i = r ; i < no_rows ; i++) {
    copy_vector(no_cols, &Temp.data[i+1], no_rows + 1, &data[i], no_rows);
  }

}

template<class Num_T>
void Mat<Num_T>::del_rows(int r1, int r2)
{
  it_assert_debug((r1 >= 0) && (r1 <= r2) && (r2 < no_rows),
                  "Mat<>::del_rows(): Indexing out of range");
  Mat<Num_T> Temp(*this);
  int no_del_rows = r2 - r1 + 1;
  set_size(no_rows - no_del_rows, no_cols, false);
  for (int i = 0; i < r1 ; ++i) {
    copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows);
  }
  for (int i = r2 + 1; i < Temp.no_rows; ++i) {
    copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i-no_del_rows],
                no_rows);
  }
}

template<class Num_T>
void Mat<Num_T>::del_col(int c)
{
  it_assert_debug(col_in_range(c), "Mat<>::del_col(): Index out of range");
  Mat<Num_T> Temp(*this);

  set_size(no_rows, no_cols - 1, false);
  copy_vector(c*no_rows, Temp.data, data);
  copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]);
}

template<class Num_T>
void Mat<Num_T>::del_cols(int c1, int c2)
{
  it_assert_debug((c1 >= 0) && (c1 <= c2) && (c2 < no_cols),
                  "Mat<>::del_cols(): Indexing out of range");
  Mat<Num_T> Temp(*this);
  int n_deleted_cols = c2 - c1 + 1;
  set_size(no_rows, no_cols - n_deleted_cols, false);
  copy_vector(c1*no_rows, Temp.data, data);
  copy_vector((no_cols - c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]);
}

template<class Num_T>
void Mat<Num_T>::ins_row(int r, const Vec<Num_T> &v)
{
  it_assert_debug((r >= 0) && (r <= no_rows),
                  "Mat<>::ins_row(): Index out of range");
  it_assert_debug((v.size() == no_cols) || (no_rows == 0),
                  "Mat<>::ins_row(): Wrong size of the input vector");

  if (no_cols == 0) {
    no_cols = v.size();
  }

  Mat<Num_T> Temp(*this);
  set_size(no_rows + 1, no_cols, false);

  for (int i = 0 ; i < r ; i++) {
    copy_vector(no_cols, &Temp.data[i], no_rows - 1, &data[i], no_rows);
  }
  copy_vector(no_cols, v._data(), 1, &data[r], no_rows);
  for (int i = r + 1 ; i < no_rows ; i++) {
    copy_vector(no_cols, &Temp.data[i-1], no_rows - 1, &data[i], no_rows);
  }
}

template<class Num_T>
void Mat<Num_T>::ins_col(int c, const Vec<Num_T> &v)
{
  it_assert_debug((c >= 0) && (c <= no_cols),
                  "Mat<>::ins_col(): Index out of range");
  it_assert_debug((v.size() == no_rows) || (no_cols == 0),
                  "Mat<>::ins_col(): Wrong size of the input vector");

  if (no_rows == 0) {
    no_rows = v.size();
  }

  Mat<Num_T> Temp(*this);
  set_size(no_rows, no_cols + 1, false);

  copy_vector(c*no_rows, Temp.data, data);
  copy_vector(no_rows, v._data(), &data[c*no_rows]);
  copy_vector((no_cols - c - 1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]);
}

template<class Num_T> inline
void Mat<Num_T>::append_row(const Vec<Num_T> &v)
{
  ins_row(no_rows, v);
}

template<class Num_T> inline
void Mat<Num_T>::append_col(const Vec<Num_T> &v)
{
  ins_col(no_cols, v);
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::transpose() const
{
  Mat<Num_T> temp(no_cols, no_rows);
  for (int i = 0; i < no_rows; ++i) {
    copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1);
  }
  return temp;
}

template<class Num_T>
Mat<Num_T> Mat<Num_T>::hermitian_transpose() const
{
  Mat<Num_T> temp(no_cols, no_rows);
  for (int i = 0; i < no_rows; ++i) {
    copy_vector(no_cols, &data[i], no_rows, &temp.data[i * no_cols], 1);
  }
  return temp;
}

template<>
ITPP_EXPORT cmat cmat::hermitian_transpose() const;

template<class Num_T>
Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  // if one of the input matrix is empty just copy the other one as a result
  if (m1.no_cols == 0)
    return m2;
  if (m2.no_cols == 0)
    return m1;
  it_assert_debug(m1.no_rows == m2.no_rows,
                  "Mat<>::concat_horizontal(): Wrong sizes");
  int no_rows = m1.no_rows;
  Mat<Num_T> temp(no_rows, m1.no_cols + m2.no_cols);
  for (int i = 0; i < m1.no_cols; ++i) {
    copy_vector(no_rows, &m1.data[i * no_rows], &temp.data[i * no_rows]);
  }
  for (int i = 0; i < m2.no_cols; ++i) {
    copy_vector(no_rows, &m2.data[i * no_rows], &temp.data[(m1.no_cols + i)
                * no_rows]);
  }
  return temp;
}

template<class Num_T>
Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  // if one of the input matrix is empty just copy the other one as a result
  if (m1.no_rows == 0)
    return m2;
  if (m2.no_rows == 0)
    return m1;
  it_assert_debug(m1.no_cols == m2.no_cols,
                  "Mat<>::concat_vertical(): Wrong sizes");
  int no_cols = m1.no_cols;
  Mat<Num_T> temp(m1.no_rows + m2.no_rows, no_cols);
  for (int i = 0; i < no_cols; ++i) {
    copy_vector(m1.no_rows, &m1.data[i * m1.no_rows],
                &temp.data[i * temp.no_rows]);
    copy_vector(m2.no_rows, &m2.data[i * m2.no_rows],
                &temp.data[i * temp.no_rows + m1.no_rows]);
  }
  return temp;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator=(Num_T t)
{
  for (int i = 0; i < datasize; i++)
    data[i] = t;
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator=(const Mat<Num_T> &m)
{
  if (this != &m) {
    set_size(m.no_rows, m.no_cols, false);
    if (m.datasize != 0)
      copy_vector(m.datasize, m.data, data);
  }
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator=(const Vec<Num_T> &v)
{
  it_assert_debug(((no_rows == 1) && (no_cols == v.size()))
                  || ((no_cols == 1) && (no_rows == v.size())),
                  "Mat<>::operator=(): Wrong size of the input vector");
  set_size(v.size(), 1, false);
  copy_vector(v.size(), v._data(), data);
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator=(const std::string &str)
{
  set(str);
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator=(const char *str)
{
  set(std::string(str));
  return *this;
}

//-------------------- Templated friend functions --------------------------

template<class Num_T>
Mat<Num_T>& Mat<Num_T>::operator+=(const Mat<Num_T> &m)
{
  if (datasize == 0)
    operator=(m);
  else {
    int i, j, m_pos = 0, pos = 0;
    it_assert_debug(m.no_rows == no_rows && m.no_cols == no_cols, "Mat<Num_T>::operator+=: wrong sizes");
    for (i = 0; i < no_cols; i++) {
      for (j = 0; j < no_rows; j++)
        data[pos+j] += m.data[m_pos+j];
      pos += no_rows;
      m_pos += m.no_rows;
    }
  }
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator+=(Num_T t)
{
  for (int i = 0; i < datasize; i++)
    data[i] += t;
  return *this;
}

template<class Num_T>
Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  Mat<Num_T> r(m1.no_rows, m1.no_cols);
  int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0;

  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::operator+(): Wrong sizes");

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++)
      r.data[r_pos+j] = m1.data[m1_pos+j] + m2.data[m2_pos+j];
    // next column
    m1_pos += m1.no_rows;
    m2_pos += m2.no_rows;
    r_pos += r.no_rows;
  }

  return r;
}


template<class Num_T>
Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);

  for (int i = 0; i < r.datasize; i++)
    r.data[i] = m.data[i] + t;

  return r;
}

template<class Num_T>
Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);

  for (int i = 0; i < r.datasize; i++)
    r.data[i] = t + m.data[i];

  return r;
}

template<class Num_T>
Mat<Num_T>& Mat<Num_T>::operator-=(const Mat<Num_T> &m)
{
  int i, j, m_pos = 0, pos = 0;

  if (datasize == 0) {
    set_size(m.no_rows, m.no_cols, false);
    for (i = 0; i < no_cols; i++) {
      for (j = 0; j < no_rows; j++)
        data[pos+j] = -m.data[m_pos+j];
      // next column
      m_pos += m.no_rows;
      pos += no_rows;
    }
  }
  else {
    it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols),
                    "Mat<>::operator-=(): Wrong sizes");
    for (i = 0; i < no_cols; i++) {
      for (j = 0; j < no_rows; j++)
        data[pos+j] -= m.data[m_pos+j];
      // next column
      m_pos += m.no_rows;
      pos += no_rows;
    }
  }
  return *this;
}

template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  Mat<Num_T> r(m1.no_rows, m1.no_cols);
  int i, j, m1_pos = 0, m2_pos = 0, r_pos = 0;
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::operator-(): Wrong sizes");

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++)
      r.data[r_pos+j] = m1.data[m1_pos+j] - m2.data[m2_pos+j];
    // next column
    m1_pos += m1.no_rows;
    m2_pos += m2.no_rows;
    r_pos += r.no_rows;
  }

  return r;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator-=(Num_T t)
{
  for (int i = 0; i < datasize; i++)
    data[i] -= t;
  return *this;
}

template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);
  int i, j, m_pos = 0, r_pos = 0;

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++)
      r.data[r_pos+j] = m.data[m_pos+j] - t;
    // next column
    m_pos += m.no_rows;
    r_pos += r.no_rows;
  }

  return r;
}

template<class Num_T>
Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);
  int i, j, m_pos = 0, r_pos = 0;

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++)
      r.data[r_pos+j] = t - m.data[m_pos+j];
    // next column
    m_pos += m.no_rows;
    r_pos += r.no_rows;
  }

  return r;
}

template<class Num_T>
Mat<Num_T> operator-(const Mat<Num_T> &m)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);
  int i, j, m_pos = 0, r_pos = 0;

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++)
      r.data[r_pos+j] = -m.data[m_pos+j];
    // next column
    m_pos += m.no_rows;
    r_pos += r.no_rows;
  }

  return r;
}

template<class Num_T>
Mat<Num_T>& Mat<Num_T>::operator*=(const Mat<Num_T> &m)
{
  it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
  Mat<Num_T> r(no_rows, m.no_cols);

  Num_T tmp;

  int i, j, k, r_pos = 0, pos = 0, m_pos = 0;

  for (i = 0; i < r.no_cols; i++) {
    for (j = 0; j < r.no_rows; j++) {
      tmp = Num_T(0);
      pos = 0;
      for (k = 0; k < no_cols; k++) {
        tmp += data[pos+j] * m.data[m_pos+k];
        pos += no_rows;
      }
      r.data[r_pos+j] = tmp;
    }
    r_pos += r.no_rows;
    m_pos += m.no_rows;
  }
  operator=(r); // time consuming
  return *this;
}

template<> ITPP_EXPORT mat& mat::operator*=(const mat &m);
template<> ITPP_EXPORT cmat& cmat::operator*=(const cmat &m);

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator*=(Num_T t)
{
  scal_vector(datasize, t, data);
  return *this;
}

template<class Num_T>
Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  it_assert_debug(m1.cols() == m2.rows(),
                  "Mat<>::operator*(): Wrong sizes");
  Mat<Num_T> r(m1.rows(), m2.cols());

  Num_T tmp;
  int i, j, k;
  Num_T *tr = r._data();
  const Num_T *t1; const Num_T *t2 = m2._data();

  for (i = 0; i < r.cols(); i++) {
    for (j = 0; j < r.rows(); j++) {
      tmp = Num_T(0);
      t1 = m1._data() + j;
      for (k = m1.cols(); k > 0; k--) {
        tmp += *(t1) * *(t2++);
            t1 += m1.rows();
      }
      *(tr++) = tmp;
          t2 -= m2.rows();
    }
    t2 += m2.rows();
  }

  return r;
}

template<> ITPP_EXPORT mat operator*(const mat &m1, const mat &m2);
template<> ITPP_EXPORT cmat operator*(const cmat &m1, const cmat &m2);

template<class Num_T>
Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v)
{
  it_assert_debug(m.cols() == v.size(),
                  "Mat<>::operator*(): Wrong sizes");
  Vec<Num_T> r(m.rows());
  int i, k, m_pos;

  for (i = 0; i < m.rows(); i++) {
    r(i) = Num_T(0);
    m_pos = 0;
    for (k = 0; k < m.cols(); k++) {
      r(i) += m._data()[m_pos+i] * v(k);
      m_pos += m.rows();
    }
  }

  return r;
}

template<> ITPP_EXPORT vec operator*(const mat &m, const vec &v);
template<> ITPP_EXPORT cvec operator*(const cmat &m, const cvec &v);

template<class Num_T>
Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t)
{
  Mat<Num_T> r(m.rows(), m.cols());

  const Num_T* m_data = m._data();
  Num_T* r_data = r._data();
  for (int i = 0; i < r._datasize(); i++)
    r_data[i] = m_data[i] * t;

  return r;
}

template<class Num_T> inline
Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m)
{
  return operator*(m, t);
}

template<class Num_T> inline
Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  Mat<Num_T> out;
  elem_mult_out(m1, m2, out);
  return out;
}

template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   Mat<Num_T> &out)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::elem_mult_out(): Wrong sizes");
  out.set_size(m1.no_rows, m1.no_cols);
  for (int i = 0; i < out.datasize; i++)
    out.data[i] = m1.data[i] * m2.data[i];
}

template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   const Mat<Num_T> &m3, Mat<Num_T> &out)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows)
                  && (m1.no_cols == m2.no_cols) && (m1.no_cols == m3.no_cols),
                  "Mat<>::elem_mult_out(): Wrong sizes");
  out.set_size(m1.no_rows, m1.no_cols);
  for (int i = 0; i < out.datasize; i++)
    out.data[i] = m1.data[i] * m2.data[i] * m3.data[i];
}

template<class Num_T>
void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                   const Mat<Num_T> &m3, const Mat<Num_T> &m4,
                   Mat<Num_T> &out)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_rows == m3.no_rows)
                  && (m1.no_rows == m4.no_rows) && (m1.no_cols == m2.no_cols)
                  && (m1.no_cols == m3.no_cols) && (m1.no_cols == m4.no_cols),
                  "Mat<>::elem_mult_out(): Wrong sizes");
  out.set_size(m1.no_rows, m1.no_cols);
  for (int i = 0; i < out.datasize; i++)
    out.data[i] = m1.data[i] * m2.data[i] * m3.data[i] * m4.data[i];
}

template<class Num_T>
#ifndef _MSC_VER
inline
#endif
void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::elem_mult_inplace(): Wrong sizes");
  for (int i = 0; i < m2.datasize; i++)
    m2.data[i] *= m1.data[i];
}

template<class Num_T> inline
Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::elem_mult_sum(): Wrong sizes");
  Num_T acc = 0;

  for (int i = 0; i < m1.datasize; i++)
    acc += m1.data[i] * m2.data[i];

  return acc;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator/=(Num_T t)
{
  for (int i = 0; i < datasize; i++)
    data[i] /= t;
  return *this;
}

template<class Num_T> inline
Mat<Num_T>& Mat<Num_T>::operator/=(const Mat<Num_T> &m)
{
  it_assert_debug((m.no_rows == no_rows) && (m.no_cols == no_cols),
                  "Mat<>::operator/=(): Wrong sizes");
  for (int i = 0; i < datasize; i++)
    data[i] /= m.data[i];
  return *this;
}

template<class Num_T>
Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);
  for (int i = 0; i < r.datasize; ++i)
    r.data[i] = m.data[i] / t;
  return r;
}

template<class Num_T>
Mat<Num_T> operator/(Num_T t, const Mat<Num_T> &m)
{
  Mat<Num_T> r(m.no_rows, m.no_cols);
  for (int i = 0; i < r.datasize; ++i)
    r.data[i] = t / m.data[i];
  return r;
}

template<class Num_T> inline
Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  Mat<Num_T> out;
  elem_div_out(m1, m2, out);
  return out;
}

template<class Num_T>
void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                  Mat<Num_T> &out)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::elem_div_out(): Wrong sizes");

  if ((out.no_rows != m1.no_rows) || (out.no_cols != m1.no_cols))
    out.set_size(m1.no_rows, m1.no_cols);

  for (int i = 0; i < out.datasize; i++)
    out.data[i] = m1.data[i] / m2.data[i];
}

template<class Num_T> inline
Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
{
  it_assert_debug((m1.no_rows == m2.no_rows) && (m1.no_cols == m2.no_cols),
                  "Mat<>::elem_div_sum(): Wrong sizes");
  Num_T acc = 0;

  for (int i = 0; i < m1.datasize; i++)
    acc += m1.data[i] / m2.data[i];

  return acc;
}

template<class Num_T>
bool Mat<Num_T>::operator==(const Mat<Num_T> &m) const
{
  if (no_rows != m.no_rows || no_cols != m.no_cols) return false;
  for (int i = 0;i < datasize;i++) {
    if (data[i] != m.data[i]) return false;
  }
  return true;
}

template<class Num_T>
bool Mat<Num_T>::operator!=(const Mat<Num_T> &m) const
{
  if (no_rows != m.no_rows || no_cols != m.no_cols) return true;
  for (int i = 0;i < datasize;i++) {
    if (data[i] != m.data[i]) return true;
  }
  return false;
}

template <class Num_T>
std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m)
{
  int i;

  switch (m.rows()) {
  case 0 :
    os << "[]";
    break;
  case 1 :
    os << '[' << m.get_row(0) << ']';
    break;
  default:
    os << '[' << m.get_row(0) << std::endl;
    for (i = 1; i < m.rows() - 1; i++)
      os << ' ' << m.get_row(i) << std::endl;
    os << ' ' << m.get_row(m.rows() - 1) << ']';
  }

  return os;
}

template <class Num_T>
std::istream &operator>>(std::istream &is, Mat<Num_T> &m)
{
  std::ostringstream buffer;
  bool started = false;
  bool finished = false;
  bool brackets = false;
  bool within_double_brackets = false;
  char c;

  while (!finished) {
    if (is.eof()) {
      finished = true;
    }
    else {
      is.get(c);

      if (is.eof() || (c == '\n')) {
        if (brackets) {
          // Right bracket missing
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else if (!((c == '\n') && !started)) {
          finished = true;
        }
      }
      else if ((c == ' ') || (c == '\t')) {
        if (started) {
          buffer << ' ';
        }
      }
      else if (c == '[') {
        if ((started && !brackets) || within_double_brackets) {
          // Unexpected left bracket
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else if (!started) {
          started = true;
          brackets = true;
        }
        else {
          within_double_brackets = true;
        }
      }
      else if (c == ']') {
        if (!started || !brackets) {
          // Unexpected right bracket
          is.setstate(std::ios_base::failbit);
          finished = true;
        }
        else if (within_double_brackets) {
          within_double_brackets = false;
          buffer << ';';
        }
        else {
          finished = true;
        }
        while (!is.eof() && (((c = static_cast<char>(is.peek())) == ' ')
                             || (c == '\t'))) {
          is.get();
        }
        if (!is.eof() && (c == '\n')) {
          is.get();
        }
      }
      else {
        started = true;
        buffer << c;
      }
    }
  }

  if (!started) {
    m.set_size(0, false);
  }
  else {
    m.set(buffer.str());
  }

  return is;
}


// ---------------------------------------------------------------------
// Instantiations
// ---------------------------------------------------------------------

// class instantiations

ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Mat<double>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Mat<std::complex<double> >;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Mat<int>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Mat<short int>;
ITPP_EXPORT_TEMPLATE template class ITPP_EXPORT Mat<bin>;

// addition operators

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator+(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator+(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator+(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator+(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator+(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator+(const mat &m, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator+(const cmat &m, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator+(const imat &m, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator+(const smat &m, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator+(const bmat &m, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator+(double t, const mat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator+(std::complex<double> t, const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator+(int t, const imat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator+(short t, const smat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator+(bin t, const bmat &m);

// subtraction operators

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator-(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator-(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator-(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator-(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator-(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator-(const mat &m, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator-(const cmat &m, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator-(const imat &m, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator-(const smat &m, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator-(const bmat &m, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator-(double t, const mat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator-(std::complex<double> t, const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator-(int t, const imat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator-(short t, const smat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator-(bin t, const bmat &m);

// unary minus

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator-(const mat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator-(const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator-(const imat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator-(const smat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator-(const bmat &m);

// multiplication operators

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator*(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator*(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator*(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  ivec operator*(const imat &m, const ivec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  svec operator*(const smat &m, const svec &v);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bvec operator*(const bmat &m, const bvec &v);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator*(const mat &m, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator*(const cmat &m, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator*(const imat &m, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator*(const smat &m, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator*(const bmat &m, bin t);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator*(double t, const mat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator*(std::complex<double> t, const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator*(int t, const imat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator*(short t, const smat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator*(bin t, const bmat &m);

// element-wise multiplication

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat elem_mult(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat elem_mult(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat elem_mult(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat elem_mult(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat elem_mult(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const mat &m1, const mat &m2, mat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const cmat &m1, const cmat &m2,
                                     cmat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const imat &m1, const imat &m2,
                                     imat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const smat &m1, const smat &m2,
                                     smat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const bmat &m1, const bmat &m2,
                                     bmat &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const mat &m1, const mat &m2,
                                     const mat &m3, mat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const cmat &m1, const cmat &m2,
                                     const cmat &m3, cmat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const imat &m1, const imat &m2,
                                     const imat &m3, imat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const smat &m1, const smat &m2,
                                     const smat &m3, smat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const bmat &m1, const bmat &m2,
                                     const bmat &m3, bmat &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const mat &m1, const mat &m2,
                                     const mat &m3, const mat &m4, mat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const cmat &m1, const cmat &m2,
                                     const cmat &m3, const cmat &m4,
                                     cmat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const imat &m1, const imat &m2,
                                     const imat &m3, const imat &m4,
                                     imat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const smat &m1, const smat &m2,
                                     const smat &m3, const smat &m4,
                                     smat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_out(const bmat &m1, const bmat &m2,
                                     const bmat &m3, const bmat &m4,
                                     bmat &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_inplace(const mat &m1, mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_inplace(const cmat &m1, cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_inplace(const imat &m1, imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_inplace(const smat &m1, smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_mult_inplace(const bmat &m1, bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  double elem_mult_sum(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::complex<double> elem_mult_sum(const cmat &m1,
    const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  int elem_mult_sum(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  short elem_mult_sum(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bin elem_mult_sum(const bmat &m1, const bmat &m2);

// division operator

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator/(double t, const mat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator/(std::complex<double> t, const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator/(int t, const imat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator/(short t, const smat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator/(bin t, const bmat &m);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat operator/(const mat &m, double t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat operator/(const cmat &m, std::complex<double> t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat operator/(const imat &m, int t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat operator/(const smat &m, short t);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat operator/(const bmat &m, bin t);

// element-wise division

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat elem_div(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat elem_div(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat elem_div(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat elem_div(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat elem_div(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_div_out(const mat &m1, const mat &m2, mat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_div_out(const cmat &m1, const cmat &m2, cmat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_div_out(const imat &m1, const imat &m2, imat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_div_out(const smat &m1, const smat &m2, smat &out);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  void elem_div_out(const bmat &m1, const bmat &m2, bmat &out);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  double elem_div_sum(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::complex<double> elem_div_sum(const cmat &m1,
    const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  int elem_div_sum(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  short elem_div_sum(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bin elem_div_sum(const bmat &m1, const bmat &m2);

// concatenation

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat concat_horizontal(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat concat_horizontal(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat concat_horizontal(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat concat_horizontal(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat concat_horizontal(const bmat &m1, const bmat &m2);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  mat concat_vertical(const mat &m1, const mat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  cmat concat_vertical(const cmat &m1, const cmat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  imat concat_vertical(const imat &m1, const imat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  smat concat_vertical(const smat &m1, const smat &m2);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  bmat concat_vertical(const bmat &m1, const bmat &m2);

// I/O streams

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::ostream &operator<<(std::ostream &os, const mat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::ostream &operator<<(std::ostream &os, const cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::ostream &operator<<(std::ostream &os, const imat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::ostream &operator<<(std::ostream &os, const smat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::ostream &operator<<(std::ostream &os, const bmat  &m);

ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::istream &operator>>(std::istream &is, mat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::istream &operator>>(std::istream &is, cmat &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::istream &operator>>(std::istream &is, imat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::istream &operator>>(std::istream &is, smat  &m);
ITPP_EXPORT_TEMPLATE template ITPP_EXPORT  std::istream &operator>>(std::istream &is, bmat  &m);


} // namespace itpp

#endif // #ifndef MAT_H