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, 8 months ago) by tim
File size:	18729 byte(s)
Log Message:	Adding Hydrodynamics Module
#	Content
1	/*
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