/*!
 * \file
 * \brief Matrix Class Definitions
 * \author Tony Ottosson, Tobias Ringstrom, Adam Piatyszek and Conrad Sanderson
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2008  (see AUTHORS file for a list of contributors)
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 */

#ifndef MAT_H
#define MAT_H

#ifndef _MSC_VER
#  include <itpp/config.h>
#else
#  include <itpp/config_msvc.h>
#endif

#include <itpp/base/itassert.h>
#include <itpp/base/math/misc.h>
#include <itpp/base/factory.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;

  //! Horizontal concatenation of two matrices
  template<class Num_T>
  Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Vertical concatenation of two matrices
  template<class Num_T>
  Mat<Num_T> concat_vertical(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

  //! Addition of two matrices
  template<class Num_T>
  Mat<Num_T> operator+(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Addition of a matrix and a scalar
  template<class Num_T>
  Mat<Num_T> operator+(const Mat<Num_T> &m, Num_T t);
  //! Addition of a scalar and a matrix
  template<class Num_T>
  Mat<Num_T> operator+(Num_T t, const Mat<Num_T> &m);

  //! Subtraction of two matrices
  template<class Num_T>
  Mat<Num_T> operator-(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Subtraction of matrix and scalar
  template<class Num_T>
  Mat<Num_T> operator-(const Mat<Num_T> &m, Num_T t);
  //! Subtraction of scalar and matrix
  template<class Num_T>
  Mat<Num_T> operator-(Num_T t, const Mat<Num_T> &m);
  //! Negation of matrix
  template<class Num_T>
  Mat<Num_T> operator-(const Mat<Num_T> &m);

  //! Multiplication of two matrices
  template<class Num_T>
  Mat<Num_T> operator*(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Multiplication of matrix and vector
  template<class Num_T>
  Vec<Num_T> operator*(const Mat<Num_T> &m, const Vec<Num_T> &v);
  //! Multiplication of vector and matrix (matrix must be a row vector)
  template<class Num_T>
  Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m);
  //! Multiplication of matrix and scalar
  template<class Num_T>
  Mat<Num_T> operator*(const Mat<Num_T> &m, Num_T t);
  //! Multiplication of scalar and matrix
  template<class Num_T>
  Mat<Num_T> operator*(Num_T t, const Mat<Num_T> &m);

  //! Element wise multiplication of two matrices
  template<class Num_T>
  Mat<Num_T> elem_mult(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Element wise multiplication of two matrices, storing the result in matrix \c out
  template<class Num_T>
  void elem_mult_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                     Mat<Num_T> &out);
  //! Element wise multiplication of three matrices, storing the result in matrix \c 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);
  //! Element wise multiplication of four matrices, storing the result in matrix \c 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);
  //! In-place element wise multiplication of two matrices. Fast version of B = elem_mult(A, B).
  template<class Num_T>
  void elem_mult_inplace(const Mat<Num_T> &m1, Mat<Num_T> &m2);
  //! Element wise multiplication of two matrices, followed by summation of the resultant elements. Fast version of sumsum(elem_mult(A, B)).
  template<class Num_T>
  Num_T elem_mult_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

  //! Division of matrix and scalar
  template<class Num_T>
  Mat<Num_T> operator/(const Mat<Num_T> &m, Num_T t);

  //! Element wise division of two matrices
  template<class Num_T>
  Mat<Num_T> elem_div(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
  //! Element wise division of two matrices, storing the result in matrix \c out
  template<class Num_T>
  void elem_div_out(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                    Mat<Num_T> &out);
  //! Element wise division of two matrices, followed by summation of the resultant elements. Fast version of sumsum(elem_div(A, B)).
  template<class Num_T>
  Num_T elem_div_sum(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

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

  /*!
    \ingroup arr_vec_mat
    \brief Matrix Class (Templated)
    \author Tony Ottosson, Tobias Ringstrom, Adam Piatyszek and Conrad Sanderson

    Matrices can be of arbitrarily types, but conversions and functions are
    prepared for \c bin, \c short, \c int, \c double, and \c complex<double>
    vectors and these are predefined as: \c bmat, \c smat, \c imat, \c mat,
    and \c cmat. \c double and \c complex<double> are usually \c double and
    \c complex<double> respectively. However, this can be changed when
    compiling the it++ (see installation notes for more details). (Note: for
    binary matrices, an alternative to the bmat class is \c GF2mat and
    \c GF2mat_dense, which offer a more memory efficient representation and
    additional functions for linear algebra.)

    Examples:

    Matrix Constructors:
    When constructing a matrix without dimensions (memory) use
    \code mat temp; \endcode
    For construction of a matrix of a given size use
    \code mat temp(rows, cols); \endcode
    It is also possible to assign the constructed matrix the value and dimension
    of another matrix by
    \code vec temp(inmatrix); \endcode
    If you have explicit values you would like to assign to the matrix it is
    possible to do this using strings as:
    \code
    mat a("0 0.7;5 9.3"); // that is a = [0, 0.7; 5, 9.3]
    mat a="0 0.7;5 9.3";  // the constructor are called implicitly
    \endcode
    It is also possible to change dimension by
    \code temp.set_size(new_rows, new_cols, false); \endcode
    where \c false is used to indicate that the old values in \c temp
    is not copied. If you like to preserve the values use \c true.

    There are a number of methods to access parts of a matrix. Examples are
    \code
    a(5,3);     // Element number (5,3)
    a(5,9,3,5);  // Sub-matrix from rows 5, 6, 7, 8, 9 the columns 3, 4, and 5
    a.get_row(10);  // Row 10
    a.get_col(10); // Column 10
    \endcode

    It is also possible to modify parts of a vector as e.g. in
    \code
    a.set_row(5, invector);    // Set row 5 to \c invector
    a.set_col(3, invector); // Set column 3 to \c invector
    a.copy_col(1, 5); // Copy column 5 to column 1
    a.swap_cols(1, 5); // Swap the contents of columns 1 and 5
    \endcode

    It is of course also possible to perform the common linear algebra
    methods such as addition, subtraction, and matrix multiplication. Observe
    though, that vectors are assumed to be column-vectors in operations with
    matrices.

    Most elementary functions such as sin(), cosh(), log(), abs(), ..., are
    available as operations on the individual elements of the matrices. Please
    see the individual functions for more details.

    By default, the Mat elements are created using the default constructor for
    the element type. This can be changed by specifying a suitable Factory in
    the Mat constructor call; see Detailed Description for Factory.
  */
  template<class Num_T>
  class Mat {
  public:
    //! The type of the matrix values
    typedef Num_T value_type;

    //! Default constructor. An element factory \c f can be specified
    explicit Mat(const Factory &f = DEFAULT_FACTORY);
    //! Create a matrix of size (rows, cols). An element factory \c f can be specified.
    Mat(int rows, int cols, const Factory &f = DEFAULT_FACTORY);
    //! Copy constructor
    Mat(const Mat<Num_T> &m);
    //! Constructor, similar to the copy constructor, but also takes an element factory \c f as argument
    Mat(const Mat<Num_T> &m, const Factory &f);
    //! Construct a matrix from a column vector \c v. An element factory \c f can be specified.
    Mat(const Vec<Num_T> &v, const Factory &f = DEFAULT_FACTORY);
    //! Set matrix equal to values in string \c str. An element factory \c f can be specified.
    Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY);
    //! Set matrix equal to values in string \c str. An element factory \c f can be specified.
    Mat(const char *str, const Factory &f = DEFAULT_FACTORY);
    /*!
     * \brief Constructor taking a C-array as input. An element factory \c f
     * can be specified.
     *
     * By default the matrix is stored as a row-major matrix (i.e. listing
     * elements in sequence beginning with the first column).
     */
    Mat(const Num_T *c_array, int rows, int cols, bool row_major = true,
        const Factory &f = DEFAULT_FACTORY);

    //! Destructor
    ~Mat();

    //! The number of columns
    int cols() const { return no_cols; }
    //! The number of rows
    int rows() const { return no_rows; }
    //! The number of elements
    int size() const { return datasize; }
    //! Set size of matrix. If copy = true then keep the data before resizing.
    void set_size(int rows, int cols, bool copy = false);
    //! Set matrix equal to the all zero matrix
    void zeros();
    //! Set matrix equal to the all zero matrix
    void clear() { zeros(); }
    //! Set matrix equal to the all one matrix
    void ones();
    //! Set matrix equal to values in \c values
    void set(const char *str);
    //! Set matrix equal to values in the string \c str
    void set(const std::string &str);

    //! Get element (r,c) from matrix
    const Num_T &operator()(int r, int c) const;
    //! Get element (r,c) from matrix
    Num_T &operator()(int r, int c);
    //! Get element \c i using linear addressing (by rows)
    const Num_T &operator()(int i) const;
    //! Get element \c i using linear addressing (by rows)
    Num_T &operator()(int i);
    //! Get element (r,c) from matrix
    const Num_T &get(int r, int c) const;
    //! Set element (r,c) of matrix
    void set(int r, int c, Num_T t);

    /*!
      \brief Sub-matrix from row \c r1 to row \c r2 and columns \c c1 to \c c2.

      Value -1 indicates the last row and column, respectively.
    */
    Mat<Num_T> operator()(int r1, int r2, int c1, int c2) const;
    /*!
      \brief Sub-matrix from row \c r1 to row \c r2 and columns \c c1 to \c c2.

      Value -1 indicates the last row and column, respectively.
    */
    Mat<Num_T> get(int r1, int r2, int c1, int c2) const;

    //! Get row \c r
    Vec<Num_T> get_row(int r) const;
    //! Get rows \c r1 through \c r2
    Mat<Num_T> get_rows(int r1, int r2) const;
    //! Get the rows specified by \c indexlist
    Mat<Num_T> get_rows(const Vec<int> &indexlist) const;
    //! Get column \c c
    Vec<Num_T> get_col(int c) const;
    //! Get columns \c c1 through \c c2
    Mat<Num_T> get_cols(int c1, int c2) const;
    //! Get the columns specified by \c indexlist
    Mat<Num_T> get_cols(const Vec<int> &indexlist) const;
    //! Set row \c r to vector \c v
    void set_row(int r, const Vec<Num_T> &v);
    //! Set column \c c to vector \c v
    void set_col(int c, const Vec<Num_T> &v);
    //! Set rows to matrix \c m, staring from row \c r
    void set_rows(int r, const Mat<Num_T> &m);
    //! Set columns to matrix \c m, starting from column \c c
    void set_cols(int c, const Mat<Num_T> &m);
    //! Copy row \c from onto row \c to
    void copy_row(int to, int from);
    //! Copy column \c from onto column \c to
    void copy_col(int to, int from);
    //! Swap the rows \c r1 and \c r2
    void swap_rows(int r1, int r2);
    //! Swap the columns \c c1 and \c c2
    void swap_cols(int c1, int c2);

    //! Set submatrix defined by rows r1,r2 and columns c1,c2 to matrix m
    void set_submatrix(int r1, int r2, int c1, int c2, const Mat<Num_T> &m);
    //! Set submatrix defined by upper-left element (r,c) and the size of matrix m to m
    void set_submatrix(int r, int c, const Mat<Num_T> &m);
    //! Set all elements of submatrix defined by rows r1,r2 and columns c1,c2 to value t
    void set_submatrix(int r1, int r2, int c1, int c2, Num_T t);

    //! Delete row number \c r
    void del_row(int r);
    //! Delete rows from \c r1 to \c r2
    void del_rows(int r1, int r2);
    //! Delete column number \c c
    void del_col(int c);
    //! Delete columns from \c c1 to \c c2
    void del_cols(int c1, int c2);
    //! Insert vector \c v at row number \c r. The matrix can be empty.
    void ins_row(int r, const Vec<Num_T> &v);
    //! Insert vector \c v at column number \c c. The matrix can be empty.
    void ins_col(int c, const Vec<Num_T> &v);
    //! Append vector \c v to the bottom of the matrix. The matrix can be empty.
    void append_row(const Vec<Num_T> &v);
    //! Append vector \c v to the right side of the matrix. The matrix can be empty.
    void append_col(const Vec<Num_T> &v);

    //! Matrix transpose
    Mat<Num_T> transpose() const;
    //! Matrix transpose
    Mat<Num_T> T() const { return this->transpose(); }
    //! Hermitian matrix transpose (conjugate transpose)
    Mat<Num_T> hermitian_transpose() const;
    //! Hermitian matrix transpose (conjugate transpose)
    Mat<Num_T> H() const { return this->hermitian_transpose(); }

    //! Concatenate the matrices \c m1 and \c m2 horizontally
    friend Mat<Num_T> concat_horizontal<>(const Mat<Num_T> &m1,
                                          const Mat<Num_T> &m2);
    //! Concatenate the matrices \c m1 and \c m2 vertically
    friend Mat<Num_T> concat_vertical<>(const Mat<Num_T> &m1,
                                        const Mat<Num_T> &m2);

    //! Set all elements of the matrix equal to \c t
    Mat<Num_T>& operator=(Num_T t);
    //! Set matrix equal to \c m
    Mat<Num_T>& operator=(const Mat<Num_T> &m);
    //! Set matrix equal to the vector \c v, assuming column vector
    Mat<Num_T>& operator=(const Vec<Num_T> &v);
    //! Set matrix equal to values in the string
    Mat<Num_T>& operator=(const char *str);

    //! Addition of matrices
    Mat<Num_T>& operator+=(const Mat<Num_T> &m);
    //! Addition of scalar to matrix
    Mat<Num_T>& operator+=(Num_T t);
    //! Addition of two matrices
    friend Mat<Num_T> operator+<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
    //! Addition of matrix and scalar
    friend Mat<Num_T> operator+<>(const Mat<Num_T> &m, Num_T t);
    //! Addition of scalar and matrix
    friend Mat<Num_T> operator+<>(Num_T t, const Mat<Num_T> &m);

    //! Subtraction of matrix
    Mat<Num_T>& operator-=(const Mat<Num_T> &m);
    //! Subtraction of scalar from matrix
    Mat<Num_T>& operator-=(Num_T t);
    //! Subtraction of \c m2 from \c m1
    friend Mat<Num_T> operator-<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
    //! Subtraction of scalar from matrix
    friend Mat<Num_T> operator-<>(const Mat<Num_T> &m, Num_T t);
    //! Subtract matrix from scalar
    friend Mat<Num_T> operator-<>(Num_T t, const Mat<Num_T> &m);
    //! Subtraction of matrix
    friend Mat<Num_T> operator-<>(const Mat<Num_T> &m);

    //! Matrix multiplication
    Mat<Num_T>& operator*=(const Mat<Num_T> &m);
    //! Multiplication by a scalar
    Mat<Num_T>& operator*=(Num_T t);
    //! Multiplication of two matrices
    friend Mat<Num_T> operator*<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
    //! Multiplication of matrix \c m and vector \c v (column vector)
    friend Vec<Num_T> operator*<>(const Mat<Num_T> &m, const Vec<Num_T> &v);
    /*!
     * \brief Multiplication of vector \c v and matrix \c m with only one row
     *
     * This operator multiplies a column vector \c v times matrix \c m that
     * consists of only one row. Thus, the result of this operator is
     * exactly the same as the result of the outer product of two vectors,
     * i.e.: <tt>outer_product(v, m.get_col(0))</tt>.
     *
     * \note This operator is deprecated and might be removed or changed in
     * future releases of IT++.
     */
    friend Mat<Num_T> operator*<>(const Vec<Num_T> &v, const Mat<Num_T> &m);
    //! Multiplication of matrix and scalar
    friend Mat<Num_T> operator*<>(const Mat<Num_T> &m, Num_T t);
    //! Multiplication of scalar and matrix
    friend Mat<Num_T> operator*<>(Num_T t, const Mat<Num_T> &m);

    //! Element wise multiplication of two matrices
    friend Mat<Num_T> elem_mult<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
    //! Element wise multiplication of two matrices, storing the result in matrix \c out
    friend void elem_mult_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                                Mat<Num_T> &out);
    //! Element wise multiplication of three matrices, storing the result in matrix \c 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);
    //! Element wise multiplication of four matrices, storing the result in matrix \c 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);
    //! In-place element wise multiplication of two matrices. Fast version of B = elem_mult(A, B).
    friend void elem_mult_inplace<>(const Mat<Num_T> &m1, Mat<Num_T> &m2);
    //! Element wise multiplication of two matrices, followed by summation of the resultant elements. Fast version of sumsum(elem_mult(A, B)).
    friend Num_T elem_mult_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

    //! Division by a scalar
    Mat<Num_T>& operator/=(Num_T t);
    //! Division of matrix with scalar
    friend Mat<Num_T> operator/<>(const Mat<Num_T> &m, Num_T t);
    //! Element-wise division with the current matrix
    Mat<Num_T>& operator/=(const Mat<Num_T> &m);

    //! Element wise division of two matrices
    friend Mat<Num_T> elem_div<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);
    //! Element wise division of two matrices, storing the result in matrix \c out
    friend void elem_div_out<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2,
                               Mat<Num_T> &out);
    //! Element wise division of two matrices, followed by summation of the resultant elements. Fast version of sumsum(elem_div(A, B)).
    friend Num_T elem_div_sum<>(const Mat<Num_T> &m1, const Mat<Num_T> &m2);

    //! Compare two matrices. False if wrong sizes or different values
    bool operator==(const Mat<Num_T> &m) const;
    //! Compare two matrices. True if different
    bool operator!=(const Mat<Num_T> &m) const;

    //! Get element (r,c) from matrix without boundary check (not recommended to use)
    Num_T &_elem(int r, int c) { return data[r+c*no_rows]; }
    //! Get element (r,c) from matrix without boundary check (not recommended to use)
    const Num_T &_elem(int r, int c) const { return data[r+c*no_rows]; }
    //! Get element \c i using linear addressing (by rows) without boundary check (not recommended to use)
    Num_T &_elem(int i) { return data[i]; }
    //! Get element \c i using linear addressing (by rows) without boundary check (not recommended to use)
    const Num_T &_elem(int i) const { return data[i]; }

    //! Access of the internal data structure (not recommended to use)
    Num_T *_data() { return data; }
    //! Access to the internal data structure (not recommended to use)
    const Num_T *_data() const { return data; }
    //! Access to the internal data structure (not recommended to use)
    int _datasize() const { return datasize; }

  protected:
    //! Allocate memory for the matrix
    void alloc(int rows, int cols);
    //! Free the memory space of the matrix
    void free();

    /*! Protected integer variables
     * @{ */
    int datasize, no_rows, no_cols;
    /*! @} */
    //! Protected data pointer
    Num_T *data;
    //! Element factory (set to DEFAULT_FACTORY to use Num_T default constructors only)
    const Factory &factory;

  private:
    //! Check whether element (r,c) is within the matrix
    bool in_range(int r, int c) const
    {
      return ((r >=0) && (r < no_rows) && (c >= 0) && (c < no_cols));
    }
    //! Check whether row \c r is in the allowed range
    bool row_in_range(int r) const { return ((r >= 0) && (r < no_rows)); }
    //! Check whether column \c c is in the allowed range
    bool col_in_range(int c) const { return ((c >= 0) && (c < no_cols)); }
    //! Check whether element \c i is in the allowed range
    bool in_range(int i) const { return ((i >= 0) && (i < datasize)); }
  };

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

  /*!
    \relates Mat
    \brief Default Matrix Type
  */
  typedef Mat<double> mat;

  /*!
    \relates Mat
    \brief Default Complex Matrix Type
  */
  typedef Mat<std::complex<double> > cmat;

  /*!
    \relates Mat
    \brief Integer matrix
  */
  typedef Mat<int> imat;

  /*!
    \relates Mat
    \brief short int matrix
  */
  typedef Mat<short int> smat;

  /*!
    \relates Mat
    \relates GF2mat
    \relates GF2mat_sparse
    \brief bin matrix
  */
  typedef Mat<bin> bmat;

} //namespace itpp


#include <itpp/base/vec.h>

namespace itpp {

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

  /*!
    \relatesalso Mat
    \brief Output stream for matrices
  */
  template <class Num_T>
  std::ostream &operator<<(std::ostream &os, const Mat<Num_T> &m);

  /*!
    \relatesalso Mat
    \brief Input stream for matrices

    The input can be on the form "1 2 3; 4 5 6" or "[[1 2 3][4 5 6]]", i.e. with
    brackets or semicolons as row delimiters. The first form is compatible with
    the set method, while the second form is compatible with the ostream
    operator. The elements on a row can be separated by blank space or commas.
    Rows that are shorter than the longest row are padded with zero elements.
    "[]" means an empty matrix.
  */
  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(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;
    // 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
  {
    it_assert_debug(in_range(r, c), "Mat<>::get(): Indexing out of range");
    return data[r+c*no_rows];
  }

  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>
  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 r2, int c1, int c2,
                                 const Mat<Num_T> &m)
  {

    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 && r2>=0 && r1<no_rows && r2<no_rows &&
	       c1>=0 && c2>=0 && c1<no_cols && c2<no_cols, "Mat<Num_T>::set_submatrix(): index out of range");

    it_assert_debug(r2>=r1 && c2>=c1, "Mat<Num_T>::set_submatrix: r2<r1 or c2<c1");
    it_assert_debug(m.no_rows == r2-r1+1 && m.no_cols == c2-c1+1, "Mat<Num_T>::set_submatrix(): sizes don't match");

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



  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");

    int i, j, pos, rows = r2-r1+1;

    for (i=c1; i<=c2; i++) {
      pos = i*no_rows+r1;
      for (j=0; j<rows; 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;
  }

  //! \cond
  template<>
  cmat cmat::hermitian_transpose() const;
  //! \endcond

  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<class Num_T>
  Mat<Num_T> concat_horizontal(const Mat<Num_T> &m1, const Mat<Num_T> &m2)
  {
    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)
  {
    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 char *str)
  {
    set(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;
  }

#if defined(HAVE_BLAS)
  template<> mat& mat::operator*=(const mat &m);
  template<> cmat& cmat::operator*=(const cmat &m);
#endif

  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<class Num_T> inline
  Mat<Num_T>& Mat<Num_T>::operator*=(Num_T t)
  {
    scal_vector(datasize, t, data);
    return *this;
  }

#if defined(HAVE_BLAS)
  template<> mat operator*(const mat &m1, const mat &m2);
  template<> cmat operator*(const cmat &m1, const cmat &m2);
#endif


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

    Num_T tmp;
    int i, j, k;
    Num_T *tr=r.data, *t1, *t2=m2.data;

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

    return r;
  }

#if defined(HAVE_BLAS)
  template<> vec operator*(const mat &m, const vec &v);
  template<> cvec operator*(const cmat &m, const cvec &v);
#endif

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

    for (i=0; i<m.no_rows; i++) {
      r(i) = Num_T(0);
      m_pos = 0;
      for (k=0; k<m.no_cols; k++) {
	r(i) += m.data[m_pos+i] * v(k);
	m_pos += m.no_rows;
      }
    }

    return r;
  }

  template<class Num_T>
  Mat<Num_T> operator*(const Vec<Num_T> &v, const Mat<Num_T> &m)
  {
    it_assert((m.no_rows == 1),"Mat<Num_T>::operator*(): wrong sizes");
    it_warning("Mat<Num_T>::operator*(v, m): This operator is deprecated. "
               "Please use outer_product(v, m.get_row(0)) instead.");
    return outer_product(v, m.get_row(0));
  }

  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> 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> inline
  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>
  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> 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> 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 {
        c = is.get();

        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 = 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;
  }

  //! \cond

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

#ifdef HAVE_EXTERN_TEMPLATE

  // class instantiations

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

  // addition operators

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

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

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

  // subtraction operators

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

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

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

  // unary minus

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

  // multiplication operators

#if !defined(HAVE_BLAS)
  extern template mat operator*(const mat &m1, const mat &m2);
  extern template cmat operator*(const cmat &m1, const cmat &m2);
#endif
  extern template imat operator*(const imat &m1, const imat &m2);
  extern template smat operator*(const smat &m1, const smat &m2);
  extern template bmat operator*(const bmat &m1, const bmat &m2);

#if !defined(HAVE_BLAS)
  extern template vec operator*(const mat &m, const vec &v);
  extern template cvec operator*(const cmat &m, const cvec &v);
#endif
  extern template ivec operator*(const imat &m, const ivec &v);
  extern template svec operator*(const smat &m, const svec &v);
  extern template bvec operator*(const bmat &m, const bvec &v);

  extern template mat operator*(const vec &v, const mat &m);
  extern template cmat operator*(const cvec &v, const cmat &m);
  extern template imat operator*(const ivec &v, const imat &m);
  extern template smat operator*(const svec &v, const smat &m);
  extern template bmat operator*(const bvec &v, const bmat &m);

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

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

  // element-wise multiplication

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

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

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

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

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

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

  // division operator

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

  // element-wise division

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

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

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

  // concatenation

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

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

  // I/O streams

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

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

#endif // HAVE_EXTERN_TEMPLATE

  //! \endcond

} // namespace itpp

#endif // #ifndef MAT_H