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