| 1 | < | /* | 
| 2 | < | * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project | 
| 3 | < | * | 
| 4 | < | * Contact: oopse@oopse.org | 
| 5 | < | * | 
| 6 | < | * This program is free software; you can redistribute it and/or | 
| 7 | < | * modify it under the terms of the GNU Lesser General Public License | 
| 8 | < | * as published by the Free Software Foundation; either version 2.1 | 
| 9 | < | * of the License, or (at your option) any later version. | 
| 10 | < | * All we ask is that proper credit is given for our work, which includes | 
| 11 | < | * - but is not limited to - adding the above copyright notice to the beginning | 
| 12 | < | * of your source code files, and to any copyright notice that you may distribute | 
| 13 | < | * with programs based on this work. | 
| 14 | < | * | 
| 15 | < | * This program is distributed in the hope that it will be useful, | 
| 16 | < | * but WITHOUT ANY WARRANTY; without even the implied warranty of | 
| 17 | < | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the | 
| 18 | < | * GNU Lesser General Public License for more details. | 
| 19 | < | * | 
| 20 | < | * You should have received a copy of the GNU Lesser General Public License | 
| 21 | < | * along with this program; if not, write to the Free Software | 
| 22 | < | * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA. | 
| 1 | > | /* | 
| 2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  | * | 
| 4 | + | * The University of Notre Dame grants you ("Licensee") a | 
| 5 | + | * non-exclusive, royalty free, license to use, modify and | 
| 6 | + | * redistribute this software in source and binary code form, provided | 
| 7 | + | * that the following conditions are met: | 
| 8 | + | * | 
| 9 | + | * 1. Acknowledgement of the program authors must be made in any | 
| 10 | + | *    publication of scientific results based in part on use of the | 
| 11 | + | *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | + | *    the article in which the program was described (Matthew | 
| 13 | + | *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | + | *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | + | *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | + | *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | + | * | 
| 18 | + | * 2. Redistributions of source code must retain the above copyright | 
| 19 | + | *    notice, this list of conditions and the following disclaimer. | 
| 20 | + | * | 
| 21 | + | * 3. Redistributions in binary form must reproduce the above copyright | 
| 22 | + | *    notice, this list of conditions and the following disclaimer in the | 
| 23 | + | *    documentation and/or other materials provided with the | 
| 24 | + | *    distribution. | 
| 25 | + | * | 
| 26 | + | * This software is provided "AS IS," without a warranty of any | 
| 27 | + | * kind. All express or implied conditions, representations and | 
| 28 | + | * warranties, including any implied warranty of merchantability, | 
| 29 | + | * fitness for a particular purpose or non-infringement, are hereby | 
| 30 | + | * excluded.  The University of Notre Dame and its licensors shall not | 
| 31 | + | * be liable for any damages suffered by licensee as a result of | 
| 32 | + | * using, modifying or distributing the software or its | 
| 33 | + | * derivatives. In no event will the University of Notre Dame or its | 
| 34 | + | * licensors be liable for any lost revenue, profit or data, or for | 
| 35 | + | * direct, indirect, special, consequential, incidental or punitive | 
| 36 | + | * damages, however caused and regardless of the theory of liability, | 
| 37 | + | * arising out of the use of or inability to use software, even if the | 
| 38 | + | * University of Notre Dame has been advised of the possibility of | 
| 39 | + | * such damages. | 
| 40 |  | */ | 
| 41 | < |  | 
| 41 | > |  | 
| 42 |  | /** | 
| 43 |  | * @file Vector.hpp | 
| 44 |  | * @author Teng Lin | 
| 52 |  | #include <cassert> | 
| 53 |  | #include <cmath> | 
| 54 |  | #include <iostream> | 
| 55 | < |  | 
| 55 | > | #include <math.h> | 
| 56 |  | namespace oopse { | 
| 57 |  |  | 
| 58 | < | const double epsilon = 0.000001; | 
| 58 | > | static const double epsilon = 0.000001; | 
| 59 |  |  | 
| 60 |  | template<typename T> | 
| 61 |  | inline bool equal(T e1, T e2) { | 
| 71 |  | inline bool equal(double e1, double e2) { | 
| 72 |  | return fabs(e1 - e2) < epsilon; | 
| 73 |  | } | 
| 74 | + |  | 
| 75 |  |  | 
| 76 |  | /** | 
| 77 |  | * @class Vector Vector.hpp "math/Vector.hpp" | 
| 81 |  | class Vector{ | 
| 82 |  | public: | 
| 83 |  |  | 
| 84 | + | typedef Real ElemType; | 
| 85 | + | typedef Real* ElemPoinerType; | 
| 86 | + |  | 
| 87 |  | /** default constructor */ | 
| 88 |  | inline Vector(){ | 
| 89 |  | for (unsigned int i = 0; i < Dim; i++) | 
| 90 | < | data_[i] = 0.0; | 
| 90 | > | this->data_[i] = 0; | 
| 91 |  | } | 
| 92 |  |  | 
| 93 |  | /** Constructs and initializes a Vector from a vector */ | 
| 101 |  | return *this; | 
| 102 |  |  | 
| 103 |  | for (unsigned int i = 0; i < Dim; i++) | 
| 104 | < | data_[i] = v[i]; | 
| 104 | > | this->data_[i] = v[i]; | 
| 105 |  |  | 
| 106 |  | return *this; | 
| 107 |  | } | 
| 108 | + |  | 
| 109 | + | template<typename T> | 
| 110 | + | inline Vector(const T& s){ | 
| 111 | + | for (unsigned int i = 0; i < Dim; i++) | 
| 112 | + | this->data_[i] = s; | 
| 113 | + | } | 
| 114 |  |  | 
| 115 |  | /** Constructs and initializes a Vector from an array */ | 
| 116 | < | inline Vector( double* v) { | 
| 117 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 118 | < | data_[i] = v[i]; | 
| 116 | > | inline Vector( Real* v) { | 
| 117 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 118 | > | this->data_[i] = v[i]; | 
| 119 |  | } | 
| 120 |  |  | 
| 121 |  | /** | 
| 123 |  | * @return reference of ith element | 
| 124 |  | * @param i index | 
| 125 |  | */ | 
| 126 | < | inline double& operator[](unsigned int  i) { | 
| 126 | > | inline Real& operator[](unsigned int  i) { | 
| 127 |  | assert( i < Dim); | 
| 128 | < | return data_[i]; | 
| 128 | > | return this->data_[i]; | 
| 129 |  | } | 
| 130 |  |  | 
| 131 |  | /** | 
| 133 |  | * @return reference of ith element | 
| 134 |  | * @param i index | 
| 135 |  | */ | 
| 136 | < | inline double& operator()(unsigned int  i) { | 
| 136 | > | inline Real& operator()(unsigned int  i) { | 
| 137 |  | assert( i < Dim); | 
| 138 | < | return data_[i]; | 
| 138 | > | return this->data_[i]; | 
| 139 |  | } | 
| 140 |  |  | 
| 141 |  | /** | 
| 143 |  | * @return reference of ith element | 
| 144 |  | * @param i index | 
| 145 |  | */ | 
| 146 | < | inline  const double& operator[](unsigned int i) const { | 
| 146 | > | inline  const Real& operator[](unsigned int i) const { | 
| 147 |  | assert( i < Dim); | 
| 148 | < | return data_[i]; | 
| 148 | > | return this->data_[i]; | 
| 149 |  | } | 
| 150 |  |  | 
| 151 |  | /** | 
| 153 |  | * @return reference of ith element | 
| 154 |  | * @param i index | 
| 155 |  | */ | 
| 156 | < | inline  const double& operator()(unsigned int i) const { | 
| 156 | > | inline  const Real& operator()(unsigned int i) const { | 
| 157 |  | assert( i < Dim); | 
| 158 | < | return data_[i]; | 
| 158 | > | return this->data_[i]; | 
| 159 |  | } | 
| 160 |  |  | 
| 161 | + | /** Copy the internal data to an array*/ | 
| 162 | + | void getArray(Real* array) { | 
| 163 | + | for (unsigned int i = 0; i < Dim; i ++) { | 
| 164 | + | array[i] = this->data_[i]; | 
| 165 | + | } | 
| 166 | + | } | 
| 167 | + |  | 
| 168 | + | /** Returns the pointer of internal array */ | 
| 169 | + | Real* getArrayPointer() { | 
| 170 | + | return this->data_; | 
| 171 | + | } | 
| 172 | + |  | 
| 173 |  | /** | 
| 174 |  | * Tests if this vetor is equal to other vector | 
| 175 |  | * @return true if equal, otherwise return false | 
| 178 |  | inline bool operator ==(const Vector<Real, Dim>& v) { | 
| 179 |  |  | 
| 180 |  | for (unsigned int i = 0; i < Dim; i ++) { | 
| 181 | < | if (!equal(data_[i], v[i])) { | 
| 181 | > | if (!equal(this->data_[i], v[i])) { | 
| 182 |  | return false; | 
| 183 |  | } | 
| 184 |  | } | 
| 197 |  |  | 
| 198 |  | /** Negates the value of this vector in place. */ | 
| 199 |  | inline void negate() { | 
| 200 | < | data_[0] = -data_[0]; | 
| 201 | < | data_[1] = -data_[1]; | 
| 164 | < | data_[2] = -data_[2]; | 
| 200 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 201 | > | this->data_[i] = -this->data_[i]; | 
| 202 |  | } | 
| 203 |  |  | 
| 204 |  | /** | 
| 207 |  | */ | 
| 208 |  | inline void negate(const Vector<Real, Dim>& v1) { | 
| 209 |  | for (unsigned int i = 0; i < Dim; i++) | 
| 210 | < | data_[i] = -v1.data_[i]; | 
| 210 | > | this->data_[i] = -v1.data_[i]; | 
| 211 |  |  | 
| 212 |  | } | 
| 213 |  |  | 
| 216 |  | * @param v1 the other vector | 
| 217 |  | */ | 
| 218 |  | inline void add( const Vector<Real, Dim>& v1 ) { | 
| 219 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 220 | < | data_[i] += v1.data_[i]; | 
| 221 | < | } | 
| 219 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 220 | > | this->data_[i] += v1.data_[i]; | 
| 221 | > | } | 
| 222 |  |  | 
| 223 |  | /** | 
| 224 |  | * Sets the value of this vector to the sum of v1 and v2 (*this = v1 + v2). | 
| 226 |  | * @param v2 the second vector | 
| 227 |  | */ | 
| 228 |  | inline void add( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) { | 
| 229 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 230 | < | data_[i] = v1.data_[i] + v2.data_[i]; | 
| 229 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 230 | > | this->data_[i] = v1.data_[i] + v2.data_[i]; | 
| 231 |  | } | 
| 232 |  |  | 
| 233 |  | /** | 
| 236 |  | */ | 
| 237 |  | inline void sub( const Vector<Real, Dim>& v1 ) { | 
| 238 |  | for (unsigned int i = 0; i < Dim; i++) | 
| 239 | < | data_[i] -= v1.data_[i]; | 
| 239 | > | this->data_[i] -= v1.data_[i]; | 
| 240 |  | } | 
| 241 |  |  | 
| 242 |  | /** | 
| 246 |  | */ | 
| 247 |  | inline void sub( const Vector<Real, Dim>& v1, const Vector  &v2 ){ | 
| 248 |  | for (unsigned int i = 0; i < Dim; i++) | 
| 249 | < | data_[i] = v1.data_[i] - v2.data_[i]; | 
| 249 | > | this->data_[i] = v1.data_[i] - v2.data_[i]; | 
| 250 |  | } | 
| 251 |  |  | 
| 252 |  | /** | 
| 253 |  | * Sets the value of this vector to the scalar multiplication of itself (*this *= s). | 
| 254 |  | * @param s the scalar value | 
| 255 |  | */ | 
| 256 | < | inline void mul( double s ) { | 
| 257 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 258 | < | data_[i] *= s; | 
| 256 | > | inline void mul( Real s ) { | 
| 257 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 258 | > | this->data_[i] *= s; | 
| 259 |  | } | 
| 260 |  |  | 
| 261 |  | /** | 
| 262 |  | * Sets the value of this vector to the scalar multiplication of vector v1 | 
| 263 |  | * (*this = s * v1). | 
| 264 | + | * @param v1 the vector | 
| 265 |  | * @param s the scalar value | 
| 228 | – | * @param v1 the vector | 
| 266 |  | */ | 
| 267 | < | inline void mul( double s, const Vector<Real, Dim>& v1 ) { | 
| 268 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 269 | < | data_[i] = s * v1.data_[i]; | 
| 267 | > | inline void mul( const Vector<Real, Dim>& v1, Real s) { | 
| 268 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 269 | > | this->data_[i] = s * v1.data_[i]; | 
| 270 |  | } | 
| 271 |  |  | 
| 272 |  | /** | 
| 273 |  | * Sets the value of this vector to the scalar division of itself  (*this /= s ). | 
| 274 |  | * @param s the scalar value | 
| 275 |  | */ | 
| 276 | < | inline void div( double s) { | 
| 277 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 278 | < | data_[i] /= s; | 
| 276 | > | inline void div( Real s) { | 
| 277 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 278 | > | this->data_[i] /= s; | 
| 279 |  | } | 
| 280 |  |  | 
| 281 |  | /** | 
| 283 |  | * @param v1 the source vector | 
| 284 |  | * @param s the scalar value | 
| 285 |  | */ | 
| 286 | < | inline void div( const Vector<Real, Dim>& v1, double s ) { | 
| 287 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 288 | < | data_[i] = v1.data_[i] / s; | 
| 286 | > | inline void div( const Vector<Real, Dim>& v1, Real s ) { | 
| 287 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 288 | > | this->data_[i] = v1.data_[i] / s; | 
| 289 |  | } | 
| 290 |  |  | 
| 291 |  | /** @see #add */ | 
| 301 |  | } | 
| 302 |  |  | 
| 303 |  | /** @see #mul */ | 
| 304 | < | inline Vector<Real, Dim>& operator *=( double s) { | 
| 304 | > | inline Vector<Real, Dim>& operator *=( Real s) { | 
| 305 |  | mul(s); | 
| 306 |  | return *this; | 
| 307 |  | } | 
| 308 |  |  | 
| 309 |  | /** @see #div */ | 
| 310 | < | inline Vector<Real, Dim>& operator /=( double s ) { | 
| 310 | > | inline Vector<Real, Dim>& operator /=( Real s ) { | 
| 311 |  | div(s); | 
| 312 |  | return *this; | 
| 313 |  | } | 
| 316 |  | * Returns the length of this vector. | 
| 317 |  | * @return the length of this vector | 
| 318 |  | */ | 
| 319 | < | inline double length() { | 
| 320 | < | return sqrt(lengthSquared()); | 
| 319 | > | inline Real length() { | 
| 320 | > | return sqrt(lengthSquare()); | 
| 321 |  | } | 
| 322 |  |  | 
| 323 |  | /** | 
| 324 |  | * Returns the squared length of this vector. | 
| 325 |  | * @return the squared length of this vector | 
| 326 |  | */ | 
| 327 | < | inline double lengthSquared() { | 
| 327 | > | inline Real lengthSquare() { | 
| 328 |  | return dot(*this, *this); | 
| 329 |  | } | 
| 330 |  |  | 
| 331 |  | /** Normalizes this vector in place */ | 
| 332 |  | inline void normalize() { | 
| 333 | < | double len; | 
| 333 | > | Real len; | 
| 334 |  |  | 
| 335 |  | len = length(); | 
| 336 | + |  | 
| 337 | + | //if (len < oopse:epsilon) | 
| 338 | + | //  throw(); | 
| 339 | + |  | 
| 340 |  | *this /= len; | 
| 341 |  | } | 
| 342 | + |  | 
| 343 | + | /** | 
| 344 | + | * Tests if this vector is normalized | 
| 345 | + | * @return true if this vector is normalized, otherwise return false | 
| 346 | + | */ | 
| 347 | + | inline bool isNormalized() { | 
| 348 | + | return equal(lengthSquare(), 1.0); | 
| 349 | + | } | 
| 350 |  |  | 
| 351 |  | protected: | 
| 352 | < | double data_[3]; | 
| 352 | > | Real data_[Dim]; | 
| 353 |  |  | 
| 354 |  | }; | 
| 355 |  |  | 
| 356 |  | /** unary minus*/ | 
| 357 |  | template<typename Real, unsigned int Dim> | 
| 358 |  | inline Vector<Real, Dim> operator -(const Vector<Real, Dim>& v1){ | 
| 359 | < | Vector tmp(v1); | 
| 360 | < | return tmp.negate(); | 
| 359 | > | Vector<Real, Dim> tmp(v1); | 
| 360 | > | tmp.negate(); | 
| 361 | > | return tmp; | 
| 362 |  | } | 
| 363 |  |  | 
| 364 |  | /** | 
| 395 |  | * @param s the scalar value | 
| 396 |  | */ | 
| 397 |  | template<typename Real, unsigned int Dim> | 
| 398 | < | Vector<Real, Dim> operator * ( const Vector<Real, Dim>& v1, double s) { | 
| 398 | > | Vector<Real, Dim> operator * ( const Vector<Real, Dim>& v1, Real s) { | 
| 399 |  | Vector<Real, Dim> result; | 
| 400 | < | result.mul(s, v1); | 
| 400 | > | result.mul(v1,s); | 
| 401 |  | return result; | 
| 402 |  | } | 
| 403 |  |  | 
| 408 |  | * @param v1 the source vector | 
| 409 |  | */ | 
| 410 |  | template<typename Real, unsigned int Dim> | 
| 411 | < | Vector<Real, Dim> operator * ( double s, const Vector<Real, Dim>& v1 ) { | 
| 411 | > | Vector<Real, Dim> operator * ( Real s, const Vector<Real, Dim>& v1 ) { | 
| 412 |  | Vector<Real, Dim> result; | 
| 413 | < | result.mul(s, v1); | 
| 413 | > | result.mul(v1, s); | 
| 414 |  | return result; | 
| 415 |  | } | 
| 416 |  |  | 
| 421 |  | * @param s the scalar value | 
| 422 |  | */ | 
| 423 |  | template<typename Real, unsigned int Dim> | 
| 424 | < | Vector<Real, Dim> operator / ( const Vector<Real, Dim>& v1, double s) { | 
| 424 | > | Vector<Real, Dim> operator / ( const Vector<Real, Dim>& v1, Real s) { | 
| 425 |  | Vector<Real, Dim> result; | 
| 426 |  | result.div( v1,s); | 
| 427 |  | return result; | 
| 428 |  | } | 
| 429 |  |  | 
| 430 |  | /** | 
| 381 | – | * Returns the  value of division of a vector by a scalar. | 
| 382 | – | * @return  the vaule of scalar division of this vector | 
| 383 | – | * @param s the scalar value | 
| 384 | – | * @param v1 the source vector | 
| 385 | – | */ | 
| 386 | – | template<typename Real, unsigned int Dim> | 
| 387 | – | inline Vector<Real, Dim> operator /( double s, const Vector<Real, Dim>& v1 ) { | 
| 388 | – | Vector<Real, Dim> result; | 
| 389 | – | result.div( v1,s); | 
| 390 | – | return result; | 
| 391 | – | } | 
| 392 | – |  | 
| 393 | – | /** fuzzy comparson */ | 
| 394 | – | template<typename Real, unsigned int Dim> | 
| 395 | – | inline bool epsilonEqual( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) { | 
| 396 | – |  | 
| 397 | – | } | 
| 398 | – |  | 
| 399 | – |  | 
| 400 | – | /** | 
| 431 |  | * Returns the dot product of two Vectors | 
| 432 |  | * @param v1 first vector | 
| 433 |  | * @param v2 second vector | 
| 435 |  | */ | 
| 436 |  | template<typename Real, unsigned int Dim> | 
| 437 |  | inline Real dot( const Vector<Real, Dim>& v1, const Vector<Real, Dim>& v2 ) { | 
| 438 | < | Real tmp; | 
| 439 | < | tmp = 0; | 
| 438 | > | Real tmp; | 
| 439 | > | tmp = 0; | 
| 440 |  |  | 
| 441 | < | for (unsigned int i = 0; i < Dim; i++) | 
| 442 | < | tmp += v1[i] + v2[i]; | 
| 443 | < |  | 
| 444 | < | return tmp; | 
| 441 | > | for (unsigned int i = 0; i < Dim; i++) | 
| 442 | > | tmp += v1[i] * v2[i]; | 
| 443 | > |  | 
| 444 | > | return tmp; | 
| 445 |  | } | 
| 446 |  |  | 
| 447 |  | /** | 
| 472 |  | * Write to an output stream | 
| 473 |  | */ | 
| 474 |  | template<typename Real, unsigned int Dim> | 
| 475 | < | std::ostream &operator<< ( std::ostream& o, const Vector<Real, Dim>& v1 ) { | 
| 475 | > | std::ostream &operator<< ( std::ostream& o, const Vector<Real, Dim>& v) { | 
| 476 | > |  | 
| 477 | > | o << "[ "; | 
| 478 |  |  | 
| 479 | + | for (unsigned int i = 0 ; i< Dim; i++) { | 
| 480 | + | o << v[i]; | 
| 481 | + |  | 
| 482 | + | if (i  != Dim -1) { | 
| 483 | + | o<< ", "; | 
| 484 | + | } | 
| 485 | + | } | 
| 486 | + |  | 
| 487 | + | o << " ]"; | 
| 488 |  | return o; | 
| 489 |  | } | 
| 490 |  |  |