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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
/**
 * @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
     * @param 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
     * @param 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_ = (unsigned int) nrow;
      ncol_ = (unsigned int) ncol;
      data_ = new Real*[nrow_];
      for (unsigned int i = 0; i < nrow_; ++i)
        data_[i] = new Real[ncol_];
    }
    
    void deallocate() {
      for (unsigned 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:	1808
Committed:	Mon Oct 22 20:42:10 2012 UTC (12 years, 6 months ago) by gezelter
File size:	18984 byte(s)
Log Message:	A bug fix in the electric field for the new electrostatic code. Also comment fixes for Doxygen
#	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	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	tim	891	*/
42
43			/**
44			* @file DynamicRectMatrix.hpp
45			* @author Teng Lin
46			* @date 10/11/2004
47			* @version 1.0
48			*/
49
50			#ifndef MATH_DYNAMICRECTMATRIX_HPP
51			#define MATH_DYNAMICRECTMATRIX_HPP
52			#include <math.h>
53			#include <cmath>
54			#include "math/DynamicVector.hpp"
55
56	gezelter	1390	namespace OpenMD {
57	tim	891
58			/**
59			* @class DynamicRectMatrix DynamicRectMatrix.hpp "math/DynamicRectMatrix.hpp"
60			* @brief rectangular matrix class
61			*/
62			template<typename Real>
63			class DynamicRectMatrix {
64			public:
65			typedef Real ElemType;
66			typedef Real* ElemPoinerType;
67			typedef DynamicRectMatrix<Real> SelfType;
68
69			/** default constructor */
70	cli2	1360	DynamicRectMatrix(){
71			nrow_ = 0;
72			ncol_ = 0;
73			data_ = NULL;
74			}
75
76	tim	891	DynamicRectMatrix(int nrow, int ncol) {
77			allocate(nrow, ncol);
78
79			for (unsigned int i = 0; i < nrow_; i++)
80			for (unsigned int j = 0; j < ncol_; j++)
81			this->data_[i][j] = 0.0;
82			}
83
84			/** Constructs and initializes every element of this matrix to a scalar */
85			DynamicRectMatrix(int nrow, int ncol, Real s) {
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] = s;
90			}
91
92			DynamicRectMatrix(int nrow, int ncol, Real* array) {
93			allocate(nrow, ncol);
94			for (unsigned int i = 0; i < nrow_; i++)
95			for (unsigned int j = 0; j < ncol_; j++)
96			this->data_[i][j] = array[i * nrow_ + j];
97			}
98
99			/** copy constructor */
100			DynamicRectMatrix(const SelfType& m) {
101			allocate(m.getNRow(), m.getNCol());
102
103			for (unsigned int i = 0; i < nrow_; i++)
104			for (unsigned int j = 0; j < ncol_; j++)
105			this->data_[i][j] = m.data_[i][j];
106			}
107
108			/** destructor*/
109			~DynamicRectMatrix() { deallocate();}
110
111			/** copy assignment operator */
112			DynamicRectMatrix<Real> operator =(const DynamicRectMatrix<Real> m) {
113			if (this == &m)
114			return *this;
115			if (nrow_ != m.getNRow() \|\| ncol_ != m.getNCol()) {
116			deallocate();
117			allocate(m.getNRow(), m.getNCol());
118			}
119
120			for (unsigned int i = 0; i < nrow_; i++)
121			for (unsigned int j = 0; j < ncol_; j++)
122			this->data_[i][j] = m.data_[i][j];
123			return *this;
124			}
125
126			/**
127			* Return the reference of a single element of this matrix.
128			* @return the reference of a single element of this matrix
129			* @param i row index
130			* @param j Column index
131			*/
132			Real& operator()(unsigned int i, unsigned int j) {
133			return this->data_[i][j];
134			}
135
136			/**
137			* Return the value of a single element of this matrix.
138			* @return the value of a single element of this matrix
139			* @param i row index
140			* @param j Column index
141			*/
142			Real operator()(unsigned int i, unsigned int j) const {
143
144			return this->data_[i][j];
145			}
146
147			/**
148			* Copy the internal data to an array
149			* @param array the pointer of destination array
150			*/
151			void getArray(Real* array) {
152			for (unsigned int i = 0; i < nrow_; i++) {
153			for (unsigned int j = 0; j < ncol_; j++) {
154			array[i * nrow_ + j] = this->data_[i][j];
155			}
156			}
157			}
158
159			/**
160			* Returns a row of this matrix as a vector.
161			* @return a row of this matrix as a vector
162			* @param row the row index
163			*/
164			DynamicVector<Real> getRow(unsigned int row) {
165			DynamicVector<Real> v;
166
167			for (unsigned int i = 0; i < ncol_; i++)
168			v[i] = this->data_[row][i];
169
170			return v;
171			}
172
173			/**
174			* Sets a row of this matrix
175			* @param row the row index
176			* @param v the vector to be set
177			*/
178			void setRow(unsigned int row, const DynamicVector<Real>& v) {
179			assert(v.size() == nrow_);
180			for (unsigned int i = 0; i < ncol_; i++)
181			this->data_[row][i] = v[i];
182			}
183
184			/**
185			* Returns a column of this matrix as a vector.
186			* @return a column of this matrix as a vector
187			* @param col the column index
188			*/
189			DynamicVector<Real> getColumn(unsigned int col) {
190			DynamicVector<Real> v(ncol_);
191
192			for (unsigned int j = 0; j < nrow_; j++)
193			v[j] = this->data_[j][col];
194
195			return v;
196			}
197
198			/**
199			* Sets a column of this matrix
200			* @param col the column index
201			* @param v the vector to be set
202			*/
203			void setColumn(unsigned int col, const DynamicVector<Real>& v){
204
205			for (unsigned int j = 0; j < nrow_; j++)
206			this->data_[j][col] = v[j];
207			}
208
209			/**
210			* swap two rows of this matrix
211			* @param i the first row
212			* @param j the second row
213			*/
214			void swapRow(unsigned int i, unsigned int j){
215			assert(i < nrow_ && j < nrow_);
216
217			for (unsigned int k = 0; k < ncol_; k++)
218			std::swap(this->data_[i][k], this->data_[j][k]);
219			}
220
221			/**
222			* swap two Columns of this matrix
223			* @param i the first Column
224			* @param j the second Column
225			*/
226			void swapColumn(unsigned int i, unsigned int j){
227			assert(i < ncol_ && j < ncol_);
228
229			for (unsigned int k = 0; k < nrow_; k++)
230			std::swap(this->data_[k][i], this->data_[k][j]);
231			}
232
233			/**
234			* Tests if this matrix is identical to matrix m
235			* @return true if this matrix is equal to the matrix m, return false otherwise
236	gezelter	1808	* @param m matrix to be compared
237	tim	891	*
238			* @todo replace operator == by template function equal
239			*/
240			bool operator ==(const DynamicRectMatrix<Real> m) {
241			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
242			for (unsigned int i = 0; i < nrow_; i++)
243			for (unsigned int j = 0; j < ncol_; j++)
244			if (!equal(this->data_[i][j], m.data_[i][j]))
245			return false;
246
247			return true;
248			}
249
250			/**
251			* Tests if this matrix is not equal to matrix m
252			* @return true if this matrix is not equal to the matrix m, return false otherwise
253	gezelter	1808	* @param m matrix to be compared
254	tim	891	*/
255			bool operator !=(const DynamicRectMatrix<Real> m) {
256			return !(*this == m);
257			}
258
259			/** Negates the value of this matrix in place. */
260			inline void negate() {
261			for (unsigned int i = 0; i < nrow_; i++)
262			for (unsigned int j = 0; j < ncol_; j++)
263			this->data_[i][j] = -this->data_[i][j];
264			}
265
266			/**
267			* Sets the value of this matrix to the negation of matrix m.
268			* @param m the source matrix
269			*/
270			inline void negate(const DynamicRectMatrix<Real> m) {
271			for (unsigned int i = 0; i < nrow_; i++)
272			for (unsigned int j = 0; j < ncol_; j++)
273			this->data_[i][j] = -m.data_[i][j];
274			}
275
276			/**
277			* Sets the value of this matrix to the sum of itself and m (*this += m).
278			* @param m the other matrix
279			*/
280			inline void add( const DynamicRectMatrix<Real> m ) {
281			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
282			for (unsigned int i = 0; i < nrow_; i++)
283			for (unsigned int j = 0; j < ncol_; j++)
284			this->data_[i][j] += m.data_[i][j];
285			}
286
287			/**
288			* Sets the value of this matrix to the sum of m1 and m2 (*this = m1 + m2).
289			* @param m1 the first matrix
290			* @param m2 the second matrix
291			*/
292			inline void add( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2 ) {
293			assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
294			for (unsigned int i = 0; i < nrow_; i++)
295			for (unsigned int j = 0; j < ncol_; j++)
296			this->data_[i][j] = m1.data_[i][j] + m2.data_[i][j];
297			}
298
299			/**
300			* Sets the value of this matrix to the difference of itself and m (*this -= m).
301			* @param m the other matrix
302			*/
303			inline void sub( const DynamicRectMatrix<Real> m ) {
304			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
305			for (unsigned int i = 0; i < nrow_; i++)
306			for (unsigned int j = 0; j < ncol_; j++)
307			this->data_[i][j] -= m.data_[i][j];
308			}
309
310			/**
311			* Sets the value of this matrix to the difference of matrix m1 and m2 (*this = m1 - m2).
312			* @param m1 the first matrix
313			* @param m2 the second matrix
314			*/
315			inline void sub( const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2){
316			assert(m1.getNRow() == m2.getNRow() && m1.getNCol() == m2.getNCol());
317			for (unsigned int i = 0; i < nrow_; i++)
318			for (unsigned int j = 0; j < ncol_; j++)
319			this->data_[i][j] = m1.data_[i][j] - m2.data_[i][j];
320			}
321
322			/**
323			* Sets the value of this matrix to the scalar multiplication of itself (this = s).
324			* @param s the scalar value
325			*/
326			inline void mul( Real s ) {
327			for (unsigned int i = 0; i < nrow_; i++)
328			for (unsigned int j = 0; j < ncol_; j++)
329			this->data_[i][j] *= s;
330			}
331
332			/**
333			* Sets the value of this matrix to the scalar multiplication of matrix m (this = s m).
334			* @param s the scalar value
335			* @param m the matrix
336			*/
337			inline void mul( Real s, const DynamicRectMatrix<Real> m ) {
338			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
339			for (unsigned int i = 0; i < nrow_; i++)
340			for (unsigned int j = 0; j < ncol_; j++)
341			this->data_[i][j] = s * m.data_[i][j];
342			}
343
344			/**
345			* Sets the value of this matrix to the scalar division of itself (*this /= s ).
346			* @param s the scalar value
347			*/
348			inline void div( Real s) {
349			for (unsigned int i = 0; i < nrow_; i++)
350			for (unsigned int j = 0; j < ncol_; j++)
351			this->data_[i][j] /= s;
352			}
353
354			/**
355			* Sets the value of this matrix to the scalar division of matrix m (*this = m /s).
356			* @param s the scalar value
357			* @param m the matrix
358			*/
359			inline void div( Real s, const DynamicRectMatrix<Real> m ) {
360			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
361			for (unsigned int i = 0; i < nrow_; i++)
362			for (unsigned int j = 0; j < ncol_; j++)
363			this->data_[i][j] = m.data_[i][j] / s;
364			}
365
366			/**
367			* Multiples a scalar into every element of this matrix.
368			* @param s the scalar value
369			*/
370			DynamicRectMatrix<Real> operator *=(const Real s) {
371			this->mul(s);
372			return *this;
373			}
374
375			/**
376			* Divides every element of this matrix by a scalar.
377			* @param s the scalar value
378			*/
379			DynamicRectMatrix<Real> operator /=(const Real s) {
380			this->div(s);
381			return *this;
382			}
383
384			/**
385			* Sets the value of this matrix to the sum of the other matrix and itself (*this += m).
386			* @param m the other matrix
387			*/
388			DynamicRectMatrix<Real> operator += (const DynamicRectMatrix<Real> m) {
389			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
390			add(m);
391			return *this;
392			}
393
394			/**
395			* Sets the value of this matrix to the differerence of itself and the other matrix (*this -= m)
396			* @param m the other matrix
397			*/
398			DynamicRectMatrix<Real> operator -= (const DynamicRectMatrix<Real> m){
399			assert(nrow_ == m.getNRow() && ncol_ == m.getNCol());
400			sub(m);
401			return *this;
402			}
403
404			/** Return the transpose of this matrix */
405			DynamicRectMatrix<Real> transpose() const{
406			DynamicRectMatrix<Real> result(ncol_,nrow_);
407
408			for (unsigned int i = 0; i < nrow_; i++)
409			for (unsigned int j = 0; j < ncol_; j++)
410			result(j, i) = this->data_[i][j];
411
412			return result;
413			}
414
415			unsigned int getNRow() const {return nrow_;}
416			unsigned int getNCol() const {return ncol_;}
417
418			template<class MatrixType>
419			void setSubMatrix(unsigned int beginRow, unsigned int beginCol, const MatrixType& m) {
420			assert(beginRow + m.getNRow() -1 <= nrow_);
421			assert(beginCol + m.getNCol() -1 <= ncol_);
422
423			for (unsigned int i = 0; i < m.getNRow(); ++i)
424			for (unsigned int j = 0; j < m.getNCol(); ++j)
425			this->data_[beginRow+i][beginCol+j] = m(i, j);
426			}
427
428			template<class MatrixType>
429			void getSubMatrix(unsigned int beginRow, unsigned int beginCol, MatrixType& m) {
430			assert(beginRow + m.getNRow() -1 <= nrow_);
431			assert(beginCol + m.getNCol() - 1 <= ncol_);
432
433			for (unsigned int i = 0; i < m.getNRow(); ++i)
434			for (unsigned int j = 0; j < m.getNCol(); ++j)
435			m(i, j) = this->data_[beginRow+i][beginCol+j];
436			}
437
438			protected:
439			Real** data_;
440			unsigned int nrow_;
441			unsigned int ncol_;
442			private:
443	gezelter	1767	void allocate( int nrow, int ncol ) {
444			nrow_ = (unsigned int) nrow;
445			ncol_ = (unsigned int) ncol;
446			data_ = new Real*[nrow_];
447			for (unsigned int i = 0; i < nrow_; ++i)
448			data_[i] = new Real[ncol_];
449	tim	891	}
450	gezelter	1767
451	tim	891	void deallocate() {
452	gezelter	1767	for (unsigned int i = 0; i < nrow_; ++i)
453	tim	891	delete data_[i];
454			delete []data_;
455	gezelter	1767
456	tim	891	nrow_ = 0;
457			ncol_ = 0;
458			data_ = NULL;
459			}
460
461			};
462
463			/** Negate the value of every element of this matrix. */
464			template<typename Real>
465			inline DynamicRectMatrix<Real> operator -(const DynamicRectMatrix<Real> m) {
466			DynamicRectMatrix<Real> result(m);
467
468			result.negate();
469
470			return result;
471			}
472
473			/**
474			* Return the sum of two matrixes (m1 + m2).
475			* @return the sum of two matrixes
476			* @param m1 the first matrix
477			* @param m2 the second matrix
478			*/
479			template<typename Real>
480			inline DynamicRectMatrix<Real> operator + (const DynamicRectMatrix<Real> m1,const DynamicRectMatrix<Real> m2) {
481
482			DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());
483
484			result.add(m1, m2);
485
486			return result;
487			}
488
489			/**
490			* Return the difference of two matrixes (m1 - m2).
491			* @return the sum of two matrixes
492			* @param m1 the first matrix
493			* @param m2 the second matrix
494			*/
495			template<typename Real>
496			inline DynamicRectMatrix<Real> operator - (const DynamicRectMatrix<Real> m1, const DynamicRectMatrix<Real> m2) {
497			DynamicRectMatrix<Real> result(m1.getNRow(), m1.getNCol());
498
499			result.sub(m1, m2);
500
501			return result;
502			}
503
504			/**
505			* Return the multiplication of scalra and matrix (m * s).
506			* @return the multiplication of a scalra and a matrix
507			* @param m the matrix
508			* @param s the scalar
509			*/
510			template<typename Real>
511			inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real> m, Real s) {
512			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
513
514			result.mul(s, m);
515
516			return result;
517			}
518
519			/**
520			* Return the multiplication of a scalra and a matrix (s * m).
521			* @return the multiplication of a scalra and a matrix
522			* @param s the scalar
523			* @param m the matrix
524			*/
525			template<typename Real>
526			inline DynamicRectMatrix<Real> operator *(Real s, const DynamicRectMatrix<Real> m) {
527			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
528
529			result.mul(s, m);
530
531			return result;
532			}
533
534			/**
535			* Return the multiplication of two matrixes (m1 * m2).
536			* @return the multiplication of two matrixes
537			* @param m1 the first matrix
538			* @param m2 the second matrix
539			*/
540			template<typename Real>
541			inline DynamicRectMatrix<Real> operator *(const DynamicRectMatrix<Real>& m1, const DynamicRectMatrix<Real>& m2) {
542			assert(m1.getNCol() == m2.getNRow());
543			unsigned int sameDim = m1.getNCol();
544			int nrow = m1.getNRow();
545			int ncol = m2.getNCol();
546			DynamicRectMatrix<Real> result(nrow, ncol );
547			for (unsigned int i = 0; i < nrow; i++)
548			for (unsigned int j = 0; j < ncol; j++)
549			for (unsigned int k = 0; k < sameDim; k++)
550			result(i, j) += m1(i, k) * m2(k, j);
551
552			return result;
553			}
554
555			/**
556			* Return the multiplication of a matrix and a vector (m * v).
557			* @return the multiplication of a matrix and a vector
558			* @param m the matrix
559			* @param v the vector
560			*/
561			template<typename Real>
562			inline DynamicVector<Real> operator *(const DynamicRectMatrix<Real> m, const DynamicVector<Real>& v) {
563			int nrow = m.getNRow();
564			int ncol = m.getNCol();
565			assert(ncol = v.size());
566			DynamicVector<Real> result(nrow);
567
568			for (unsigned int i = 0; i < nrow ; i++)
569			for (unsigned int j = 0; j < ncol ; j++)
570			result[i] += m(i, j) * v[j];
571
572			return result;
573			}
574
575			/**
576			* Return the scalar division of matrix (m / s).
577			* @return the scalar division of matrix
578			* @param m the matrix
579			* @param s the scalar
580			*/
581			template<typename Real>
582			inline DynamicRectMatrix<Real> operator /(const DynamicRectMatrix<Real> m, Real s) {
583			DynamicRectMatrix<Real> result(m.getNRow(), m.getNCol());
584
585			result.div(s, m);
586
587			return result;
588			}
589
590			/**
591			* Write to an output stream
592			*/
593			template<typename Real>
594			std::ostream &operator<< ( std::ostream& o, const DynamicRectMatrix<Real> m) {
595			for (unsigned int i = 0; i < m.getNRow() ; i++) {
596			o << "(";
597			for (unsigned int j = 0; j < m.getNCol() ; j++) {
598			o << m(i, j);
599			if (j != m.getNCol() -1)
600			o << "\t";
601			}
602			o << ")" << std::endl;
603			}
604			return o;
605			}
606			}
607			#endif //MATH_RECTMATRIX_HPP
608