--- trunk/src/math/RectMatrix.hpp 2004/10/14 23:28:09 76 +++ trunk/src/math/RectMatrix.hpp 2005/03/01 20:10:14 385 @@ -1,29 +1,44 @@ -/* - * Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project - * - * Contact: oopse@oopse.org - * - * This program is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public License - * as published by the Free Software Foundation; either version 2.1 - * of the License, or (at your option) any later version. - * All we ask is that proper credit is given for our work, which includes - * - but is not limited to - adding the above copyright notice to the beginning - * of your source code files, and to any copyright notice that you may distribute - * with programs based on this work. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + /* + * 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 RectMatrix.hpp * @author Teng Lin @@ -33,7 +48,7 @@ #ifndef MATH_RECTMATRIX_HPP #define MATH_RECTMATRIX_HPP - +#include #include #include "Vector.hpp" @@ -46,284 +61,334 @@ namespace oopse { template 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; + } - /** default constructor */ - RectMatrix() { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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; + } - /** 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++) - 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& m) { - *this = m; - } - - /** destructor*/ - ~RectMatrix() {} + /** copy constructor */ + RectMatrix(const RectMatrix& m) { + *this = m; + } + + /** destructor*/ + ~RectMatrix() {} - /** copy assignment operator */ - RectMatrix& operator =(const RectMatrix& m) { - if (this == &m) + /** copy assignment operator */ + RectMatrix& operator =(const RectMatrix& 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; + } - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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 colum index - */ - double& operator()(unsigned int i, unsigned int j) { - //assert( i < Row && j < Col); - return data_[i][j]; - } + /** + * 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 colum index - */ - double operator()(unsigned int i, unsigned int j) const { - - return 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]; + } - /** - * Returns a row of this matrix as a vector. - * @return a row of this matrix as a vector - * @param row the row index - */ - Vector getRow(unsigned int row) { - Vector v; + /** + * 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]; + } + } + } - for (unsigned int i = 0; i < Row; i++) - v[i] = data_[row][i]; - return v; - } + /** Returns the pointer of internal array */ + Real* getArrayPointer() { + return &this->data_[0][0]; + } - /** - * 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& v) { + /** + * Returns a row of this matrix as a vector. + * @return a row of this matrix as a vector + * @param row the row index + */ + Vector getRow(unsigned int row) { + Vector v; - for (unsigned int i = 0; i < Row; i++) - data_[row][i] = v[i]; - } + for (unsigned int i = 0; i < Row; i++) + v[i] = this->data_[row][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 getColum(unsigned int col) { - Vector v; + return v; + } - for (unsigned int j = 0; j < Col; j++) - v[j] = data_[j][col]; + /** + * 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& v) { - return v; - } + for (unsigned int i = 0; i < Row; i++) + this->data_[row][i] = v[i]; + } - /** - * Sets a column of this matrix - * @param col the column index - * @param v the vector to be set - */ - void setColum(unsigned int col, const Vector& v){ + /** + * Returns a column of this matrix as a vector. + * @return a column of this matrix as a vector + * @param col the column index + */ + Vector getColumn(unsigned int col) { + Vector v; - for (unsigned int j = 0; j < Col; j++) - data_[j][col] = v[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 RectMatrix& m) { - for (unsigned int i = 0; i < Row; i++) for (unsigned int j = 0; j < Col; j++) - if (!equal(data_[i][j], m.data_[i][j])) - return false; + v[j] = this->data_[j][col]; - return true; - } + return v; + } - /** - * 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 RectMatrix& m) { - return !(*this == m); - } + /** + * 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& v){ - /** 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++) - data_[i][j] = -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& m) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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& m ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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& m1, const RectMatrix& m2 ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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& m ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - 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& m1, const RectMatrix& m2){ - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - data_[i][j] = m1.data_[i][j] - m2.data_[i][j]; - } + this->data_[j][col] = v[j]; + } - /** - * Sets the value of this matrix to the scalar multiplication of itself (*this *= s). - * @param s the scalar value - */ - inline void mul( double s ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - data_[i][j] *= s; - } + /** + * 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); - /** - * 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( double s, const RectMatrix& m ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - data_[i][j] = s * m.data_[i][j]; - } + for (unsigned int k = 0; k < Col; k++) + std::swap(this->data_[i][k], this->data_[j][k]); + } - /** - * Sets the value of this matrix to the scalar division of itself (*this /= s ). - * @param s the scalar value - */ - inline void div( double s) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - data_[i][j] /= s; - } + /** + * 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]); + } - /** - * 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( double s, const RectMatrix& m ) { - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - data_[i][j] = m.data_[i][j] / s; - } + /** + * 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 RectMatrix& 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; - /** - * Multiples a scalar into every element of this matrix. - * @param s the scalar value - */ - RectMatrix& operator *=(const double s) { - this->mul(s); - return *this; - } + return true; + } - /** - * Divides every element of this matrix by a scalar. - * @param s the scalar value - */ - RectMatrix& operator /=(const double s) { - this->div(s); - return *this; - } + /** + * 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 RectMatrix& m) { + return !(*this == m); + } - /** - * Sets the value of this matrix to the sum of the other matrix and itself (*this += m). - * @param m the other matrix - */ - RectMatrix& operator += (const RectMatrix& 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& operator -= (const RectMatrix& m){ - sub(m); - return *this; - } - - /** Return the transpose of this matrix */ - RectMatrix transpose(){ - RectMatrix result; + /** 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]; + } - for (unsigned int i = 0; i < Row; i++) - for (unsigned int j = 0; j < Col; j++) - result(j, i) = 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& 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& 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& m1, const RectMatrix& 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& 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& m1, const RectMatrix& 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]; + } - return result; - } + /** + * 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& 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& 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& operator *=(const Real s) { + this->mul(s); + return *this; + } + + /** + * Divides every element of this matrix by a scalar. + * @param s the scalar value + */ + RectMatrix& 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& operator += (const RectMatrix& 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& operator -= (const RectMatrix& m){ + sub(m); + return *this; + } + + /** Return the transpose of this matrix */ + RectMatrix transpose() const{ + RectMatrix 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; + } protected: Real data_[Row][Col]; @@ -448,5 +513,22 @@ namespace oopse { return result; } + + /** + * Write to an output stream + */ + template + std::ostream &operator<< ( std::ostream& o, const RectMatrix& 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