| 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, 24107 (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 | * | 
| 98 | *  This constructor fills the %DynamicVector with @a n copies of @a value. | 
| 99 | */ | 
| 100 | DynamicVector(size_type n, const value_type& value, | 
| 101 | const allocator_type& alloc = allocator_type()) | 
| 102 | : std::vector<Real, Alloc>(n, value, alloc){ } | 
| 103 |  | 
| 104 | /** | 
| 105 | *  @brief  Create a %DynamicVector with default elements. | 
| 106 | *  @param  n  The number of elements to initially create. | 
| 107 | * | 
| 108 | *  This constructor fills the %DynamicVector with @a n copies of a | 
| 109 | *  default-constructed element. | 
| 110 | */ | 
| 111 | explicit | 
| 112 | DynamicVector(size_type n) : std::vector<Real, Alloc>(n) { } | 
| 113 |  | 
| 114 | /** | 
| 115 | *  @brief  %Vector copy constructor. | 
| 116 | *  @param  x  A %DynamicVector of identical element and allocator types. | 
| 117 | * | 
| 118 | *  The newly-created %DynamicVector uses a copy of the allocation | 
| 119 | *  object used by @a x.  All the elements of @a x are copied, | 
| 120 | *  but any extra memory in | 
| 121 | *  @a x (for fast expansion) will not be copied. | 
| 122 | */ | 
| 123 | DynamicVector(const DynamicVector& x) | 
| 124 | : std::vector<Real, Alloc>(x) {} | 
| 125 |  | 
| 126 | template<typename _InputIterator> | 
| 127 | DynamicVector(_InputIterator first, _InputIterator last, | 
| 128 | const allocator_type& alloc = allocator_type()) | 
| 129 | : std::vector<Real, Alloc>(first, last, alloc) {} | 
| 130 |  | 
| 131 | inline Real operator()(unsigned int i) const{ | 
| 132 | return (*this)[i]; | 
| 133 | } | 
| 134 |  | 
| 135 | inline Real& operator()(unsigned int i){ | 
| 136 | return (*this)[i]; | 
| 137 | } | 
| 138 | /** | 
| 139 | * Tests if this vetor is equal to other vector | 
| 140 | * @return true if equal, otherwise return false | 
| 141 | * @param v vector to be compared | 
| 142 | */ | 
| 143 | inline bool operator ==(const DynamicVector<Real>& v) { | 
| 144 |  | 
| 145 | for (unsigned int i = 0; i < this->size(); i ++) { | 
| 146 | if (!equal((*this)[i], v[i])) { | 
| 147 | return false; | 
| 148 | } | 
| 149 | } | 
| 150 |  | 
| 151 | return true; | 
| 152 | } | 
| 153 |  | 
| 154 | /** | 
| 155 | * Tests if this vetor is not equal to other vector | 
| 156 | * @return true if equal, otherwise return false | 
| 157 | * @param v vector to be compared | 
| 158 | */ | 
| 159 | inline bool operator !=(const DynamicVector<Real>& v) { | 
| 160 | return !(*this == v); | 
| 161 | } | 
| 162 |  | 
| 163 | /** Negates the value of this vector in place. */ | 
| 164 | inline void negate() { | 
| 165 | for (unsigned int i = 0; i < this->size(); i++) | 
| 166 | (*this)[i] = -(*this)[i]; | 
| 167 | } | 
| 168 |  | 
| 169 | /** | 
| 170 | * Sets the value of this vector to the negation of vector v1. | 
| 171 | * @param v1 the source vector | 
| 172 | */ | 
| 173 | inline void negate(const DynamicVector<Real>& v1) { | 
| 174 | for (unsigned int i = 0; i < this->size(); i++) | 
| 175 | (*this)[i] = -v1[i]; | 
| 176 |  | 
| 177 | } | 
| 178 |  | 
| 179 | /** | 
| 180 | * Sets the value of this vector to the sum of itself and v1 (*this += v1). | 
| 181 | * @param v1 the other vector | 
| 182 | */ | 
| 183 | inline void add( const DynamicVector<Real>& v1 ) { | 
| 184 | std::transform(this->begin(), this->end(), v1.begin(),this->begin(),std::plus<Real>()); | 
| 185 | } | 
| 186 |  | 
| 187 | /** | 
| 188 | * Sets the value of this vector to the sum of v1 and v2 (*this = v1 + v2). | 
| 189 | * @param v1 the first vector | 
| 190 | * @param v2 the second vector | 
| 191 | */ | 
| 192 | inline void add( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 193 | std::transform(v1.begin(), v1.end(), v2.begin(),this->begin(),std::plus<Real>()); | 
| 194 | } | 
| 195 |  | 
| 196 | /** | 
| 197 | * Sets the value of this vector to the difference  of itself and v1 (*this -= v1). | 
| 198 | * @param v1 the other vector | 
| 199 | */ | 
| 200 | inline void sub( const DynamicVector<Real>& v1 ) { | 
| 201 | std::transform(this->begin(), this->end(), v1.begin(),this->begin(),std::minus<Real>()); | 
| 202 | } | 
| 203 |  | 
| 204 | /** | 
| 205 | * Sets the value of this vector to the difference of vector v1 and v2 (*this = v1 - v2). | 
| 206 | * @param v1 the first vector | 
| 207 | * @param v2 the second vector | 
| 208 | */ | 
| 209 | inline void sub( const DynamicVector<Real>& v1, const DynamicVector  &v2 ){ | 
| 210 | std::transform(v1.begin(), v1.end(), v2.begin(),this->begin(),std::minus<Real>()); | 
| 211 | } | 
| 212 |  | 
| 213 | /** | 
| 214 | * Sets the value of this vector to the scalar multiplication of itself (*this *= s). | 
| 215 | * @param s the scalar value | 
| 216 | */ | 
| 217 | inline void mul( Real s ) { | 
| 218 | for (unsigned int i = 0; i < this->size(); i++) | 
| 219 | (*this)[i] *= s; | 
| 220 | } | 
| 221 |  | 
| 222 | /** | 
| 223 | * Sets the value of this vector to the scalar multiplication of vector v1 | 
| 224 | * (*this = s * v1). | 
| 225 | * @param v1 the vector | 
| 226 | * @param s the scalar value | 
| 227 | */ | 
| 228 | inline void mul( const DynamicVector<Real>& v1, Real s) { | 
| 229 | this->resize(v1.size()); | 
| 230 | for (unsigned int i = 0; i < this->size(); i++) | 
| 231 | (*this)[i] = s * v1[i]; | 
| 232 | } | 
| 233 |  | 
| 234 | /** | 
| 235 | * Sets the value of this vector to the scalar division of itself  (*this /= s ). | 
| 236 | * @param s the scalar value | 
| 237 | */ | 
| 238 | inline void div( Real s) { | 
| 239 | for (unsigned int i = 0; i < this->size(); i++) | 
| 240 | (*this)[i] /= s; | 
| 241 | } | 
| 242 |  | 
| 243 | /** | 
| 244 | * Sets the value of this vector to the scalar division of vector v1  (*this = v1 / s ). | 
| 245 | * @param v1 the source vector | 
| 246 | * @param s the scalar value | 
| 247 | */ | 
| 248 | inline void div( const DynamicVector<Real>& v1, Real s ) { | 
| 249 | for (unsigned int i = 0; i < this->size(); i++) | 
| 250 | (*this)[i] = v1[i] / s; | 
| 251 | } | 
| 252 |  | 
| 253 | /** @see #add */ | 
| 254 | inline DynamicVector<Real>& operator +=( const DynamicVector<Real>& v1 ) { | 
| 255 | add(v1); | 
| 256 | return *this; | 
| 257 | } | 
| 258 |  | 
| 259 | /** @see #sub */ | 
| 260 | inline DynamicVector<Real>& operator -=( const DynamicVector<Real>& v1 ) { | 
| 261 | sub(v1); | 
| 262 | return *this; | 
| 263 | } | 
| 264 |  | 
| 265 | /** @see #mul */ | 
| 266 | inline DynamicVector<Real>& operator *=( Real s) { | 
| 267 | mul(s); | 
| 268 | return *this; | 
| 269 | } | 
| 270 |  | 
| 271 | /** @see #div */ | 
| 272 | inline DynamicVector<Real>& operator /=( Real s ) { | 
| 273 | div(s); | 
| 274 | return *this; | 
| 275 | } | 
| 276 |  | 
| 277 | /** zero out the vector */ | 
| 278 | inline void setZero( ) { | 
| 279 | for (unsigned int i = 0; i < this->size(); i++) | 
| 280 | (*this)[i] = 0; | 
| 281 | } | 
| 282 |  | 
| 283 | /** | 
| 284 | * Returns the length of this vector. | 
| 285 | * @return the length of this vector | 
| 286 | */ | 
| 287 | inline Real length() { | 
| 288 | return sqrt(lengthSquare()); | 
| 289 | } | 
| 290 |  | 
| 291 | /** | 
| 292 | * Returns the squared length of this vector. | 
| 293 | * @return the squared length of this vector | 
| 294 | */ | 
| 295 | inline Real lengthSquare() { | 
| 296 | return dot(*this, *this); | 
| 297 | } | 
| 298 |  | 
| 299 | /** Normalizes this vector in place */ | 
| 300 | inline void normalize() { | 
| 301 | Real len; | 
| 302 |  | 
| 303 | len = length(); | 
| 304 |  | 
| 305 | //if (len < OpenMD::NumericConstant::epsilon) | 
| 306 | //  throw(); | 
| 307 |  | 
| 308 | *this /= len; | 
| 309 | } | 
| 310 |  | 
| 311 | /** | 
| 312 | * Tests if this vector is normalized | 
| 313 | * @return true if this vector is normalized, otherwise return false | 
| 314 | */ | 
| 315 | inline bool isNormalized() { | 
| 316 | return equal(lengthSquare(), 1.0); | 
| 317 | } | 
| 318 |  | 
| 319 | template<class VectorType> | 
| 320 | void getSubVector(unsigned int beginning, VectorType& v) { | 
| 321 | assert(beginning + v.size() -1 <= this->size()); | 
| 322 |  | 
| 323 | for (unsigned int i = 0; i < v.size(); ++i) | 
| 324 | v(i) = (*this)[beginning+i]; | 
| 325 | } | 
| 326 |  | 
| 327 |  | 
| 328 | }; | 
| 329 |  | 
| 330 | /** unary minus*/ | 
| 331 | template<typename Real> | 
| 332 | inline DynamicVector<Real> operator -(const DynamicVector<Real>& v1){ | 
| 333 | DynamicVector<Real> tmp(v1); | 
| 334 | tmp.negate(); | 
| 335 | return tmp; | 
| 336 | } | 
| 337 |  | 
| 338 | /** | 
| 339 | * Return the sum of two vectors  (v1 - v2). | 
| 340 | * @return the sum of two vectors | 
| 341 | * @param v1 the first vector | 
| 342 | * @param v2 the second vector | 
| 343 | */ | 
| 344 | template<typename Real> | 
| 345 | inline DynamicVector<Real> operator +(const DynamicVector<Real>& v1, const DynamicVector<Real>& v2) { | 
| 346 | assert(v1.size() == v2.size()); | 
| 347 | DynamicVector<Real>result(v1.size()); | 
| 348 | result.add(v1, v2); | 
| 349 | return result; | 
| 350 | } | 
| 351 |  | 
| 352 | /** | 
| 353 | * Return the difference of two vectors  (v1 - v2). | 
| 354 | * @return the difference of two vectors | 
| 355 | * @param v1 the first vector | 
| 356 | * @param v2 the second vector | 
| 357 | */ | 
| 358 | template<typename Real> | 
| 359 | DynamicVector<Real> operator -(const DynamicVector<Real>& v1, const DynamicVector<Real>& v2) { | 
| 360 | assert(v1.size() == v2.size()); | 
| 361 | DynamicVector<Real> result(v1.size()); | 
| 362 | result.sub(v1, v2); | 
| 363 | return result; | 
| 364 | } | 
| 365 |  | 
| 366 | /** | 
| 367 | * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). | 
| 368 | * @return  the vaule of scalar multiplication of this vector | 
| 369 | * @param v1 the source vector | 
| 370 | * @param s the scalar value | 
| 371 | */ | 
| 372 | template<typename Real> | 
| 373 | DynamicVector<Real> operator *( const DynamicVector<Real>& v1, Real s) { | 
| 374 | DynamicVector<Real> result(v1.size()); | 
| 375 | result.mul(v1,s); | 
| 376 | return result; | 
| 377 | } | 
| 378 |  | 
| 379 | /** | 
| 380 | * Returns the vaule of scalar multiplication of this vector v1 (v1 * r). | 
| 381 | * @return  the vaule of scalar multiplication of this vector | 
| 382 | * @param s the scalar value | 
| 383 | * @param v1 the source vector | 
| 384 | */ | 
| 385 | template<typename Real> | 
| 386 | DynamicVector<Real> operator *( Real s, const DynamicVector<Real>& v1 ) { | 
| 387 | DynamicVector<Real> result(v1.size()); | 
| 388 | result.mul(v1, s); | 
| 389 | return result; | 
| 390 | } | 
| 391 |  | 
| 392 | /** | 
| 393 | * Returns the  value of division of a vector by a scalar. | 
| 394 | * @return  the vaule of scalar division of this vector | 
| 395 | * @param v1 the source vector | 
| 396 | * @param s the scalar value | 
| 397 | */ | 
| 398 | template<typename Real> | 
| 399 | DynamicVector<Real> operator / ( const DynamicVector<Real>& v1, Real s) { | 
| 400 | DynamicVector<Real> result(v1.size()); | 
| 401 | result.div( v1,s); | 
| 402 | return result; | 
| 403 | } | 
| 404 |  | 
| 405 | /** | 
| 406 | * Returns the dot product of two DynamicVectors | 
| 407 | * @param v1 first vector | 
| 408 | * @param v2 second vector | 
| 409 | * @return the dot product of v1 and v2 | 
| 410 | */ | 
| 411 | template<typename Real> | 
| 412 | inline Real dot( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 413 | Real tmp; | 
| 414 | tmp = 0; | 
| 415 | assert(v1.size() == v2.size()); | 
| 416 | for (unsigned int i = 0; i < v1.size(); i++) | 
| 417 | tmp += v1[i] * v2[i]; | 
| 418 |  | 
| 419 | return tmp; | 
| 420 | } | 
| 421 |  | 
| 422 | /** | 
| 423 | * Returns the distance between  two DynamicVectors | 
| 424 | * @param v1 first vector | 
| 425 | * @param v2 second vector | 
| 426 | * @return the distance between v1 and v2 | 
| 427 | */ | 
| 428 | template<typename Real> | 
| 429 | inline Real distance( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 430 | DynamicVector<Real> tempDynamicVector = v1 - v2; | 
| 431 | return tempDynamicVector.length(); | 
| 432 | } | 
| 433 |  | 
| 434 | /** | 
| 435 | * Returns the squared distance between  two DynamicVectors | 
| 436 | * @param v1 first vector | 
| 437 | * @param v2 second vector | 
| 438 | * @return the squared distance between v1 and v2 | 
| 439 | */ | 
| 440 | template<typename Real> | 
| 441 | inline Real distanceSquare( const DynamicVector<Real>& v1, const DynamicVector<Real>& v2 ) { | 
| 442 | DynamicVector<Real> tempDynamicVector = v1 - v2; | 
| 443 | return tempDynamicVector.lengthSquare(); | 
| 444 | } | 
| 445 |  | 
| 446 | /** | 
| 447 | * Write to an output stream | 
| 448 | */ | 
| 449 | template<typename Real> | 
| 450 | std::ostream &operator<< ( std::ostream& o, const DynamicVector<Real>& v) { | 
| 451 |  | 
| 452 | o << "[ "; | 
| 453 |  | 
| 454 | for (unsigned int i = 0 ; i< v.size(); i++) { | 
| 455 | o << v[i]; | 
| 456 |  | 
| 457 | if (i  != v.size() -1) { | 
| 458 | o<< ", "; | 
| 459 | } | 
| 460 | } | 
| 461 |  | 
| 462 | o << " ]"; | 
| 463 | return o; | 
| 464 | } | 
| 465 |  | 
| 466 | } | 
| 467 | #endif | 
| 468 |  |