src/math/RectMatrix.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 RectMatrix.hpp
 * @author Teng Lin
 * @date 10/11/2004
 * @version 1.0
 */

#ifndef MATH_RECTMATRIX_HPP
#define MATH_RECTMATRIX_HPP
#include <math.h>
#include <cmath>
#include "Vector.hpp"

namespace OpenMD {

  /**
   * @class RectMatrix RectMatrix.hpp "math/RectMatrix.hpp"
   * @brief rectangular matrix class
   */
  template<typename Real, unsigned int Row, unsigned int Col> 
  class RectMatrix {
  public:
    typedef Real ElemType;
    typedef Real* ElemPoinerType;
            
    /** default constructor */
    RectMatrix() {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; j++)
          this->data_[i][j] = 0.0;
    }

    /** Constructs and initializes every element of this matrix to a scalar */ 
    RectMatrix(Real s) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; j++)
          this->data_[i][j] = s;
    }

    RectMatrix(Real* array) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; j++)
          this->data_[i][j] = array[i * Row + j];
    }

    /** copy constructor */
    RectMatrix(const RectMatrix<Real, Row, Col>& m) {
      *this = m;
    }
            
    /** destructor*/
    ~RectMatrix() {}

    /** copy assignment operator */
    RectMatrix<Real, Row, Col>& operator =(const RectMatrix<Real, Row, Col>& m) {
      if (this == &m)
        return *this;
                
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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) {
      //assert( i < Row && j < Col);
      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 < Row; i++) {
        for (unsigned int j = 0; j < Col; j++) {
          array[i * Row + j] = this->data_[i][j];
        }
      }
    }


    /** Returns the pointer of internal array */
    Real* getArrayPointer() {
      return &this->data_[0][0];
    }

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

      for (unsigned int i = 0; i < Col; 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 Vector<Real, Row>& v) {

      for (unsigned int i = 0; i < Col; 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
     */                
    Vector<Real, Col> getColumn(unsigned int col) {
      Vector<Real, Col> v;

      for (unsigned int j = 0; j < Row; 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 Vector<Real, Col>& v){

      for (unsigned int j = 0; j < Row; 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 < Row && j < Row);

      for (unsigned int k = 0; k < Col; 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 < Col && j < Col);
                    
      for (unsigned int k = 0; k < Row; 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 RectMatrix<Real, Row, Col>& m) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m) {
      return !(*this == m);
    }

    /** Negates the value of this matrix in place. */           
    inline void negate() {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m ) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2 ) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m ) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2){
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m ) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; 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 < Row; i++)
        for (unsigned int j = 0; j < Col; 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 RectMatrix<Real, Row, Col>& m ) {
      for (unsigned int i = 0; i < Row; i++)
        for (unsigned int j = 0; j < Col; j++)  
          this->data_[i][j] = m.data_[i][j] / s;
    }

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

    /**
     *  Divides every element of this matrix by a scalar.
     * @param s the scalar value
     */
    RectMatrix<Real, Row, Col>& 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
     */
    RectMatrix<Real, Row, Col>& operator += (const RectMatrix<Real, Row, Col>& m) {
      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
     */
    RectMatrix<Real, Row, Col>& operator -= (const RectMatrix<Real, Row, Col>& m){
      sub(m);
      return *this;
    }

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

      return result;
    }

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

        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 <= getNRow());
        assert(beginCol + m.getNCol() - 1 <= getNCol());

        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];
    }
    
    unsigned int getNRow() const {return Row;}
    unsigned int getNCol() const {return Col;}        

  protected:
    Real data_[Row][Col];
  };

  /** Negate the value of every element of this matrix. */
  template<typename Real, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator -(const RectMatrix<Real, Row, Col>& m) {
    RectMatrix<Real, Row, Col> 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, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator + (const RectMatrix<Real, Row, Col>& m1,const RectMatrix<Real, Row, Col>& m2) {
    RectMatrix<Real, Row, Col> result;

    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, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator - (const RectMatrix<Real, Row, Col>& m1, const RectMatrix<Real, Row, Col>& m2) {
    RectMatrix<Real, Row, Col> result;

    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, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, Col>& m, Real s) {
    RectMatrix<Real, Row, Col> result;

    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, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator *(Real s, const RectMatrix<Real, Row, Col>& m) {
    RectMatrix<Real, Row, Col> result;

    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, unsigned int Row, unsigned int Col, unsigned int SameDim> 
  inline RectMatrix<Real, Row, Col> operator *(const RectMatrix<Real, Row, SameDim>& m1, const RectMatrix<Real, SameDim, Col>& m2) {
    RectMatrix<Real, Row, Col> result;

    for (unsigned int i = 0; i < Row; i++)
      for (unsigned int j = 0; j < Col; j++)
        for (unsigned int k = 0; k < SameDim; k++)
          result(i, j)  += m1(i, k) * m2(k, j);                

    return result;
  }
    
  /**
   * Returns 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, unsigned int Row, unsigned int Col>
  inline Vector<Real, Row> operator *(const RectMatrix<Real, Row, Col>& m, const Vector<Real, Col>& v) {
    Vector<Real, Row> result;

    for (unsigned int i = 0; i < Row ; i++)
      for (unsigned int j = 0; j < Col ; j++)            
        result[i] += m(i, j) * v[j];
            
    return result;                                                                 
  }

  /**
   * Returns the multiplication of a vector transpose and a matrix  (v^T * m). 
   * @return the multiplication of a vector transpose and a matrix
   * @param v the vector
   * @param m the matrix
   */
  template<typename Real, unsigned int Row, unsigned int Col>
  inline Vector<Real, Col> operator *(const Vector<Real, Row>& v, const RectMatrix<Real, Row, Col>& m) {
    Vector<Real, Row> result;
    
    for (unsigned int i = 0; i < Col ; i++)
      for (unsigned int j = 0; j < Row ; j++)            
        result[i] += v[j] * m(j, i);
            
    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, unsigned int Row, unsigned int Col> 
  inline RectMatrix<Real, Row, Col> operator /(const RectMatrix<Real, Row, Col>& m, Real s) {
    RectMatrix<Real, Row, Col> result;

    result.div(s, m);

    return result;
  }    

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