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. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. 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.
 */
 
/**
 * @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 oopse {

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