DynMatrix.h

Go to the documentation of this file.

/********************************************************************
**                Image Component Library (ICL)                    **
**                                                                 **
** Copyright (C) 2006-2013 CITEC, University of Bielefeld          **
**                         Neuroinformatics Group                  **
** Website: www.iclcv.org and                                      **
**          http://opensource.cit-ec.de/projects/icl               **
**                                                                 **
** File   : ICLMath/src/ICLMath/DynMatrix.h                        **
** Module : ICLMath                                                **
** Authors: Christof Elbrechter                                    **
**                                                                 **
**                                                                 **
** GNU LESSER GENERAL PUBLIC LICENSE                               **
** This file may be used under the terms of the GNU Lesser General **
** Public License version 3.0 as published by the                  **
**                                                                 **
** Free Software Foundation and appearing in the file LICENSE.LGPL **
** included in the packaging of this file.  Please review the      **
** following information to ensure the license requirements will   **
** be met: http://www.gnu.org/licenses/lgpl-3.0.txt                **
**                                                                 **
** The development of this software was supported by the           **
** Excellence Cluster EXC 277 Cognitive Interaction Technology.    **
** The Excellence Cluster EXC 277 is a grant of the Deutsche       **
** Forschungsgemeinschaft (DFG) in the context of the German       **
** Excellence Initiative.                                          **
**                                                                 **
********************************************************************/

#pragma once

#include <ICLUtils/Macros.h>
#include <ICLUtils/Exception.h>

#include <iterator>
#include <algorithm>
#include <numeric>
#include <functional>
#include <vector>
#include <cmath>

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

// Intel Math Kernel Library
#ifdef ICL_HAVE_MKL
#include "mkl_cblas.h"
#endif

namespace icl{
  namespace math{

    struct InvalidMatrixDimensionException :public utils::ICLException{
      InvalidMatrixDimensionException(const std::string &msg):utils::ICLException(msg){}
    };

    struct IncompatibleMatrixDimensionException :public utils::ICLException{
      IncompatibleMatrixDimensionException(const std::string &msg):utils::ICLException(msg){}
    };

    struct InvalidIndexException : public utils::ICLException{
      InvalidIndexException(const std::string &msg):utils::ICLException(msg){}
    };

    struct SingularMatrixException : public utils::ICLException{
      SingularMatrixException(const std::string &msg):utils::ICLException(msg){}
    };


    template<class T>
    struct DynMatrix{

      class DynMatrixColumn;
      DynMatrix(const DynMatrixColumn &column);

      inline DynMatrix():m_rows(0),m_cols(0),m_data(0),m_ownData(true){}

      inline DynMatrix(unsigned int cols,unsigned int rows,const  T &initValue=0)  :
      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);
      }


      inline DynMatrix(unsigned int cols,unsigned int rows, T *data, bool deepCopy=true)  :
        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;
        }
      }

      inline DynMatrix(unsigned int cols,unsigned int rows,const T *data)  :
        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());
      }

      inline DynMatrix(const DynMatrix &other):
        m_rows(other.m_rows),m_cols(other.m_cols),m_data(dim() ? new T[dim()] : 0),m_ownData(true){
        std::copy(other.begin(),other.end(),begin());
      }


      inline DynMatrix(const std::string &filename):m_rows(0),m_cols(0),m_data(0),m_ownData(true){
        *this = loadCSV(filename);
      }


      static DynMatrix<T> loadCSV(const std::string &filename) ;


      void saveCSV(const std::string &filename) ;

      inline bool isNull() const { return !m_data; }

      inline ~DynMatrix(){
        if(m_data && m_ownData) delete [] m_data;
      }


      inline DynMatrix &operator=(const DynMatrix &other){
        if(!m_data && !other.m_ownData){
          m_data = other.m_data;
          m_ownData = false;
          m_rows = other.m_rows;
          m_cols = other.m_cols;
        }else{
          if(dim() != other.dim()){
            delete[] m_data;
            m_data = other.dim() ? new T[other.dim()] : 0;
          }
          m_cols = other.m_cols;
          m_rows = other.m_rows;

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

      inline void setBounds(unsigned int cols, unsigned int rows, bool holdContent=false, const T &initializer=0) {
        if((int)cols == m_cols && (int)rows==m_rows) return;
        if(cols*rows == 0) 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;
        if(m_data && m_ownData) delete [] m_data;
        m_data = M.begin();
        m_ownData = true;
        M.set_data(0);
      }

      inline bool isSimilar(const DynMatrix &other, T tollerance=T(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>0?diff:-diff) > tollerance) return false;
        }
        return true;
      }

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

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


      inline DynMatrix operator*(T f) const{
        DynMatrix dst(cols(),rows());
        return mult(f,dst);
      }

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

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

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

      inline DynMatrix &operator/=(T f){
        return this->operator*=(1/f);
      }

      inline DynMatrix &mult(const DynMatrix &m, DynMatrix &dst) const {
        if( cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be rows(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;
      }

      inline DynMatrix &elementwise_mult(const DynMatrix &m, DynMatrix &dst) const {
        if((m.cols() != cols()) || (m.rows() != rows())) throw IncompatibleMatrixDimensionException("A.*B dimension mismatch");
        dst.setBounds(cols(),rows());
        for(unsigned int i=0;i<dim();++i){
        dst[i] = m_data[i] * m[i];
        }
        return dst;
      }

      inline DynMatrix elementwise_mult(const DynMatrix &m) const {
        DynMatrix dst(cols(),rows());
        return elementwise_mult(m,dst);
      }

      inline DynMatrix &elementwise_div(const DynMatrix &m, DynMatrix &dst) const {
        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;
      }

      inline DynMatrix elementwise_div(const DynMatrix &m) const {
        DynMatrix dst(cols(),rows());
        return elementwise_div(m,dst);
      }




      inline DynMatrix operator*(const DynMatrix &m) const {
        DynMatrix d(m.cols(),rows());
        return mult(m,d);
      }

      inline DynMatrix &operator*=(const DynMatrix &m) {
        return *this=((*this)*m);
      }

      inline DynMatrix operator/(const DynMatrix &m) const {
        return this->operator*(m.inv());
      }

      inline DynMatrix &operator/=(const DynMatrix &m) const {
        return *this = this->operator*(m.inv());
      }

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

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

      inline DynMatrix &operator+=(const T &t){
        std::transform(begin(),end(),begin(),std::bind2nd(std::plus<T>(),t));
        return *this;
      }

      inline DynMatrix &operator-=(const T &t){
        std::transform(begin(),end(),begin(),std::bind2nd(std::minus<T>(),t));
        return *this;
      }

      inline DynMatrix operator+(const DynMatrix &m) const {
        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;
      }

      inline DynMatrix operator-(const DynMatrix &m) const {
        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;
      }

      inline DynMatrix &operator+=(const DynMatrix &m) {
        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;
      }

      inline DynMatrix &operator-=(const DynMatrix &m) {
        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;
      }

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

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

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

      inline const T &at(unsigned int col,unsigned int row) const {
        return const_cast<DynMatrix*>(this)->at(col,row);
      }



      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];

      }


      inline const T &operator[](unsigned int idx) const {
        idx_check(idx);
        return m_data[idx];
      }

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

    private:
      static T diff_power_two(const T&a, const T&b){
        T d = a-b;
        return d*d;
      }
    public:
      inline T sqrDistanceTo(const DynMatrix &other) const {
        ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch"));
        return std::inner_product(begin(),end(),other.begin(),T(0), std::plus<T>(), diff_power_two);
      }

      inline T distanceTo(const DynMatrix &other) const {
        ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch"));
        return ::sqrt( distanceTo(other) );
      }


      typedef T* iterator;

      typedef const T* const_iterator;

      typedef T* row_iterator;

      typedef const T* const_row_iterator;

      unsigned int rows() const { return m_rows; }

      unsigned int cols() const { return m_cols; }

      T *data() { return m_data; }

      const T *data() const { return m_data; }

      unsigned int dim() const { return m_rows*m_cols; }

      int stride0() const { return sizeof(T) * dim(); }

      int stride1() const { return sizeof(T) * cols(); }

      int stride2() const { return sizeof(T); }

      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){}


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

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

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

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

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

        inline col_iterator &operator+=(difference_type n){
          p += n * stride;
          return *this;
        }

        inline const col_iterator &operator+=(difference_type n) const{
          p += n * stride;
          return *this;
        }


        inline col_iterator &operator-=(difference_type n){
          p -= n * stride;
          return *this;
        }

        inline const col_iterator &operator-=(difference_type n) const{
          p -= n * stride;
          return *this;
        }


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

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

        inline col_iterator operator-(difference_type n) {
          col_iterator tmp = *this;
          tmp-=n;
          return tmp;
        }

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


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

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

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

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

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

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

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

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

      typedef const col_iterator const_col_iterator;

      class ICLMath_API DynMatrixColumn{
        public:
  #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
        DynMatrix<T> *matrix;

        unsigned int column;

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

        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");
        }
        inline DynMatrixColumn(const DynMatrixColumn &c):
        matrix(c.matrix),column(c.column){}

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

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

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

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

        inline unsigned int dim() const { return matrix->rows(); }

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

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

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

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

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

        inline DynMatrixColumn &operator*=(const T&val){
          std::for_each(begin(),end(),std::bind2nd(std::multiplies<T>(),val));
          return *this;
        }
        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



      inline iterator begin() { return m_data; }

      inline iterator end() { return m_data+dim(); }

      inline const_iterator begin() const { return m_data; }

      inline const_iterator end() const { return m_data+dim(); }

      inline col_iterator col_begin(unsigned int col) {
        col_check(col);
        return col_iterator(m_data+col,cols());
      }

      inline col_iterator col_end(unsigned int col) {
        col_check(col);
        return col_iterator(m_data+col+dim(),cols());
      }

      inline const_col_iterator col_begin(unsigned int col) const {
        col_check(col);
        return col_iterator(m_data+col,cols());
      }

      inline const_col_iterator col_end(unsigned int col) const {
        col_check(col);
        return col_iterator(m_data+col+dim(),cols());
      }

      inline row_iterator row_begin(unsigned int row) {
        row_check(row);
        return m_data+row*cols();
      }

      inline row_iterator row_end(unsigned int row) {
        row_check(row);
        return m_data+(row+1)*cols();
      }

      inline const_row_iterator row_begin(unsigned int row) const {
        row_check(row);
        return m_data+row*cols();
      }

      inline const_row_iterator row_end(unsigned int row) const {
        row_check(row);
        return m_data+(row+1)*cols();
      }

      inline DynMatrix row(int row){
        row_check(row);
        return DynMatrix(m_cols,1,row_begin(row),false);
      }

      inline const DynMatrix row(int row) const{
        row_check(row);
        return DynMatrix(m_cols,1,const_cast<T*>(row_begin(row)),false);
      }

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

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

      void decompose_QR(DynMatrix &Q, DynMatrix &R) const
        ;

      void decompose_RQ(DynMatrix &R, DynMatrix &Q) const
        ;


      void decompose_LU(DynMatrix &L, DynMatrix &U, T zeroThreshold = T(1E-16)) const;

      DynMatrix solve_upper_triangular(const DynMatrix &b) const;

      DynMatrix solve_lower_triangular(const DynMatrix &b) const;


      DynMatrix solve(const DynMatrix &b, const std::string &method = "lu", T zeroThreshold = T(1E-16));


      DynMatrix inv() const ;


      void eigen(DynMatrix &eigenvectors, DynMatrix &eigenvalues) const ;


      void svd(DynMatrix &U, DynMatrix &S, DynMatrix &V) const ;


      DynMatrix pinv(bool useSVD = false, T zeroThreshold = T(1E-16)) const
        ;


      DynMatrix big_matrix_pinv(T zeroThreshold = T(1E-16)) const
        ;

  #ifdef ICL_HAVE_MKL
      typedef void(*GESDD)(const char*, const int*, const int*, T*, const int*, T*, T*, const int*, T*, const int*, T*, const int*, int*, int*);
      typedef void(*CBLAS_GEMM)(CBLAS_ORDER,CBLAS_TRANSPOSE,CBLAS_TRANSPOSE,int,int,int,T,const T*,int,const T*,int,T,T*,int);
      DynMatrix big_matrix_pinv(T zeroThreshold, GESDD gesdd, CBLAS_GEMM cblas_gemm) const
        ;
  #endif

      T det() const ;

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


      inline const DynMatrix<T> shallowTransposed() const{
        return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false);
      }

      const DynMatrix<T> shallowTransposed() {
        return DynMatrix<T>(m_rows,m_cols,const_cast<T*>(m_data),false);
      }


      inline void reshape(int newCols, int newRows) {
        if((cols() * rows()) != (newCols * newRows)){
          throw InvalidMatrixDimensionException("DynMatrix<T>::reshape: source dimension and destination dimension differs!");
        }
        m_cols = newCols;
        m_rows = newRows;
      }


      T element_wise_inner_product(const DynMatrix<T> &other) const {
        return std::inner_product(begin(),end(),other.begin(),T(0));
      }



      DynMatrix<T> dot(const DynMatrix<T> &M) const {
        return this->transp() * M;
      }


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

      T trace() const{
        ICLASSERT_RETURN_VAL(cols()==rows(),0);
        double accu = 0;
        for(unsigned int i=0;i<dim();i+=cols()+1){
          accu += m_data[i];
        }
        return accu;
      }

      static DynMatrix<T> cross(const DynMatrix<T> &x, const DynMatrix<T> &y){
        if(x.cols()==1 && y.cols()==1 && x.rows()==3 && y.rows()==3){
            DynMatrix<T> r(1,x.rows());
            r(0,0) = x(0,1)*y(0,2)-x(0,2)*y(0,1);
            r(0,1) = x(0,2)*y(0,0)-x(0,0)*y(0,2);
            r(0,2) = x(0,0)*y(0,1)-x(0,1)*y(0,0);
            return r;
        }else{
            ICLASSERT_RETURN_VAL(x.rows() == 3 && y.rows() == 3,DynMatrix<T>());
            return DynMatrix<T>();
        }
      }

      T cond(const double p=2) const {
        if(cols() == 3 && rows() == 3){
            DynMatrix<T> M_inv = (*this).inv();
            return (*this).norm(p) * M_inv.norm(p);
        } else {
            DynMatrix<T> U,S,V;
            (*this).svd(U,S,V);
            if(S[S.rows()-1]){
                return S[0]/S[S.rows()-1];
            } else {
                return S[0];
            }
        }
      }

      inline T *set_data(T *newData){
        T *old_data = m_data;
        m_data = newData;
        return old_data;
      }

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

    template<class T>
    DynMatrix<T>::DynMatrix(const typename DynMatrix<T>::DynMatrixColumn &column) :
    m_rows(column.dim()),m_cols(1),m_data(new T[column.dim()]),m_ownData(true){
      std::copy(column.begin(),column.end(),begin());
    }
    template<class T> ICLMath_IMP
    std::ostream &operator<<(std::ostream &s, const DynMatrix<T> &m);

    template<class T> ICLMath_IMP
    std::istream &operator>>(std::istream &s, DynMatrix<T> &m);


#ifdef ICL_HAVE_IPP

    template<>
    inline float DynMatrix<float>::sqrDistanceTo(const DynMatrix<float> &other) const {
      ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch"));
      float norm = 0 ;
      ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm);
      return norm*norm;
    }

    template<>
    inline double DynMatrix<double>::sqrDistanceTo(const DynMatrix<double> &other) const {
      ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::sqrDistanceTo: dimension missmatch"));
      double norm = 0 ;
      ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm);
      return norm*norm;
    }

    template<>
    inline float DynMatrix<float>::distanceTo(const DynMatrix<float> &other) const {
      ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch"));
      float norm = 0 ;
      ippsNormDiff_L2_32f(begin(), other.begin(), dim(), &norm);
      return norm;
    }

    template<>
    inline double DynMatrix<double>::distanceTo(const DynMatrix<double> &other) const {
      ICLASSERT_THROW(dim() == other.dim(), InvalidMatrixDimensionException("DynMatrix::distanceTo: dimension missmatch"));
      double norm = 0 ;
      ippsNormDiff_L2_64f(begin(), other.begin(), dim(), &norm);
      return norm;
    }

#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 \
    {                       \
      if(cols() != m.rows() ) throw IncompatibleMatrixDimensionException("A*B : cols(A) must be row(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 \
    {                       \
      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

#endif // ICL_HAVE_IPP


    template<class T>
    inline DynMatrix<T> operator,(const DynMatrix<T> &left, const DynMatrix<T> &right){
      int w = left.cols() + right.cols();
      int h = iclMax(left.rows(),right.rows());
      DynMatrix<T> result(w,h,float(0));
      for(unsigned int y=0;y<left.rows();++y){
        std::copy(left.row_begin(y), left.row_end(y), result.row_begin(y));
      }
      for(unsigned int y=0;y<right.rows();++y){
        std::copy(right.row_begin(y), right.row_end(y), result.row_begin(y) + left.cols());
      }
      return result;
    }


    template<class T>
    inline DynMatrix<T> operator%(const DynMatrix<T> &top, const DynMatrix<T> &bottom){
      int w = iclMax(top.cols(),bottom.cols());
      int h = top.rows() + bottom.rows();
      DynMatrix<T> result(w,h,float(0));
      for(unsigned int y=0;y<top.rows();++y){
        std::copy(top.row_begin(y), top.row_end(y), result.row_begin(y));
      }
      for(unsigned int y=0;y<bottom.rows();++y){
        std::copy(bottom.row_begin(y), bottom.row_end(y), result.row_begin(y+top.rows()));
      }
      return result;
    }
  } // namespace math
}