--- trunk/src/math/RectMatrix.hpp 2005/04/15 22:04:00 507 +++ trunk/src/math/RectMatrix.hpp 2013/08/23 15:59:23 1933 @@ -6,19 +6,10 @@ * 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 + * 1. 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 + * 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. @@ -37,6 +28,16 @@ * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -52,7 +53,7 @@ #include #include "Vector.hpp" -namespace oopse { +namespace OpenMD { /** * @class RectMatrix RectMatrix.hpp "math/RectMatrix.hpp" @@ -151,7 +152,7 @@ namespace oopse { Vector getRow(unsigned int row) { Vector v; - for (unsigned int i = 0; i < Row; i++) + for (unsigned int i = 0; i < Col; i++) v[i] = this->data_[row][i]; return v; @@ -164,7 +165,7 @@ namespace oopse { */ void setRow(unsigned int row, const Vector& v) { - for (unsigned int i = 0; i < Row; i++) + for (unsigned int i = 0; i < Col; i++) this->data_[row][i] = v[i]; } @@ -176,7 +177,7 @@ namespace oopse { Vector getColumn(unsigned int col) { Vector v; - for (unsigned int j = 0; j < Col; j++) + for (unsigned int j = 0; j < Row; j++) v[j] = this->data_[j][col]; return v; @@ -189,7 +190,7 @@ namespace oopse { */ void setColumn(unsigned int col, const Vector& v){ - for (unsigned int j = 0; j < Col; j++) + for (unsigned int j = 0; j < Row; j++) this->data_[j][col] = v[j]; } @@ -220,7 +221,7 @@ namespace oopse { /** * 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 + * @param m matrix to be compared * * @todo replace operator == by template function equal */ @@ -236,7 +237,7 @@ namespace oopse { /** * 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 + * @param m matrix to be compared */ bool operator !=(const RectMatrix& m) { return !(*this == m); @@ -389,7 +390,30 @@ namespace oopse { return result; } - + + template + 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 + 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]; }; @@ -483,7 +507,7 @@ namespace oopse { } /** - * Return the multiplication of a matrix and a vector (m * v). + * 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 @@ -495,6 +519,23 @@ namespace oopse { 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 + inline Vector operator *(const Vector& v, const RectMatrix& m) { + Vector 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; } @@ -514,7 +555,72 @@ namespace oopse { return result; } + + /** + * Returns the tensor contraction (double dot product) of two rank 2 + * tensors (or Matrices) + * + * \f[ \mathbf{A} \colon \! \mathbf{B} = \sum_\alpha \sum_\beta \mathbf{A}_{\alpha \beta} B_{\alpha \beta} \f] + * + * @param t1 first tensor + * @param t2 second tensor + * @return the tensor contraction (double dot product) of t1 and t2 + */ + template + inline Real doubleDot( const RectMatrix& t1, + const RectMatrix& t2 ) { + Real tmp; + tmp = 0; + + for (unsigned int i = 0; i < Row; i++) + for (unsigned int j =0; j < Col; j++) + tmp += t1(i,j) * t2(i,j); + + return tmp; + } + + + /** + * 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 + inline Vector mCross( const RectMatrix& t1, + const RectMatrix& t2 ) { + Vector 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