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(){
      nrow_ = 0;
      ncol_ = 0;
      data_ = NULL;
    }

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

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

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

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

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

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

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

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

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

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

      return v;
    }

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

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

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

      return v;
    }

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

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

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

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

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

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

      return true;
    }

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

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

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

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

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

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

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

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

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

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

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

      return result;
    }

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

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

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

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

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

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

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

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

    result.negate();

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

    result.add(m1, m2);

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

    result.sub(m1, m2);

    return result;
  }

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

    result.mul(s, m);

    return result;
  }

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

    result.mul(s, m);

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

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

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

    result.div(s, m);

    return result;
  }    

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

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