| 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. Redistributions of source code must retain the above copyright | 
| 10 | *    notice, this list of conditions and the following disclaimer. | 
| 11 | * | 
| 12 | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 | *    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 | * | 
| 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, 234107 (2008). | 
| 39 | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 | */ | 
| 42 |  | 
| 43 | /** | 
| 44 | * @file DynamicVector.hpp | 
| 45 | * @author Teng Lin | 
| 46 | * @date 09/14/2004 | 
| 47 | * @version 1.0 | 
| 48 | */ | 
| 49 |  | 
| 50 | #ifndef MATH_DYNAMICVECTOR_HPP | 
| 51 | #define MATH_DYNAMICVECTOR_HPP | 
| 52 |  | 
| 53 | #include <cassert> | 
| 54 | #include <cmath> | 
| 55 | #include <iostream> | 
| 56 | #include <math.h> | 
| 57 | #include <algorithm> | 
| 58 | #include <vector> | 
| 59 |  | 
| 60 | namespace OpenMD { | 
| 61 |  | 
| 62 | /** | 
| 63 | * @class DynamicVector DynamicVector.hpp "math/DynamicVector.hpp" | 
| 64 | * @brief Fix length vector class | 
| 65 | */ | 
| 66 | template<typename Real, typename Alloc = std::allocator<Real> > | 
| 67 | class DynamicVector : public std::vector<Real, Alloc> { | 
| 68 |  | 
| 69 | public: | 
| 70 | typedef Real value_type; | 
| 71 | typedef std::vector<Real, Alloc> VectorType; | 
| 72 | typedef typename VectorType::pointer pointer; | 
| 73 | typedef typename VectorType::const_pointer const_pointer; | 
| 74 | typedef typename VectorType::reference reference; | 
| 75 | typedef typename VectorType::const_reference const_reference; | 
| 76 | typedef typename VectorType::iterator iterator; | 
| 77 | typedef typename VectorType::const_iterator const_iterator; | 
| 78 | typedef typename VectorType::const_reverse_iterator  const_reverse_iterator; | 
| 79 | typedef typename VectorType::reverse_iterator reverse_iterator; | 
| 80 | typedef typename VectorType::size_type size_type; | 
| 81 | typedef typename VectorType::difference_type difference_type; | 
| 82 | typedef typename VectorType::allocator_type allocator_type; | 
| 83 |  | 
| 84 |  | 
| 85 | // [23.2.4.1] construct/copy/destroy | 
| 86 | // (assign() and get_allocator() are also listed in this section) | 
| 87 | /** | 
| 88 | *  @brief  Default constructor creates no elements. | 
| 89 | */      explicit | 
| 90 | DynamicVector(const allocator_type& alloc = allocator_type()) | 
| 91 | : std::vector<Real, Alloc>(alloc) { } | 
| 92 |  | 
| 93 | /** | 
| 94 | *  @brief  Create a %DynamicVector with copies of an exemplar element. | 
| 95 | *  @param  n  The number of elements to initially create. | 
| 96 | *  @param  value  An element to copy. | 
| 97 | *  @param  alloc  The allocator_type to use | 
| 98 | * | 
| 99 | *  This constructor fills the %DynamicVector with @a n copies of @a value. | 
| 100 | */ | 
| 101 | DynamicVector(size_type n, const value_type& value, | 
| 102 | const allocator_type& alloc = allocator_type()) | 
| 103 | : std::vector<Real, Alloc>(n, value, alloc){ } | 
| 104 |  | 
| 105 | /** | 
| 106 | *  @brief  Create a %DynamicVector with default elements. | 
| 107 | *  @param  n  The number of elements to initially create. | 
| 108 | * | 
| 109 | *  This constructor fills the %DynamicVector with @a n copies of a | 
| 110 | *  default-constructed element. | 
| 111 | */ | 
| 112 | explicit | 
| 113 | DynamicVector(size_type n) : std::vector<Real, Alloc>(n) { } | 
| 114 |  | 
| 115 | /** | 
| 116 | *  @brief  %Vector copy constructor. | 
| 117 | *  @param  x  A %DynamicVector of identical element and allocator types. | 
| 118 | * | 
| 119 | *  The newly-created %DynamicVector uses a copy of the allocation | 
| 120 | *  object used by @a x.  All the elements of @a x are copied, | 
| 121 | *  but any extra memory in | 
| 122 | *  @a x (for fast expansion) will not be copied. | 
| 123 | */ | 
| 124 | DynamicVector(const DynamicVector& x) | 
| 125 | : std::vector<Real, Alloc>(x) {} | 
| 126 |  | 
| 127 | template<typename _InputIterator> | 
| 128 | DynamicVector(_InputIterator first, _InputIterator last, | 
| 129 | const allocator_type& alloc = allocator_type()) | 
| 130 | : std::vector<Real, Alloc>(first, last, alloc) {} | 
| 131 |  | 
| 132 | inline Real operator()(unsigned int i) const{ | 
| 133 | return (*this)[i]; | 
| 134 | } | 
| 135 |  | 
| 136 | inline Real& operator()(unsigned int i){ | 
| 137 | return (*this)[i]; | 
| 138 | } | 
| 139 | /** | 
| 140 | * Tests if this vetor is equal to other vector | 
| 141 | * @return true if equal, otherwise return false | 
| 142 | * @param v vector to be compared | 
| 143 | */ | 
| 144 | inline bool operator ==(const DynamicVector<Real>& v) { | 
| 145 |  | 
| 146 | for (unsigned int i = 0; i < this->size(); i ++) { | 
| 147 | if (!equal((*this)[i], v[i])) { | 
| 148 | return false; | 
| 149 | } | 
| 150 | } | 
| 151 |  | 
| 152 | return true; | 
| 153 | } | 
| 154 |  | 
| 155 | /** | 
| 156 | * Tests if this vetor is not equal to other vector | 
| 157 | * @return true if equal, otherwise return false | 
| 158 | * @param v vector to be compared | 
| 159 | */ | 
| 160 | inline bool operator !=(const DynamicVector<Real>& v) { | 
| 161 | return !(*this == v); | 
| 162 | } | 
| 163 |  | 
| 164 | /** Negates the value of this vector in place. */ | 
| 165 | inline void negate() { | 
| 166 | for (unsigned int i = 0; i < this->size(); i++) | 
| 167 | (*this)[i] = -(*this)[i]; | 
| 168 | } | 
| 169 |  | 
| 170 | /** | 
| 171 | * Sets the value of this vector to the negation of vector v1. | 
| 172 | * @param v1 the source vector | 
| 173 | */ | 
| 174 | inline void negate(const DynamicVector<Real>& v1) { | 
| 175 | for (unsigned int i = 0; i < this->size(); i++) | 
| 176 | (*this)[i] = -v1[i]; | 
| 177 |  | 
| 178 | } | 
| 179 |  | 
| 180 | /** | 
| 181 | * Sets the value of this vector to the sum of itself and v1 (*this += v1). | 
| 182 | * @param v1 the other vector | 
| 183 | */ | 
| 184 | inline void add( const DynamicVector<Real>& v1 ) { | 
| 185 | std::transform(this->begin(), this->end(), v1.begin(),this->begin(),std::plus<Real>()); | 
| 186 | } | 
| 187 |  | 
| 188 | /** | 
| 189 | * Sets the value of this vector to the sum of v1 and v2 (*this = v1 + v2). | 
| 190 | * @param v1 the first vector | 
| 191 | * @param v2 the second vector | 
| 192 | */ | 
| 193 | inline void add( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 194 | std::transform(v1.begin(), v1.end(), v2.begin(),this->begin(),std::plus<Real>()); | 
| 195 | } | 
| 196 |  | 
| 197 | /** | 
| 198 | * Sets the value of this vector to the difference  of itself and v1 (*this -= v1). | 
| 199 | * @param v1 the other vector | 
| 200 | */ | 
| 201 | inline void sub( const DynamicVector<Real>& v1 ) { | 
| 202 | std::transform(this->begin(), this->end(), v1.begin(),this->begin(),std::minus<Real>()); | 
| 203 | } | 
| 204 |  | 
| 205 | /** | 
| 206 | * Sets the value of this vector to the difference of vector v1 and v2 (*this = v1 - v2). | 
| 207 | * @param v1 the first vector | 
| 208 | * @param v2 the second vector | 
| 209 | */ | 
| 210 | inline void sub( const DynamicVector<Real>& v1, const DynamicVector  &v2 ){ | 
| 211 | std::transform(v1.begin(), v1.end(), v2.begin(),this->begin(),std::minus<Real>()); | 
| 212 | } | 
| 213 |  | 
| 214 | /** | 
| 215 | * Sets the value of this vector to the scalar multiplication of itself (*this *= s). | 
| 216 | * @param s the scalar value | 
| 217 | */ | 
| 218 | inline void mul( Real s ) { | 
| 219 | for (unsigned int i = 0; i < this->size(); i++) | 
| 220 | (*this)[i] *= s; | 
| 221 | } | 
| 222 |  | 
| 223 | /** | 
| 224 | * Sets the value of this vector to the scalar multiplication of vector v1 | 
| 225 | * (*this = s * v1). | 
| 226 | * @param v1 the vector | 
| 227 | * @param s the scalar value | 
| 228 | */ | 
| 229 | inline void mul( const DynamicVector<Real>& v1, Real s) { | 
| 230 | this->resize(v1.size()); | 
| 231 | for (unsigned int i = 0; i < this->size(); i++) | 
| 232 | (*this)[i] = s * v1[i]; | 
| 233 | } | 
| 234 |  | 
| 235 | /** | 
| 236 | * Sets the value of this vector to the scalar division of itself  (*this /= s ). | 
| 237 | * @param s the scalar value | 
| 238 | */ | 
| 239 | inline void div( Real s) { | 
| 240 | for (unsigned int i = 0; i < this->size(); i++) | 
| 241 | (*this)[i] /= s; | 
| 242 | } | 
| 243 |  | 
| 244 | /** | 
| 245 | * Sets the value of this vector to the scalar division of vector v1  (*this = v1 / s ). | 
| 246 | * @param v1 the source vector | 
| 247 | * @param s the scalar value | 
| 248 | */ | 
| 249 | inline void div( const DynamicVector<Real>& v1, Real s ) { | 
| 250 | for (unsigned int i = 0; i < this->size(); i++) | 
| 251 | (*this)[i] = v1[i] / s; | 
| 252 | } | 
| 253 |  | 
| 254 | /** @see #add */ | 
| 255 | inline DynamicVector<Real>& operator +=( const DynamicVector<Real>& v1 ) { | 
| 256 | add(v1); | 
| 257 | return *this; | 
| 258 | } | 
| 259 |  | 
| 260 | /** @see #sub */ | 
| 261 | inline DynamicVector<Real>& operator -=( const DynamicVector<Real>& v1 ) { | 
| 262 | sub(v1); | 
| 263 | return *this; | 
| 264 | } | 
| 265 |  | 
| 266 | /** @see #mul */ | 
| 267 | inline DynamicVector<Real>& operator *=( Real s) { | 
| 268 | mul(s); | 
| 269 | return *this; | 
| 270 | } | 
| 271 |  | 
| 272 | /** @see #div */ | 
| 273 | inline DynamicVector<Real>& operator /=( Real s ) { | 
| 274 | div(s); | 
| 275 | return *this; | 
| 276 | } | 
| 277 |  | 
| 278 | /** zero out the vector */ | 
| 279 | inline void setZero( ) { | 
| 280 | for (unsigned int i = 0; i < this->size(); i++) | 
| 281 | (*this)[i] = 0; | 
| 282 | } | 
| 283 |  | 
| 284 | /** | 
| 285 | * Returns the length of this vector. | 
| 286 | * @return the length of this vector | 
| 287 | */ | 
| 288 | inline Real length() { | 
| 289 | return sqrt(lengthSquare()); | 
| 290 | } | 
| 291 |  | 
| 292 | /** | 
| 293 | * Returns the squared length of this vector. | 
| 294 | * @return the squared length of this vector | 
| 295 | */ | 
| 296 | inline Real lengthSquare() { | 
| 297 | return dot(*this, *this); | 
| 298 | } | 
| 299 |  | 
| 300 | /** Normalizes this vector in place */ | 
| 301 | inline void normalize() { | 
| 302 | Real len; | 
| 303 |  | 
| 304 | len = length(); | 
| 305 |  | 
| 306 | //if (len < OpenMD::NumericConstant::epsilon) | 
| 307 | //  throw(); | 
| 308 |  | 
| 309 | *this /= len; | 
| 310 | } | 
| 311 |  | 
| 312 | /** | 
| 313 | * Tests if this vector is normalized | 
| 314 | * @return true if this vector is normalized, otherwise return false | 
| 315 | */ | 
| 316 | inline bool isNormalized() { | 
| 317 | return equal(lengthSquare(), 1.0); | 
| 318 | } | 
| 319 |  | 
| 320 | template<class VectorType> | 
| 321 | void getSubVector(unsigned int beginning, VectorType& v) { | 
| 322 | assert(beginning + v.size() -1 <= this->size()); | 
| 323 |  | 
| 324 | for (unsigned int i = 0; i < v.size(); ++i) | 
| 325 | v(i) = (*this)[beginning+i]; | 
| 326 | } | 
| 327 |  | 
| 328 |  | 
| 329 | }; | 
| 330 |  | 
| 331 | /** unary minus*/ | 
| 332 | template<typename Real> | 
| 333 | inline DynamicVector<Real> operator -(const DynamicVector<Real>& v1){ | 
| 334 | DynamicVector<Real> tmp(v1); | 
| 335 | tmp.negate(); | 
| 336 | return tmp; | 
| 337 | } | 
| 338 |  | 
| 339 | /** | 
| 340 | * Return the sum of two vectors  (v1 - v2). | 
| 341 | * @return the sum of two vectors | 
| 342 | * @param v1 the first vector | 
| 343 | * @param v2 the second vector | 
| 344 | */ | 
| 345 | template<typename Real> | 
| 346 | inline DynamicVector<Real> operator +(const DynamicVector<Real>& v1, const DynamicVector<Real>& v2) { | 
| 347 | assert(v1.size() == v2.size()); | 
| 348 | DynamicVector<Real>result(v1.size()); | 
| 349 | result.add(v1, v2); | 
| 350 | return result; | 
| 351 | } | 
| 352 |  | 
| 353 | /** | 
| 354 | * Return the difference of two vectors  (v1 - v2). | 
| 355 | * @return the difference of two vectors | 
| 356 | * @param v1 the first vector | 
| 357 | * @param v2 the second vector | 
| 358 | */ | 
| 359 | template<typename Real> | 
| 360 | DynamicVector<Real> operator -(const DynamicVector<Real>& v1, const DynamicVector<Real>& v2) { | 
| 361 | assert(v1.size() == v2.size()); | 
| 362 | DynamicVector<Real> result(v1.size()); | 
| 363 | result.sub(v1, v2); | 
| 364 | return result; | 
| 365 | } | 
| 366 |  | 
| 367 | /** | 
| 368 | * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). | 
| 369 | * @return  the vaule of scalar multiplication of this vector | 
| 370 | * @param v1 the source vector | 
| 371 | * @param s the scalar value | 
| 372 | */ | 
| 373 | template<typename Real> | 
| 374 | DynamicVector<Real> operator *( const DynamicVector<Real>& v1, Real s) { | 
| 375 | DynamicVector<Real> result(v1.size()); | 
| 376 | result.mul(v1,s); | 
| 377 | return result; | 
| 378 | } | 
| 379 |  | 
| 380 | /** | 
| 381 | * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). | 
| 382 | * @return  the vaule of scalar multiplication of this vector | 
| 383 | * @param s the scalar value | 
| 384 | * @param v1 the source vector | 
| 385 | */ | 
| 386 | template<typename Real> | 
| 387 | DynamicVector<Real> operator *( Real s, const DynamicVector<Real>& v1 ) { | 
| 388 | DynamicVector<Real> result(v1.size()); | 
| 389 | result.mul(v1, s); | 
| 390 | return result; | 
| 391 | } | 
| 392 |  | 
| 393 | /** | 
| 394 | * Returns the  value of division of a vector by a scalar. | 
| 395 | * @return  the vaule of scalar division of this vector | 
| 396 | * @param v1 the source vector | 
| 397 | * @param s the scalar value | 
| 398 | */ | 
| 399 | template<typename Real> | 
| 400 | DynamicVector<Real> operator / ( const DynamicVector<Real>& v1, Real s) { | 
| 401 | DynamicVector<Real> result(v1.size()); | 
| 402 | result.div( v1,s); | 
| 403 | return result; | 
| 404 | } | 
| 405 |  | 
| 406 | /** | 
| 407 | * Returns the dot product of two DynamicVectors | 
| 408 | * @param v1 first vector | 
| 409 | * @param v2 second vector | 
| 410 | * @return the dot product of v1 and v2 | 
| 411 | */ | 
| 412 | template<typename Real> | 
| 413 | inline Real dot( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 414 | Real tmp; | 
| 415 | tmp = 0; | 
| 416 | assert(v1.size() == v2.size()); | 
| 417 | for (unsigned int i = 0; i < v1.size(); i++) | 
| 418 | tmp += v1[i] * v2[i]; | 
| 419 |  | 
| 420 | return tmp; | 
| 421 | } | 
| 422 |  | 
| 423 | /** | 
| 424 | * Returns the distance between  two DynamicVectors | 
| 425 | * @param v1 first vector | 
| 426 | * @param v2 second vector | 
| 427 | * @return the distance between v1 and v2 | 
| 428 | */ | 
| 429 | template<typename Real> | 
| 430 | inline Real distance( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 431 | DynamicVector<Real> tempDynamicVector = v1 - v2; | 
| 432 | return tempDynamicVector.length(); | 
| 433 | } | 
| 434 |  | 
| 435 | /** | 
| 436 | * Returns the squared distance between  two DynamicVectors | 
| 437 | * @param v1 first vector | 
| 438 | * @param v2 second vector | 
| 439 | * @return the squared distance between v1 and v2 | 
| 440 | */ | 
| 441 | template<typename Real> | 
| 442 | inline Real distanceSquare( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 443 | DynamicVector<Real> tempDynamicVector = v1 - v2; | 
| 444 | return tempDynamicVector.lengthSquare(); | 
| 445 | } | 
| 446 |  | 
| 447 | /** | 
| 448 | * Write to an output stream | 
| 449 | */ | 
| 450 | template<typename Real> | 
| 451 | std::ostream &operator<< ( std::ostream& o, const DynamicVector<Real>& v) { | 
| 452 |  | 
| 453 | o << "[ "; | 
| 454 |  | 
| 455 | for (unsigned int i = 0 ; i< v.size(); i++) { | 
| 456 | o << v[i]; | 
| 457 |  | 
| 458 | if (i  != v.size() -1) { | 
| 459 | o<< ", "; | 
| 460 | } | 
| 461 | } | 
| 462 |  | 
| 463 | o << " ]"; | 
| 464 | return o; | 
| 465 | } | 
| 466 |  | 
| 467 | } | 
| 468 | #endif | 
| 469 |  |