src/math/DynamicRectMatrix.hpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */
 
/**
 * @file DynamicRectMatrix.hpp
 * @author Teng Lin
 * @date 10/11/2004
 * @version 1.0
 */

#ifndef MATH_DYNAMICRECTMATRIX_HPP
#define MATH_DYNAMICRECTMATRIX_HPP
#include <math.h>
#include <cmath>
#include "math/DynamicVector.hpp"

namespace OpenMD {

  /**
   * @class DynamicRectMatrix DynamicRectMatrix.hpp "math/DynamicRectMatrix.hpp"
   * @brief rectangular matrix class
   */
  template<typename Real> 
  class DynamicRectMatrix {
  public:
    typedef Real ElemType;
    typedef Real* ElemPoinerType;
    typedef DynamicRectMatrix<Real> SelfType;
            
    /** default constructor */
    DynamicRectMatrix(){
      nrow_ = 0;
      ncol_ = 0;
      data_ = NULL;
    }

    DynamicRectMatrix(int nrow, int ncol) {
      allocate(nrow, ncol);

      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = 0.0;
    }

    /** Constructs and initializes every element of this matrix to a scalar */ 
    DynamicRectMatrix(int nrow, int ncol, Real s) {
      allocate(nrow, ncol);
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = s;
    }

    DynamicRectMatrix(int nrow, int ncol, Real* array) {
      allocate(nrow, ncol);
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = array[i * nrow_ + j];
    }

    /** copy constructor */
    DynamicRectMatrix(const SelfType& m) {
      allocate(m.getNRow(), m.getNCol());
      
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
            this->data_[i][j] = m.data_[i][j];
    }
            
    /** destructor*/
    ~DynamicRectMatrix() { deallocate();}

    /** copy assignment operator */
    DynamicRectMatrix<Real> operator =(const DynamicRectMatrix<Real> m) {
      if (this == &m)
          return *this;
      if (nrow_ != m.getNRow() || ncol_ != m.getNCol()) {
        deallocate();
        allocate(m.getNRow(), m.getNCol());
      }
      
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = m.data_[i][j];
      return *this;
    }
            
    /**
     * Return the reference of a single element of this matrix.
     * @return the reference of a single element of this matrix 
     * @param i row index
     * @param j Column index
     */
    Real& operator()(unsigned int i, unsigned int j) {
      return this->data_[i][j];
    }

    /**
     * Return the value of a single element of this matrix.
     * @return the value of a single element of this matrix 
     * @param i row index
     * @param j Column index
     */        
    Real operator()(unsigned int i, unsigned int j) const  {
                
      return this->data_[i][j];  
    }

    /** 
     * Copy the internal data to an array
     * @param array the pointer of destination array
     */
    void getArray(Real* array) {
      for (unsigned int i = 0; i < nrow_; i++) {
        for (unsigned int j = 0; j < ncol_; j++) {
          array[i * nrow_ + j] = this->data_[i][j];
        }
      }
    }

    /**
     * Returns a row of  this matrix as a vector.
     * @return a row of  this matrix as a vector 
     * @param row the row index
     */                
    DynamicVector<Real> getRow(unsigned int row) {
      DynamicVector<Real> v;

      for (unsigned int i = 0; i < ncol_; i++)
        v[i] = this->data_[row][i];

      return v;
    }

    /**
     * Sets a row of  this matrix
     * @param row the row index
     * @param v the vector to be set
     */                
    void setRow(unsigned int row, const DynamicVector<Real>& v) {
      assert(v.size() == nrow_);
      for (unsigned int i = 0; i < ncol_; i++)
        this->data_[row][i] = v[i];
    }

    /**
     * Returns a column of  this matrix as a vector.
     * @return a column of  this matrix as a vector 
     * @param col the column index
     */                
    DynamicVector<Real> getColumn(unsigned int col) {
      DynamicVector<Real> v(ncol_);

      for (unsigned int j = 0; j < nrow_; j++)
        v[j] = this->data_[j][col];

      return v;
    }

    /**
     * Sets a column of  this matrix
     * @param col the column index
     * @param v the vector to be set
     */                
    void setColumn(unsigned int col, const DynamicVector<Real>& v){

      for (unsigned int j = 0; j < nrow_; j++)
        this->data_[j][col] = v[j];
    }         

    /**
     * swap two rows of this matrix
     * @param i the first row
     * @param j the second row
     */
    void swapRow(unsigned int i, unsigned int j){
      assert(i < nrow_ && j < nrow_);

      for (unsigned int k = 0; k < ncol_; k++)
        std::swap(this->data_[i][k], this->data_[j][k]);
    }

    /**
     * swap two Columns of this matrix
     * @param i the first Column
     * @param j the second Column
     */
    void swapColumn(unsigned int i, unsigned int j){
      assert(i < ncol_ && j < ncol_);
                    
      for (unsigned int k = 0; k < nrow_; k++)
        std::swap(this->data_[k][i], this->data_[k][j]);
    }

    /**
     * Tests if this matrix is identical to matrix m
     * @return true if this matrix is equal to the matrix m, return false otherwise
     * @m matrix to be compared
     *
     * @todo replace operator == by template function equal
     */
    bool operator ==(const DynamicRectMatrix<Real> m) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          if (!equal(this->data_[i][j], m.data_[i][j]))
            return false;

      return true;
    }

    /**
     * Tests if this matrix is not equal to matrix m
     * @return true if this matrix is not equal to the matrix m, return false otherwise
     * @m matrix to be compared
     */
    bool operator !=(const DynamicRectMatrix<Real> m) {
      return !(*this == m);
    }

    /** Negates the value of this matrix in place. */           
    inline void negate() {
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = -this->data_[i][j];
    }
            
    /**
     * Sets the value of this matrix to the negation of matrix m.
     * @param m the source matrix
     */
    inline void negate(const DynamicRectMatrix<Real> m) {
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)
          this->data_[i][j] = -m.data_[i][j];        
    }
            
    /**
     * Sets the value of this matrix to the sum of itself and m (*this += m).
     * @param m the other matrix
     */
    inline void add( const DynamicRectMatrix<Real> m ) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)        
          this->data_[i][j] += m.data_[i][j];
    }
            
    /**
     * Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
     * @param m1 the first matrix
     * @param m2 the second matrix
     */
    inline void add( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2 ) {
      assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)        
          this->data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
    }
            
    /**
     * Sets the value of this matrix to the difference  of itself and m (*this -= m).
     * @param m the other matrix
     */
    inline void sub( const DynamicRectMatrix<Real> m ) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)        
          this->data_[i][j] -= m.data_[i][j];
    }
            
    /**
     * Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
     * @param m1 the first matrix
     * @param m2 the second matrix
     */
    inline void sub( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2){
      assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)        
          this->data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
    }

    /**
     * Sets the value of this matrix to the scalar multiplication of itself (*this *= s).
     * @param s the scalar value
     */
    inline void mul( Real s ) {
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)  
          this->data_[i][j] *= s;
    }

    /**
     * Sets the value of this matrix to the scalar multiplication of matrix m  (*this = s * m).
     * @param s the scalar value
     * @param m the matrix
     */
    inline void mul( Real s, const DynamicRectMatrix<Real> m ) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());    
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)  
          this->data_[i][j] = s * m.data_[i][j];
    }

    /**
     * Sets the value of this matrix to the scalar division of itself  (*this /= s ).
     * @param s the scalar value
     */             
    inline void div( Real s) {
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)  
          this->data_[i][j] /= s;
    }

    /**
     * Sets the value of this matrix to the scalar division of matrix m  (*this = m /s).
     * @param s the scalar value
     * @param m the matrix
     */
    inline void div( Real s, const DynamicRectMatrix<Real> m ) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)  
          this->data_[i][j] = m.data_[i][j] / s;
    }

    /**
     *  Multiples a scalar into every element of this matrix.
     * @param s the scalar value
     */
    DynamicRectMatrix<Real> operator *=(const Real s) {
      this->mul(s);
      return *this;
    }

    /**
     *  Divides every element of this matrix by a scalar.
     * @param s the scalar value
     */
    DynamicRectMatrix<Real> operator /=(const Real s) {
      this->div(s);
      return *this;
    }

    /**
     * Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
     * @param m the other matrix
     */
    DynamicRectMatrix<Real> operator += (const DynamicRectMatrix<Real> m) {
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      add(m);
      return *this;
    }

    /**
     * Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m) 
     * @param m the other matrix
     */
    DynamicRectMatrix<Real> operator -= (const DynamicRectMatrix<Real> m){
      assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
      sub(m);
      return *this;
    }

    /** Return the transpose of this matrix */
    DynamicRectMatrix<Real> transpose() const{
      DynamicRectMatrix<Real> result(ncol_,nrow_);
                
      for (unsigned int i = 0; i < nrow_; i++)
        for (unsigned int j = 0; j < ncol_; j++)              
          result(j, i) = this->data_[i][j];

      return result;
    }

    unsigned int getNRow() const {return nrow_;}
    unsigned int getNCol() const {return ncol_;}

    template<class MatrixType>
    void setSubMatrix(unsigned int beginRow, unsigned int beginCol, const MatrixType& m) {
        assert(beginRow + m.getNRow() -1 <= nrow_);
        assert(beginCol + m.getNCol() -1 <= ncol_);

        for (unsigned int i = 0; i < m.getNRow(); ++i)
            for (unsigned int j = 0; j < m.getNCol(); ++j)
                this->data_[beginRow+i][beginCol+j] = m(i, j);
    }

    template<class MatrixType>
    void getSubMatrix(unsigned int beginRow, unsigned int beginCol, MatrixType& m) {
        assert(beginRow + m.getNRow() -1 <= nrow_);
        assert(beginCol + m.getNCol() - 1 <= ncol_);

        for (unsigned int i = 0; i < m.getNRow(); ++i)
            for (unsigned int j = 0; j < m.getNCol(); ++j)
                m(i, j) = this->data_[beginRow+i][beginCol+j];
    }
    
  protected:
    Real** data_;
    unsigned int nrow_;
    unsigned int ncol_;
  private:
    void allocate(int nrow, int ncol) {
        nrow_ = nrow;
        ncol_ = ncol;
        data_ = new Real*[nrow_];
        for (int i = 0; i < nrow_; ++i)
            data_[i] = new Real[ncol_];
    }

    void deallocate() {
      for (int i = 0; i < nrow_; ++i)
        delete data_[i];
      delete []data_;

      nrow_ = 0;
      ncol_ = 0;
      data_ = NULL;
    }
    
  };

  /** Negate the value of every element of this matrix. */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator -(const DynamicRectMatrix<Real> m) {
    DynamicRectMatrix<Real> result(m);

    result.negate();

    return result;
  }
    
  /**
   * Return the sum of two matrixes  (m1 + m2). 
   * @return the sum of two matrixes
   * @param m1 the first matrix
   * @param m2 the second matrix
   */ 
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator + (const DynamicRectMatrix<Real> m1,const DynamicRectMatrix<Real> m2) {
    
    DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());

    result.add(m1, m2);

    return result;
  }
    
  /**
   * Return the difference of two matrixes  (m1 - m2). 
   * @return the sum of two matrixes
   * @param m1 the first matrix
   * @param m2 the second matrix
   */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator - (const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2) {
    DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());

    result.sub(m1, m2);

    return result;
  }

  /**
   * Return the multiplication of scalra and  matrix  (m * s). 
   * @return the multiplication of a scalra and  a matrix 
   * @param m the matrix
   * @param s the scalar
   */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real> m, Real s) {
    DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());

    result.mul(s, m);

    return result;
  }

  /**
   * Return the multiplication of a scalra and  a matrix  (s * m). 
   * @return the multiplication of a scalra and  a matrix 
   * @param s the scalar
   * @param m the matrix
   */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator *(Real s, const DynamicRectMatrix<Real> m) {
    DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());

    result.mul(s, m);

    return result;
  }
    
  /**
   * Return the multiplication of two matrixes  (m1 * m2). 
   * @return the multiplication of two matrixes
   * @param m1 the first matrix
   * @param m2 the second matrix
   */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real>& m1, const DynamicRectMatrix<Real>& m2) {
    assert(m1.getNCol() == m2.getNRow());
    unsigned int sameDim = m1.getNCol();
    int nrow = m1.getNRow();
    int ncol = m2.getNCol();
    DynamicRectMatrix<Real> result(nrow, ncol );
    for (unsigned int i = 0; i < nrow; i++)
      for (unsigned int j = 0; j < ncol; j++)
        for (unsigned int k = 0; k < sameDim; k++)
          result(i, j)  += m1(i, k) * m2(k, j);                

    return result;
  }
    
  /**
   * Return the multiplication of  a matrix and a vector  (m * v). 
   * @return the multiplication of a matrix and a vector
   * @param m the matrix
   * @param v the vector
   */
  template<typename Real>
  inline DynamicVector<Real> operator *(const DynamicRectMatrix<Real> m, const DynamicVector<Real>& v) {
    int nrow = m.getNRow();
    int ncol = m.getNCol();
    assert(ncol = v.size());
    DynamicVector<Real> result(nrow);
    
    for (unsigned int i = 0; i < nrow ; i++)
      for (unsigned int j = 0; j < ncol ; j++)            
        result[i] += m(i, j) * v[j];
            
    return result;                                                                 
  }

  /**
   * Return the scalar division of matrix   (m / s). 
   * @return the scalar division of matrix  
   * @param m the matrix
   * @param s the scalar
   */
  template<typename Real> 
  inline DynamicRectMatrix<Real> operator /(const DynamicRectMatrix<Real> m, Real s) {
    DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());

    result.div(s, m);

    return result;
  }    

  /**
   * Write to an output stream
   */
  template<typename Real>
  std::ostream &operator<< ( std::ostream& o, const DynamicRectMatrix<Real> m) {
    for (unsigned int i = 0; i < m.getNRow() ; i++) {
      o << "(";
      for (unsigned int j = 0; j < m.getNCol() ; j++) {
        o << m(i, j);
        if (j != m.getNCol() -1)
          o << "\t";
      }
      o << ")" << std::endl;
    }
    return o;        
  }    
}
#endif //MATH_RECTMATRIX_HPP

Revision:	1390
Committed:	Wed Nov 25 20:02:06 2009 UTC (15 years, 5 months ago) by gezelter
File size:	18857 byte(s)
Log Message:	Almost all of the changes necessary to create OpenMD out of our old project (OOPSE-4)
#	User	Rev	Content
1	tim	891	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	tim	891	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	tim	891	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			* [4] Vardeman & Gezelter, in progress (2009).
40	tim	891	*/
41
42			/**
43			* @file DynamicRectMatrix.hpp
44			* @author Teng Lin
45			* @date 10/11/2004
46			* @version 1.0
47			*/
48
49			#ifndef MATH_DYNAMICRECTMATRIX_HPP
50			#define MATH_DYNAMICRECTMATRIX_HPP
51			#include <math.h>
52			#include <cmath>
53			#include "math/DynamicVector.hpp"
54
55	gezelter	1390	namespace OpenMD {
56	tim	891
57			/**
58			* @class DynamicRectMatrix DynamicRectMatrix.hpp "math/DynamicRectMatrix.hpp"
59			* @brief rectangular matrix class
60			*/
61			template<typename Real>
62			class DynamicRectMatrix {
63			public:
64			typedef Real ElemType;
65			typedef Real* ElemPoinerType;
66			typedef DynamicRectMatrix<Real> SelfType;
67
68			/** default constructor */
69	cli2	1360	DynamicRectMatrix(){
70			nrow_ = 0;
71			ncol_ = 0;
72			data_ = NULL;
73			}
74
75	tim	891	DynamicRectMatrix(int nrow, int ncol) {
76			allocate(nrow, ncol);
77
78			for (unsigned int i = 0; i < nrow_; i++)
79			for (unsigned int j = 0; j < ncol_; j++)
80			this->data_[i][j] = 0.0;
81			}
82
83			/** Constructs and initializes every element of this matrix to a scalar */
84			DynamicRectMatrix(int nrow, int ncol, Real s) {
85			allocate(nrow, ncol);
86			for (unsigned int i = 0; i < nrow_; i++)
87			for (unsigned int j = 0; j < ncol_; j++)
88			this->data_[i][j] = s;
89			}
90
91			DynamicRectMatrix(int nrow, int ncol, Real* array) {
92			allocate(nrow, ncol);
93			for (unsigned int i = 0; i < nrow_; i++)
94			for (unsigned int j = 0; j < ncol_; j++)
95			this->data_[i][j] = array[i * nrow_ + j];
96			}
97
98			/** copy constructor */
99			DynamicRectMatrix(const SelfType& m) {
100			allocate(m.getNRow(), m.getNCol());
101
102			for (unsigned int i = 0; i < nrow_; i++)
103			for (unsigned int j = 0; j < ncol_; j++)
104			this->data_[i][j] = m.data_[i][j];
105			}
106
107			/** destructor*/
108			~DynamicRectMatrix() { deallocate();}
109
110			/** copy assignment operator */
111			DynamicRectMatrix<Real> operator =(const DynamicRectMatrix<Real> m) {
112			if (this == &m)
113			return *this;
114			if (nrow_ != m.getNRow() \|\| ncol_ != m.getNCol()) {
115			deallocate();
116			allocate(m.getNRow(), m.getNCol());
117			}
118
119			for (unsigned int i = 0; i < nrow_; i++)
120			for (unsigned int j = 0; j < ncol_; j++)
121			this->data_[i][j] = m.data_[i][j];
122			return *this;
123			}
124
125			/**
126			* Return the reference of a single element of this matrix.
127			* @return the reference of a single element of this matrix
128			* @param i row index
129			* @param j Column index
130			*/
131			Real& operator()(unsigned int i, unsigned int j) {
132			return this->data_[i][j];
133			}
134
135			/**
136			* Return the value of a single element of this matrix.
137			* @return the value of a single element of this matrix
138			* @param i row index
139			* @param j Column index
140			*/
141			Real operator()(unsigned int i, unsigned int j) const {
142
143			return this->data_[i][j];
144			}
145
146			/**
147			* Copy the internal data to an array
148			* @param array the pointer of destination array
149			*/
150			void getArray(Real* array) {
151			for (unsigned int i = 0; i < nrow_; i++) {
152			for (unsigned int j = 0; j < ncol_; j++) {
153			array[i * nrow_ + j] = this->data_[i][j];
154			}
155			}
156			}
157
158			/**
159			* Returns a row of this matrix as a vector.
160			* @return a row of this matrix as a vector
161			* @param row the row index
162			*/
163			DynamicVector<Real> getRow(unsigned int row) {
164			DynamicVector<Real> v;
165
166			for (unsigned int i = 0; i < ncol_; i++)
167			v[i] = this->data_[row][i];
168
169			return v;
170			}
171
172			/**
173			* Sets a row of this matrix
174			* @param row the row index
175			* @param v the vector to be set
176			*/
177			void setRow(unsigned int row, const DynamicVector<Real>& v) {
178			assert(v.size() == nrow_);
179			for (unsigned int i = 0; i < ncol_; i++)
180			this->data_[row][i] = v[i];
181			}
182
183			/**
184			* Returns a column of this matrix as a vector.
185			* @return a column of this matrix as a vector
186			* @param col the column index
187			*/
188			DynamicVector<Real> getColumn(unsigned int col) {
189			DynamicVector<Real> v(ncol_);
190
191			for (unsigned int j = 0; j < nrow_; j++)
192			v[j] = this->data_[j][col];
193
194			return v;
195			}
196
197			/**
198			* Sets a column of this matrix
199			* @param col the column index
200			* @param v the vector to be set
201			*/
202			void setColumn(unsigned int col, const DynamicVector<Real>& v){
203
204			for (unsigned int j = 0; j < nrow_; j++)
205			this->data_[j][col] = v[j];
206			}
207
208			/**
209			* swap two rows of this matrix
210			* @param i the first row
211			* @param j the second row
212			*/
213			void swapRow(unsigned int i, unsigned int j){
214			assert(i < nrow_ && j < nrow_);
215
216			for (unsigned int k = 0; k < ncol_; k++)
217			std::swap(this->data_[i][k], this->data_[j][k]);
218			}
219
220			/**
221			* swap two Columns of this matrix
222			* @param i the first Column
223			* @param j the second Column
224			*/
225			void swapColumn(unsigned int i, unsigned int j){
226			assert(i < ncol_ && j < ncol_);
227
228			for (unsigned int k = 0; k < nrow_; k++)
229			std::swap(this->data_[k][i], this->data_[k][j]);
230			}
231
232			/**
233			* Tests if this matrix is identical to matrix m
234			* @return true if this matrix is equal to the matrix m, return false otherwise
235			* @m matrix to be compared
236			*
237			* @todo replace operator == by template function equal
238			*/
239			bool operator ==(const DynamicRectMatrix<Real> m) {
240			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
241			for (unsigned int i = 0; i < nrow_; i++)
242			for (unsigned int j = 0; j < ncol_; j++)
243			if (!equal(this->data_[i][j], m.data_[i][j]))
244			return false;
245
246			return true;
247			}
248
249			/**
250			* Tests if this matrix is not equal to matrix m
251			* @return true if this matrix is not equal to the matrix m, return false otherwise
252			* @m matrix to be compared
253			*/
254			bool operator !=(const DynamicRectMatrix<Real> m) {
255			return !(*this == m);
256			}
257
258			/** Negates the value of this matrix in place. */
259			inline void negate() {
260			for (unsigned int i = 0; i < nrow_; i++)
261			for (unsigned int j = 0; j < ncol_; j++)
262			this->data_[i][j] = -this->data_[i][j];
263			}
264
265			/**
266			* Sets the value of this matrix to the negation of matrix m.
267			* @param m the source matrix
268			*/
269			inline void negate(const DynamicRectMatrix<Real> m) {
270			for (unsigned int i = 0; i < nrow_; i++)
271			for (unsigned int j = 0; j < ncol_; j++)
272			this->data_[i][j] = -m.data_[i][j];
273			}
274
275			/**
276			* Sets the value of this matrix to the sum of itself and m (*this += m).
277			* @param m the other matrix
278			*/
279			inline void add( const DynamicRectMatrix<Real> m ) {
280			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
281			for (unsigned int i = 0; i < nrow_; i++)
282			for (unsigned int j = 0; j < ncol_; j++)
283			this->data_[i][j] += m.data_[i][j];
284			}
285
286			/**
287			* Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
288			* @param m1 the first matrix
289			* @param m2 the second matrix
290			*/
291			inline void add( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2 ) {
292			assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
293			for (unsigned int i = 0; i < nrow_; i++)
294			for (unsigned int j = 0; j < ncol_; j++)
295			this->data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
296			}
297
298			/**
299			* Sets the value of this matrix to the difference of itself and m (*this -= m).
300			* @param m the other matrix
301			*/
302			inline void sub( const DynamicRectMatrix<Real> m ) {
303			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
304			for (unsigned int i = 0; i < nrow_; i++)
305			for (unsigned int j = 0; j < ncol_; j++)
306			this->data_[i][j] -= m.data_[i][j];
307			}
308
309			/**
310			* Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
311			* @param m1 the first matrix
312			* @param m2 the second matrix
313			*/
314			inline void sub( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2){
315			assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
316			for (unsigned int i = 0; i < nrow_; i++)
317			for (unsigned int j = 0; j < ncol_; j++)
318			this->data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
319			}
320
321			/**
322			* Sets the value of this matrix to the scalar multiplication of itself (this = s).
323			* @param s the scalar value
324			*/
325			inline void mul( Real s ) {
326			for (unsigned int i = 0; i < nrow_; i++)
327			for (unsigned int j = 0; j < ncol_; j++)
328			this->data_[i][j] *= s;
329			}
330
331			/**
332			* Sets the value of this matrix to the scalar multiplication of matrix m (this = s m).
333			* @param s the scalar value
334			* @param m the matrix
335			*/
336			inline void mul( Real s, const DynamicRectMatrix<Real> m ) {
337			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
338			for (unsigned int i = 0; i < nrow_; i++)
339			for (unsigned int j = 0; j < ncol_; j++)
340			this->data_[i][j] = s * m.data_[i][j];
341			}
342
343			/**
344			* Sets the value of this matrix to the scalar division of itself (*this /= s ).
345			* @param s the scalar value
346			*/
347			inline void div( Real s) {
348			for (unsigned int i = 0; i < nrow_; i++)
349			for (unsigned int j = 0; j < ncol_; j++)
350			this->data_[i][j] /= s;
351			}
352
353			/**
354			* Sets the value of this matrix to the scalar division of matrix m (*this = m /s).
355			* @param s the scalar value
356			* @param m the matrix
357			*/
358			inline void div( Real s, const DynamicRectMatrix<Real> m ) {
359			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
360			for (unsigned int i = 0; i < nrow_; i++)
361			for (unsigned int j = 0; j < ncol_; j++)
362			this->data_[i][j] = m.data_[i][j] / s;
363			}
364
365			/**
366			* Multiples a scalar into every element of this matrix.
367			* @param s the scalar value
368			*/
369			DynamicRectMatrix<Real> operator *=(const Real s) {
370			this->mul(s);
371			return *this;
372			}
373
374			/**
375			* Divides every element of this matrix by a scalar.
376			* @param s the scalar value
377			*/
378			DynamicRectMatrix<Real> operator /=(const Real s) {
379			this->div(s);
380			return *this;
381			}
382
383			/**
384			* Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
385			* @param m the other matrix
386			*/
387			DynamicRectMatrix<Real> operator += (const DynamicRectMatrix<Real> m) {
388			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
389			add(m);
390			return *this;
391			}
392
393			/**
394			* Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m)
395			* @param m the other matrix
396			*/
397			DynamicRectMatrix<Real> operator -= (const DynamicRectMatrix<Real> m){
398			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
399			sub(m);
400			return *this;
401			}
402
403			/** Return the transpose of this matrix */
404			DynamicRectMatrix<Real> transpose() const{
405			DynamicRectMatrix<Real> result(ncol_,nrow_);
406
407			for (unsigned int i = 0; i < nrow_; i++)
408			for (unsigned int j = 0; j < ncol_; j++)
409			result(j, i) = this->data_[i][j];
410
411			return result;
412			}
413
414			unsigned int getNRow() const {return nrow_;}
415			unsigned int getNCol() const {return ncol_;}
416
417			template<class MatrixType>
418			void setSubMatrix(unsigned int beginRow, unsigned int beginCol, const MatrixType& m) {
419			assert(beginRow + m.getNRow() -1 <= nrow_);
420			assert(beginCol + m.getNCol() -1 <= ncol_);
421
422			for (unsigned int i = 0; i < m.getNRow(); ++i)
423			for (unsigned int j = 0; j < m.getNCol(); ++j)
424			this->data_[beginRow+i][beginCol+j] = m(i, j);
425			}
426
427			template<class MatrixType>
428			void getSubMatrix(unsigned int beginRow, unsigned int beginCol, MatrixType& m) {
429			assert(beginRow + m.getNRow() -1 <= nrow_);
430			assert(beginCol + m.getNCol() - 1 <= ncol_);
431
432			for (unsigned int i = 0; i < m.getNRow(); ++i)
433			for (unsigned int j = 0; j < m.getNCol(); ++j)
434			m(i, j) = this->data_[beginRow+i][beginCol+j];
435			}
436
437			protected:
438			Real** data_;
439			unsigned int nrow_;
440			unsigned int ncol_;
441			private:
442			void allocate(int nrow, int ncol) {
443			nrow_ = nrow;
444			ncol_ = ncol;
445			data_ = new Real*[nrow_];
446			for (int i = 0; i < nrow_; ++i)
447			data_[i] = new Real[ncol_];
448			}
449
450			void deallocate() {
451			for (int i = 0; i < nrow_; ++i)
452			delete data_[i];
453			delete []data_;
454
455			nrow_ = 0;
456			ncol_ = 0;
457			data_ = NULL;
458			}
459
460			};
461
462			/** Negate the value of every element of this matrix. */
463			template<typename Real>
464			inline DynamicRectMatrix<Real> operator -(const DynamicRectMatrix<Real> m) {
465			DynamicRectMatrix<Real> result(m);
466
467			result.negate();
468
469			return result;
470			}
471
472			/**
473			* Return the sum of two matrixes (m1 + m2).
474			* @return the sum of two matrixes
475			* @param m1 the first matrix
476			* @param m2 the second matrix
477			*/
478			template<typename Real>
479			inline DynamicRectMatrix<Real> operator + (const DynamicRectMatrix<Real> m1,const DynamicRectMatrix<Real> m2) {
480
481			DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());
482
483			result.add(m1, m2);
484
485			return result;
486			}
487
488			/**
489			* Return the difference of two matrixes (m1 - m2).
490			* @return the sum of two matrixes
491			* @param m1 the first matrix
492			* @param m2 the second matrix
493			*/
494			template<typename Real>
495			inline DynamicRectMatrix<Real> operator - (const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2) {
496			DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());
497
498			result.sub(m1, m2);
499
500			return result;
501			}
502
503			/**
504			* Return the multiplication of scalra and matrix (m * s).
505			* @return the multiplication of a scalra and a matrix
506			* @param m the matrix
507			* @param s the scalar
508			*/
509			template<typename Real>
510			inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real> m, Real s) {
511			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
512
513			result.mul(s, m);
514
515			return result;
516			}
517
518			/**
519			* Return the multiplication of a scalra and a matrix (s * m).
520			* @return the multiplication of a scalra and a matrix
521			* @param s the scalar
522			* @param m the matrix
523			*/
524			template<typename Real>
525			inline DynamicRectMatrix<Real> operator *(Real s, const DynamicRectMatrix<Real> m) {
526			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
527
528			result.mul(s, m);
529
530			return result;
531			}
532
533			/**
534			* Return the multiplication of two matrixes (m1 * m2).
535			* @return the multiplication of two matrixes
536			* @param m1 the first matrix
537			* @param m2 the second matrix
538			*/
539			template<typename Real>
540			inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real>& m1, const DynamicRectMatrix<Real>& m2) {
541			assert(m1.getNCol() == m2.getNRow());
542			unsigned int sameDim = m1.getNCol();
543			int nrow = m1.getNRow();
544			int ncol = m2.getNCol();
545			DynamicRectMatrix<Real> result(nrow, ncol );
546			for (unsigned int i = 0; i < nrow; i++)
547			for (unsigned int j = 0; j < ncol; j++)
548			for (unsigned int k = 0; k < sameDim; k++)
549			result(i, j) += m1(i, k) * m2(k, j);
550
551			return result;
552			}
553
554			/**
555			* Return the multiplication of a matrix and a vector (m * v).
556			* @return the multiplication of a matrix and a vector
557			* @param m the matrix
558			* @param v the vector
559			*/
560			template<typename Real>
561			inline DynamicVector<Real> operator *(const DynamicRectMatrix<Real> m, const DynamicVector<Real>& v) {
562			int nrow = m.getNRow();
563			int ncol = m.getNCol();
564			assert(ncol = v.size());
565			DynamicVector<Real> result(nrow);
566
567			for (unsigned int i = 0; i < nrow ; i++)
568			for (unsigned int j = 0; j < ncol ; j++)
569			result[i] += m(i, j) * v[j];
570
571			return result;
572			}
573
574			/**
575			* Return the scalar division of matrix (m / s).
576			* @return the scalar division of matrix
577			* @param m the matrix
578			* @param s the scalar
579			*/
580			template<typename Real>
581			inline DynamicRectMatrix<Real> operator /(const DynamicRectMatrix<Real> m, Real s) {
582			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
583
584			result.div(s, m);
585
586			return result;
587			}
588
589			/**
590			* Write to an output stream
591			*/
592			template<typename Real>
593			std::ostream &operator<< ( std::ostream& o, const DynamicRectMatrix<Real> m) {
594			for (unsigned int i = 0; i < m.getNRow() ; i++) {
595			o << "(";
596			for (unsigned int j = 0; j < m.getNCol() ; j++) {
597			o << m(i, j);
598			if (j != m.getNCol() -1)
599			o << "\t";
600			}
601			o << ")" << std::endl;
602			}
603			return o;
604			}
605			}
606			#endif //MATH_RECTMATRIX_HPP
607