#ifndef ICL_DYN_MATRIX_H
#define ICL_DYN_MATRIX_H
#include <iterator>
#include <algorithm>
#include <numeric>
#include <functional>
#include <iostream>
#include <vector>
#include <cmath>

#include <iclException.h>
#include <iclMacros.h>

#ifdef HAVE_IPP
#include <ipp.h>
#endif


namespace icl{

  /// Special linear algebra exception type  \ingroup LINALG \ingroup EXCEPT
  struct InvalidMatrixDimensionException :public ICLException{
    InvalidMatrixDimensionException(const std::string &msg):ICLException(msg){}
  };

  /// Special linear algebra exception type  \ingroup LINALG \ingroup EXCEPT
  struct IncompatibleMatrixDimensionException :public ICLException{
    IncompatibleMatrixDimensionException(const std::string &msg):ICLException(msg){}
  };

  /// Special linear algebra exception type  \ingroup LINALG \ingroup EXCEPT
  struct  InvalidIndexException : public ICLException{
    InvalidIndexException(const std::string &msg):ICLException(msg){}
  };

  /// Special linear algebra exception type  \ingroup LINALG \ingroup EXCEPT
  struct SingularMatrixException : public ICLException{
    SingularMatrixException(const std::string &msg):ICLException(msg){}
  };

  /// Special linear algebra exception type  \ingroup LINALG \ingroup EXCEPT
  struct QRDecompException : public ICLException{
    QRDecompException(const std::string &msg):ICLException(msg){}
  };
  

  /// Highly flexible and optimized matrix class implementation  \ingroup LINALG
  /** In contrast to the FixedMatrix template class, the DynMatrix instances are dynamically sized at runtime
      The template class is instantiated for the common ICL types
      uint8_t, int16_t, int32_t, float and double
  */
  template<class T>
  struct DynMatrix{
    
    /// Default empty constructor creates a null-matrix
    inline DynMatrix():m_rows(0),m_cols(0),m_data(0),m_ownData(true){}
    
    /// Create a dyn matrix with given dimensions (and optional initialValue)
    inline DynMatrix(unsigned int cols,unsigned int rows,const  T &initValue=0) throw (InvalidMatrixDimensionException) : 
    m_rows(rows),m_cols(cols),m_ownData(true){
      if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0");
      m_data = new T[cols*rows];
      std::fill(begin(),end(),initValue);
    }

    /// Create a matrix with given data
    /** Data can be wrapped deeply or shallowly. If the latter is true, given data pointer
        will not be released in the destructor*/
    inline DynMatrix(unsigned int cols,unsigned int rows, T *data, bool deepCopy=true) throw (InvalidMatrixDimensionException) : 
      m_rows(rows),m_cols(cols),m_ownData(deepCopy){
      if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0");
      if(deepCopy){
        m_data = new T[dim()];
        std::copy(data,data+dim(),begin());
      }else{
        m_data = data;
      }
    }

    /// Create a matrix with given data (const version: deepCopy only)
    inline DynMatrix(unsigned int cols,unsigned int rows,const T *data) throw (InvalidMatrixDimensionException) : 
      m_rows(rows),m_cols(cols),m_ownData(true){
      if(!dim()) throw InvalidMatrixDimensionException("matrix dimensions must be > 0");
      m_data = new T[dim()];
      std::copy(data,data+dim(),begin());
    }

    /// Default copy constructor
    inline DynMatrix(const DynMatrix &other):
    m_rows(other.m_rows),m_cols(other.m_cols),m_data(new T[dim()]),m_ownData(true){
      std::copy(other.begin(),other.end(),begin());
    }
    
    /// returns with this matrix has a valid data pointer
    inline bool isNull() const { return !m_data; }

    /// Destructor (deletes data if no wrapped shallowly)
    inline ~DynMatrix(){
      if(m_data && m_ownData) delete [] m_data;
    }

    /// Assignment operator (using deep-copy)
    inline DynMatrix &operator=(const DynMatrix &other){
      if(dim() != other.dim()){
        delete m_data;
        m_data = new T[other.dim()];
      }
      m_cols = other.m_cols;
      m_rows = other.m_rows;

      std::copy(other.begin(),other.end(),begin());
      return *this;
    }

    /// resets matrix dimensions  
    inline void setBounds(unsigned int cols, unsigned int rows, bool holdContent=false, const T &initializer=0) throw (InvalidMatrixDimensionException){
      if((int)cols == m_cols && (int)rows==m_rows) return;
      if(!(cols*rows)) throw InvalidMatrixDimensionException("matrix dimensions must be > 0");
      DynMatrix M(cols,rows,initializer);
      if(holdContent){
        unsigned int min_cols = iclMin(cols,(unsigned int)m_cols);
        unsigned int min_rows = iclMin(rows,(unsigned int)m_rows);
        for(unsigned int i=0;i<min_cols;++i){
          for(unsigned int j=0;j<min_rows;++j){
            M(i,j) = (*this)(i,j);
          }
        }
      }
      m_cols = cols;
      m_rows = rows;
      ICL_DELETE(m_data);
      m_data = M.begin();
      M.set_data(0);
    }
  
    /// tests weather a matrix is enough similar to another matrix
    inline bool isSimilar(const DynMatrix &other, T tollerance=0.0001) const{
      if(other.cols() != cols() || other.rows() != rows()) return false;
      for(unsigned int i=0;i<dim();++i){
        T diff = m_data[i] - other.m_data[i];
        if((diff?diff:-diff) > tollerance) return false;
      }
      return true;
    }
    
    /// elementwise comparison (==)
    inline bool operator==(const DynMatrix &other) const{
      if(other.cols() != cols() || other.rows() != rows()) return false;
      for(unsigned int i=0;i<dim();++i){
        if(m_data[i] !=  other.m_data[i]) return false;
      }
      return true;
    }

    /// elementwise comparison (!=)
    inline bool operator!=(const DynMatrix &other) const{
      if(other.cols() != cols() || other.rows() != rows()) return false;
      for(unsigned int i=0;i<dim();++i){
        if(m_data[i] !=  other.m_data[i]) return true;
      }
      return false;
    }


    /// Multiply elements with scalar
    inline DynMatrix operator*(T f) const{
      DynMatrix dst(cols(),rows());
      return mult(f,dst);
    }

    /// Multiply elements with scalar (in source destination fashion)
    inline DynMatrix &mult(T f, DynMatrix &dst) const{
      dst.setBounds(cols(),rows());
      std::transform(begin(),end(),dst.begin(),std::bind2nd(std::multiplies<T>(),f));
      return dst;      
    }

    /// Multiply elements with scalar (inplace)
    inline DynMatrix &operator*=(T f){
      std::transform(begin(),end(),begin(),std::bind2nd(std::multiplies<T>(),f));
      return *this;
    }

    /// Device elements by scalar
    inline DynMatrix operator/(T f) const{
      return this->operator*(1/f);
    }

    /// Device elements by scalar (inplace)
    inline DynMatrix &operator/=(T f){
      return this->operator*=(1/f);
    }

    /// Matrix multiplication (in source destination fashion) [IPP-Supported]
    inline DynMatrix &mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){
      if(m.rows() != cols()) throw IncompatibleMatrixDimensionException("A*B : rows(A) must be cols(B)");
      dst.setBounds(m.cols(),rows());
      for(unsigned int c=0;c<dst.cols();++c){
        for(unsigned int r=0;r<dst.rows();++r){
          dst(c,r) = std::inner_product(row_begin(r),row_end(r),m.col_begin(c),T(0));
        }
      }
      return dst;
    }

    /// Elementwise matrix multiplication (in source destination fashion) [IPP-Supported]
    inline DynMatrix &elementwise_mult(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){
      if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A.*B dimension mismatch");
      dst.setBounds(cols(),rows());
      for(int i=0;i<dim();++i){
	dst[i] = m_data[i] * m[i];
      }
      return dst;
    }

    /// Elementwise matrix multiplication (without destination matrix) [IPP-Supported]
    inline DynMatrix elementwise_mult(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){
      DynMatrix dst(cols(),rows());
      return elementwise_mult(m,dst);
    }

    /// Elementwise division (in source destination fashion) [IPP-Supported]
    inline DynMatrix &elementwise_div(const DynMatrix &m, DynMatrix &dst) const throw (IncompatibleMatrixDimensionException){
      if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A./B dimension mismatch");
      dst.setBounds(cols(),rows());
      for(int i=0;i<dim();++i){
	dst[i] = m_data[i] / m[i];
      }
      return dst;
    }

    /// Elementwise matrix multiplication (without destination matrix) [IPP-Supported]
    inline DynMatrix elementwise_div(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){
      DynMatrix dst(cols(),rows());
      return elementwise_div(m,dst);
    }

    
    
    
    /// Essential matrix multiplication [IPP-Supported]
    inline DynMatrix operator*(const DynMatrix &m) const throw (IncompatibleMatrixDimensionException){
      DynMatrix d(m.cols(),rows());
      return mult(m,d);
    }
    
    /// inplace matrix multiplication applying this = this*m [IPP-Supported]
    inline DynMatrix &operator*=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){
      return *this=((*this)*m);
    }
    
    /// inplace matrix devision (calling this/m.inv()) [IPP-Supported]
    inline DynMatrix operator/(const DynMatrix &m) const 
      throw (IncompatibleMatrixDimensionException,
             InvalidMatrixDimensionException,
             SingularMatrixException){
      return this->operator*(m.inv());
    }

    /// inplace matrix devision (calling this/m.inv()) (inplace)
    inline DynMatrix &operator/=(const DynMatrix &m) const 
      throw (IncompatibleMatrixDimensionException,
             InvalidMatrixDimensionException,
             SingularMatrixException){
      return *this = this->operator*(m.inv());
    }

    /// adds a scalar to each element
    inline DynMatrix operator+(const T &t){
      DynMatrix d(cols(),rows());
      std::transform(begin(),end(),d.begin(),std::bind2nd(std::plus<T>(),t));
      return d;
    }

    /// substacts a scalar from each element
    inline DynMatrix operator-(const T &t){
      DynMatrix d(cols(),rows());
      std::transform(begin(),end(),d.begin(),std::bind2nd(std::minus<T>(),t));
      return d;
    }

    /// adds a scalar to each element (inplace)
    inline DynMatrix &operator+=(const T &t){
      std::transform(begin(),end(),begin(),std::bind2nd(std::plus<T>(),t));
      return *this;
    }

    /// substacts a scalar from each element (inplace)
    inline DynMatrix &operator-=(const T &t){
      std::transform(begin(),end(),begin(),std::bind2nd(std::minus<T>(),t));
      return *this;
    }

    /// Matrix addition
    inline DynMatrix operator+(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){
      if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)");
      DynMatrix d(cols(),rows());
      std::transform(begin(),end(),m.begin(),d.begin(),std::plus<T>());
      return d;
    }

    /// Matrix substraction
    inline DynMatrix operator-(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){
      if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)");
      DynMatrix d(cols(),rows());
      std::transform(begin(),end(),m.begin(),d.begin(),std::minus<T>());
      return d;
    }

    /// Matrix addition (inplace)
    inline DynMatrix &operator+=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){
      if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)");
      std::transform(begin(),end(),m.begin(),begin(),std::plus<T>());
      return *this;
    }

    /// Matrix substraction (inplace)
    inline DynMatrix &operator-=(const DynMatrix &m) throw (IncompatibleMatrixDimensionException){
      if(cols() != m.cols() || rows() != m.rows()) throw IncompatibleMatrixDimensionException("A+B size(A) must be size(B)");
      std::transform(begin(),end(),m.begin(),begin(),std::minus<T>());
      return *this;
    }
  
    /// element access operator (x,y)-access index begin 0!
    inline T &operator()(unsigned int col,unsigned int row){
#ifdef DYN_MATRIX_INDEX_CHECK
      if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")");
#endif
      return m_data[col+cols()*row];
    }

    /// element access operator (x,y)-access index begin 0! (const)
    inline const T &operator() (unsigned int col,unsigned int row) const{
#ifdef DYN_MATRIX_INDEX_CHECK
      if((int)col >= m_cols || (int)row >= m_rows) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index (" << col << "," << row << ")");
#endif
      return m_data[col+cols()*row];
    }

    /// element access with index check
    inline T &at(unsigned int col,unsigned int row) throw (InvalidIndexException){
      if(col>=cols() || row >= rows()) throw InvalidIndexException("row or col index too large");
      return m_data[col+cols()*row];
    }

    /// element access with index check (const)
    inline const T &at(unsigned int col,unsigned int row) const throw (InvalidIndexException){
      return const_cast<DynMatrix*>(this)->at(col,row);
    }
    

    
    /// linear access to actual data array
    inline T &operator[](unsigned int idx) { 
      idx_check(idx);
      if(idx >= dim()) ERROR_LOG("access to "<<m_cols<<'x'<<m_rows<<"-matrix index [" << idx<< "]");

      return m_data[idx]; 
      
    }
    

    /// linear access to actual data array (const)
    inline const T &operator[](unsigned int idx) const { 
      idx_check(idx);
      return m_data[idx]; 
    }
    
    /// applies an L_l norm on the matrix elements (all elements are treated as vector)
    inline T norm(double l=2) const{
      double accu = 0;
      for(unsigned int i=0;i<dim();++i){
        accu += ::pow(double(m_data[i]),l);
      }
      return ::pow(double(accu),1.0/l);
    }

    

    /// default iterator type (just a data-pointer) 
    typedef T* iterator;
    
    /// dafault const_iterator type (just a data-pointer)
    typedef const T* const_iterator;

    /// comples row_iterator type
    typedef T* row_iterator;
    
    /// complex const_row_iterator type
    typedef const T* const_row_iterator;

    /// height of the matrix (number of rows)
    unsigned int rows() const { return m_rows; }
  
    /// width of the matrix (number of columns)
    unsigned int cols() const { return m_cols; }

    /// internal data pointer
    T *data() { return m_data; }

    /// internal data pointer (const)
    const T *data() const { return m_data; }

    /// matrix dimension (width*height) or (cols*rows)
    unsigned int dim() const { return m_rows*m_cols; }

    /// returns sizeof (T)*dim()
    int stride0() const { return sizeof(T) * dim(); }
    
    /// returns sizeof(T)*cols()
    int stride1() const { return sizeof(T) * cols(); }

    /// returns sizeof (T)
    int stride2() const { return sizeof(T); }
  
    /// Internal column iterator struct (using height-stride) \ingroup LINALG
    struct col_iterator : public std::iterator<std::random_access_iterator_tag,T>{
      typedef unsigned int difference_type;
      mutable T *p;
      unsigned int stride;
      inline col_iterator(T *col_begin,unsigned int stride):p(col_begin),stride(stride){}


    /// prefix increment operator
      inline col_iterator &operator++(){
        p+=stride;
        return *this;
      }
      /// prefix increment operator (const)
      inline const col_iterator &operator++() const{
        p+=stride;
        return *this;
      }
      /// postfix increment operator
      inline col_iterator operator++(int){
        col_iterator tmp = *this;
        ++(*this);
        return tmp;
      }
      /// postfix increment operator (const)
      inline const col_iterator operator++(int) const{
        col_iterator tmp = *this;
        ++(*this);
        return tmp;
      }

      /// prefix decrement operator
      inline col_iterator &operator--(){
        p-=stride;
        return *this;
      }

      /// prefix decrement operator (const)
      inline const col_iterator &operator--() const{
        p-=stride;
        return *this;
      }

      /// postfix decrement operator
      inline col_iterator operator--(int){
        col_iterator tmp = *this;
        --(*this);
        return tmp;
      }

      /// postfix decrement operator (const)
      inline const col_iterator operator--(int) const{
        col_iterator tmp = *this;
        --(*this);
        return tmp;
      }

      /// jump next n elements (inplace)
      inline col_iterator &operator+=(difference_type n){
        p += n * stride;
        return *this;
      }

      /// jump next n elements (inplace) (const)
      inline const col_iterator &operator+=(difference_type n) const{
        p += n * stride;
        return *this;
      }


      /// jump backward next n elements (inplace)
      inline col_iterator &operator-=(difference_type n){
        p -= n * stride;
        return *this;
      }

      /// jump backward next n elements (inplace) (const)
      inline const col_iterator &operator-=(difference_type n) const{
        p -= n * stride;
        return *this;
      }


      /// jump next n elements
      inline col_iterator operator+(difference_type n) {
        col_iterator tmp = *this;
        tmp+=n;
        return tmp;
      }

      /// jump next n elements (const)
      inline const col_iterator operator+(difference_type n) const{
        col_iterator tmp = *this;
        tmp+=n;
        return tmp;
      }

      /// jump backward next n elements
      inline col_iterator operator-(difference_type n) {
        col_iterator tmp = *this;
        tmp-=n;
        return tmp;
      }

      /// jump backward next n elements (const)
      inline const col_iterator operator-(difference_type n) const {
        col_iterator tmp = *this;
        tmp-=n;
        return tmp;
      }

      inline difference_type operator-(const col_iterator &other) const{
        return (p-other.p)/stride;
      }

      
      /// Dereference operator
      inline T &operator*(){
        return *p;
      }

      /// const Dereference operator
      inline T operator*() const{
        return *p;
      }

      /// comparison operator ==
      inline bool operator==(const col_iterator &i) const{ return p == i.p; }

      /// comparison operator !=
      inline bool operator!=(const col_iterator &i) const{ return p != i.p; }

      /// comparison operator <
      inline bool operator<(const col_iterator &i) const{ return p < i.p; }

      /// comparison operator <=
      inline bool operator<=(const col_iterator &i) const{ return p <= i.p; }

      /// comparison operator >=
      inline bool operator>=(const col_iterator &i) const{ return p >= i.p; }

      /// comparison operator >
      inline bool operator>(const col_iterator &i) const{ return p > i.p; }
    };
    
    /// const column iterator typedef
    typedef const col_iterator const_col_iterator;

    /// Internally used Utility structure referencing a matrix column shallowly
    struct DynMatrixColumn{
#ifdef DYN_MATRIX_INDEX_CHECK
#define DYN_MATRIX_COLUMN_CHECK(C,E) if(C) ERROR_LOG(E)
#else
#define DYN_MATRIX_COLUMN_CHECK(C,E)
#endif
      /// Matrix reference
      DynMatrix<T> *matrix;
      
      /// referenced column in matrix
      unsigned int column;
      
      /// create from source matrix and column index
      inline DynMatrixColumn(const DynMatrix<T> *matrix, unsigned int column):
      matrix(const_cast<DynMatrix<T>*>(matrix)),column(column){
        DYN_MATRIX_COLUMN_CHECK(column >= matrix->cols(),"invalid column index");
      }
      
      /// Create from source matrix (only works if matrix has only single column = column-vector)
      inline DynMatrixColumn(const DynMatrix<T> &matrix):
      matrix(const_cast<DynMatrix<T>*>(&matrix)),column(0){
        DYN_MATRIX_COLUMN_CHECK(matrix->cols() != 1,"source matrix must have exactly ONE column");
      }
      /// Shallow copy from another matrix column reference
      inline DynMatrixColumn(const DynMatrixColumn &c):
      matrix(c.matrix),column(c.column){}
      
      /// returns column begin
      inline col_iterator begin() { return matrix->col_begin(column); }
      
      /// returns column end
      inline col_iterator end() { return matrix->col_end(column); }

      /// returns column begin (const)
      inline const col_iterator begin() const { return matrix->col_begin(column); }

      /// returns column end (const)
      inline const col_iterator end() const { return matrix->col_end(column); }

      /// returns column length (matrix->rows())
      inline unsigned int dim() const { return matrix->rows(); }
      
      /// assignment by another column
      inline DynMatrixColumn &operator=(const DynMatrixColumn &c){
        DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch");
        std::copy(c.begin(),c.end(),begin());
        return *this;
      }
        
      /// assigne dyn matrix to matrix columns
      inline DynMatrixColumn &operator=(const DynMatrix &src){
        DYN_MATRIX_COLUMN_CHECK(dim() != src.dim(),"dimension missmatch");
        std::copy(src.begin(),src.end(),begin());
        return *this;
      }   

      /// operator += for other columns
      inline DynMatrixColumn &operator+=(const DynMatrixColumn &c){
        DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch");
        std::transform(c.begin(),c.end(),begin(),begin(),std::plus<T>());
        return *this;
      }   
      
      /// operator += for other columns
      inline DynMatrixColumn &operator-=(const DynMatrixColumn &c){
        DYN_MATRIX_COLUMN_CHECK(dim() != c.dim(),"dimension missmatch");
        std::transform(c.begin(),c.end(),begin(),begin(),std::minus<T>());
        return *this;
      }   

      /// operator += for DynMatrices
      inline DynMatrixColumn &operator+=(const DynMatrix &m){
        DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch");
        std::transform(m.begin(),m.end(),begin(),begin(),std::plus<T>());
        return *this;
      }   
      /// operator -= for DynMatrices
      inline DynMatrixColumn &operator-=(const DynMatrix &m){
        DYN_MATRIX_COLUMN_CHECK(dim() != m.dim(),"dimension missmatch");
        std::transform(m.begin(),m.end(),begin(),begin(),std::minus<T>());
        return *this;
      }   
      
      /// operator *= for scalars
      inline DynMatrixColumn &operator*=(const T&val){
        std::for_each(begin(),end(),std::bind2nd(std::multiplies<T>(),val));
        return *this;
      }
      /// operator /= for scalars
      inline DynMatrixColumn &operator/=(const T&val){
        std::for_each(begin(),end(),std::bind2nd(std::divides<T>(),val));
        return *this;
      }
      
    };

    inline DynMatrix &operator=(const DynMatrixColumn &col){
      DYN_MATRIX_COLUMN_CHECK(dim() != col.dim(),"dimension missmatch");
      std::copy(col.begin(),col.end(),begin());
      return *this;
    }

#undef DYN_MATRIX_COLUMN_CHECK
    
    

    /// returns an iterator to the begin of internal data array
    inline iterator begin() { return m_data; }

    /// returns an iterator to the end of internal data array
    inline iterator end() { return m_data+dim(); }
  
    /// returns an iterator to the begin of internal data array (const)
    inline const_iterator begin() const { return m_data; }

    /// returns an iterator to the end of internal data array (const)
    inline const_iterator end() const { return m_data+dim(); }
  
    /// returns an iterator running through a certain matrix column 
    inline col_iterator col_begin(unsigned int col) {       
      col_check(col);
      return col_iterator(m_data+col,cols()); 
    }

    /// returns an iterator end of a certain matrix column 
    inline col_iterator col_end(unsigned int col) {
      col_check(col);
      return col_iterator(m_data+col+dim(),cols()); 
    }
  
    /// returns an iterator running through a certain matrix column (const)
    inline const_col_iterator col_begin(unsigned int col) const { 
      col_check(col);
      return col_iterator(m_data+col,cols()); 
    }

    /// returns an iterator end of a certain matrix column (const)
    inline const_col_iterator col_end(unsigned int col) const { 
      col_check(col);
      return col_iterator(m_data+col+dim(),cols()); 
    }

    /// returns an iterator running through a certain matrix row 
    inline row_iterator row_begin(unsigned int row) { 
      row_check(row);
      return m_data+row*cols(); 
    }

    /// returns an iterator of a certains row's end
    inline row_iterator row_end(unsigned int row) { 
      row_check(row);
      return m_data+(row+1)*cols(); 
    }

    /// returns an iterator running through a certain matrix row  (const)
    inline const_row_iterator row_begin(unsigned int row) const { 
      row_check(row);
      return m_data+row*cols(); 
    }

    /// returns an iterator of a certains row's end (const)
    inline const_row_iterator row_end(unsigned int row) const {
      row_check(row);
      return m_data+(row+1)*cols(); 
    }
    
    /// Extracts a shallow copied matrix row
    inline DynMatrix row(int row){
      row_check(row);
      return DynMatrix(m_cols,1,row_begin(row),false);
    }

    /// Extracts a shallow copied matrix row (const)
    inline const DynMatrix row(int row) const{
      row_check(row);
      return DynMatrix(m_cols,1,const_cast<T*>(row_begin(row)),false);
    }

    /// Extracts a shallow copied matrix column
    inline DynMatrixColumn col(int col){
      return DynMatrixColumn(this,col);
    }

    inline const DynMatrixColumn col(int col) const{
      return DynMatrixColumn(this,col);
    }

    /// applies QR-decomposition (only for icl32f and icl64f)
    void decompose_QR(DynMatrix &Q, DynMatrix &R) const throw (QRDecompException);
    
    /// invert the matrix (only for icl32f and icl64f)
    DynMatrix inv() const throw (InvalidMatrixDimensionException,SingularMatrixException);
    
    /// calculates the Moore-Penrose pseudo-inverse (only implemented with IPP_OPTIMIZATION and only for icl32f and icl64f)
    /** Internally, this functions uses a QR-decomposition based approach, which is much more stable than
        the naiv approach pinv(X) * X*(X*X')^(-1)
        \code
        DynMatrix Q,R;
        decompose_QR(Q,R);
        return R.inv() * Q.transp();
        \endcode
    */
    DynMatrix pinv() const throw (InvalidMatrixDimensionException,SingularMatrixException,QRDecompException);


    /// matrix determinant (only for icl32f and icl64f)
    T det() const throw (InvalidMatrixDimensionException);

    /// matrix transposed
    inline DynMatrix transp() const{
      DynMatrix d(rows(),cols());
      for(unsigned int x=0;x<cols();++x){
        for(unsigned int y=0;y<rows();++y){
          d(y,x) = (*this)(x,y);
        }
      }
      return d;
    }

    /// inner product of data pointers (not matrix-mulitiplication)
    /** computes the inner-product of data vectors */
    T inner_product(const DynMatrix<T> &other) const {
      return std::inner_product(begin(),end(),other.begin(),T(0));
    }
    
    /// returns diagonal-elements as column-vector
    DynMatrix<T> diag() const{
      ICLASSERT_RETURN_VAL(cols()==rows(),DynMatrix<T>());
      DynMatrix<T> d(1,rows());
      for(int i=0;i<rows();++i){
        d[i] = (*this)(i,i);
      }
      return d;
    }
    
    /// computes the sum of all diagonal elements
    T trace() const{
      ICLASSERT_RETURN_VAL(cols()==rows(),0);
      double accu = 0;
      for(int i=0;i<dim();i+=cols()+1){
        accu += m_data[i];
      }
      return accu;
    }
    
    /// sets new data internally and returns old data pointer (for experts only!)
    inline T *set_data(T *newData){
      T *old_data = m_data;
      m_data = newData;
      return old_data;
    }
    
    /// creates a dim-D identity Matrix
    static inline DynMatrix id(unsigned int dim) {
      DynMatrix M(dim,dim);
      for(unsigned int i=0;i<dim;++i) M(i,i) = 1;
      return M;
    }
    
  private:
    inline void row_check(unsigned int row) const{
#ifdef DYN_MATRIX_INDEX_CHECK
      if((int)row >= m_rows) ERROR_LOG("access to row index " << row << " on a "<<m_cols<<'x'<<m_rows<<"-matrix");
#else
      (void)row;
#endif
    }
    inline void col_check(unsigned int col) const{
#ifdef DYN_MATRIX_INDEX_CHECK
      if((int)col >= m_cols) ERROR_LOG("access to column index " << col << " on a "<<m_cols<<'x'<<m_rows<<"-matrix");
#else
      (void)col;
#endif
    }
    inline void idx_check(unsigned int col, unsigned int row) const{
      col_check(col);
      row_check(row);
    }
    
    inline void idx_check(unsigned int idx) const{
#ifdef DYN_MATRIX_INDEX_CHECK
      if(idx >= dim()) ERROR_LOG("access to linear index " << idx << " on a "<<m_cols<<'x'<<m_rows<<"-matrix");
#else
      (void)idx;
#endif
    }

    int m_rows;
    int m_cols;
    T *m_data;
    bool m_ownData;
  };

  /// ostream operator implemented for uchar, short, int, float and double matrices  \ingroup LINALG
  template<class T>
  std::ostream &operator<<(std::ostream &s,const DynMatrix<T> &m);

  /// istream operator implemented for uchar, short, int, float and double matrices  \ingroup LINALG
  template<class T>
  std::istream &operator>>(std::istream &s,DynMatrix<T> &m);


#ifdef HAVE_IPP
  /** \cond */
  
#define DYN_MATRIX_MULT_SPECIALIZE(IPPT)                                                                \
  template<>	    							                                \
    inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult(                                            \
      const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const	                                \
    throw (IncompatibleMatrixDimensionException){			                                \
    if(m.rows() != cols()) throw IncompatibleMatrixDimensionException("A*B : rows(A) must be cols(B)");	\
    dst.setBounds(m.cols(),rows());					                                \
    ippmMul_mm_##IPPT(data(),sizeof(Ipp##IPPT)*cols(),sizeof(Ipp##IPPT),cols(),rows(),                  \
		      m.data(),sizeof(Ipp##IPPT)*m.cols(),sizeof(Ipp##IPPT),m.cols(),m.rows(),          \
		      dst.data(),m.cols()*sizeof(Ipp##IPPT),sizeof(Ipp##IPPT));	                        \
    return dst;					\
  }

  DYN_MATRIX_MULT_SPECIALIZE(32f)
  DYN_MATRIX_MULT_SPECIALIZE(64f)
#undef DYN_MATRIX_MULT_SPECIALIZE


#define DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(IPPT)	                      \
    template<>								      \
    inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::elementwise_div(      \
          const DynMatrix<Ipp##IPPT> &m, DynMatrix<Ipp##IPPT> &dst) const     \
      throw (IncompatibleMatrixDimensionException){                           \
    if((m.cols() != cols()) || (m.rows() != rows())){                         \
      throw IncompatibleMatrixDimensionException("A./B dimension mismatch");  \
    }                                                                         \
    dst.setBounds(cols(),rows());                                             \
    ippsDiv_##IPPT(data(),m.data(),dst.data(),dim());                         \
    return dst;                                                               \
  }

  DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(32f)
  DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE(64f)

#undef DYN_MATRIX_ELEMENT_WISE_DIV_SPECIALIZE






#define DYN_MATRIX_MULT_BY_CONSTANT(IPPT)	             \
    template<>						     \
    inline DynMatrix<Ipp##IPPT> &DynMatrix<Ipp##IPPT>::mult( \
      Ipp##IPPT f, DynMatrix<Ipp##IPPT> &dst) const{         \
    dst.setBounds(cols(),rows());                            \
    ippsMulC_##IPPT(data(),f, dst.data(),dim());             \
    return dst;                                              \
  }

DYN_MATRIX_MULT_BY_CONSTANT(32f)
DYN_MATRIX_MULT_BY_CONSTANT(64f)

#undef DYN_MATRIX_MULT_BY_CONSTANT
 
#define DYN_MATRIX_NORM_SPECIALZE(T,IPPT)                  \
  template<>                                               \
  inline T DynMatrix<T> ::norm(double l) const{            \
    if(l==1){                                              \
      T val;                                               \
      ippsNorm_L1_##IPPT(m_data,dim(),&val);               \
      return val;                                          \
    }else if(l==2){                                        \
      T val;                                               \
      ippsNorm_L2_##IPPT(m_data,dim(),&val);               \
      return val;                                          \
    }                                                      \
    double accu = 0;                                       \
    for(unsigned int i=0;i<dim();++i){                     \
      accu += ::pow(double(m_data[i]),l);                  \
    }                                                      \
    return ::pow(accu,1.0/l);                              \
  }
 
 DYN_MATRIX_NORM_SPECIALZE(float,32f)
 // DYN_MATRIX_NORM_SPECIALZE(double,64f)
 
#undef DYN_MATRIX_NORM_SPECIALZE

 /** \endcond */

#endif // HAVE_IPP

}

#endif