src/gromacs/math/vectypes.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2016,2017,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifndef GMX_MATH_VECTYPES_H
  39 #define GMX_MATH_VECTYPES_H
  40
  41 #include <cassert>
  42 #include <cmath>
  43
  44 #include <algorithm>
  45 #include <type_traits>
  46
  47 #include "gromacs/utility/real.h"
  48
  49 #define XX 0 /* Defines for indexing in */
  50 #define YY 1 /* vectors                 */
  51 #define ZZ 2
  52 #define DIM 3 /* Dimension of vectors    */
  53
  54 typedef real rvec[DIM];
  55
  56 typedef double dvec[DIM];
  57
  58 typedef real matrix[DIM][DIM];
  59
  60 typedef real tensor[DIM][DIM];
  61
  62 typedef int ivec[DIM];
  63
  64 namespace gmx
  65 {
  66
  67 /*! \brief
  68  * C++ class for 3D vectors.
  69  *
  70  * \tparam ValueType  Type
  71  *
  72  * This class provides a C++ version of rvec/dvec/ivec that can be put into STL
  73  * containers etc.  It is more or less a drop-in replacement for `rvec` and
  74  * friends: it can be used in most contexts that accept the equivalent C type.
  75  * However, there is one case where explicit conversion is necessary:
  76  *  - An array of these objects needs to be converted with as_vec_array() (or
  77  *    convenience methods like as_rvec_array()).
  78  *
  79  * For the array conversion to work, the compiler should not add any extra
  80  * alignment/padding in the layout of this class;  that this actually works as
  81  * intended is tested in the unit tests.
  82  *
  83  * \inpublicapi
  84  */
  85 template<typename ValueType>
  86 class BasicVector
  87 {
  88 public:
  89     //! Underlying raw C array type (rvec/dvec/ivec).
  90     using RawArray = ValueType[DIM];
  91
  92     // The code here assumes ValueType has been deduced as a data type like int
  93     // and not a pointer like int*. If there is a use case for a 3-element array
  94     // of pointers, the implementation will be different enough that the whole
  95     // template class should have a separate partial specialization. We try to avoid
  96     // accidental matching to pointers, but this assertion is a no-cost extra check.
  97     //
  98     // TODO: Use std::is_pointer_v when CUDA 11 is a requirement.
  99     static_assert(!std::is_pointer<std::remove_cv_t<ValueType>>::value,
 100                   "BasicVector value type must not be a pointer.");
 101
 102     //! Constructs default (uninitialized) vector.
 103     BasicVector() {}
 104     //! Constructs a vector from given values.
 105     BasicVector(ValueType x, ValueType y, ValueType z) : x_{ x, y, z } {}
 106     /*! \brief
 107      * Constructs a vector from given values.
 108      *
 109      * This constructor is not explicit to support implicit conversions
 110      * that allow, e.g., calling `std::vector<RVec>:``:push_back()` directly
 111      * with an `rvec` parameter.
 112      */
 113     BasicVector(const RawArray x) : x_{ x[XX], x[YY], x[ZZ] } {}
 114     //! Default copy constructor.
 115     BasicVector(const BasicVector& src) = default;
 116     //! Default copy assignment operator.
 117     BasicVector& operator=(const BasicVector& v) = default;
 118     //! Default move constructor.
 119     BasicVector(BasicVector&& src) noexcept = default;
 120     //! Default move assignment operator.
 121     BasicVector& operator=(BasicVector&& v) noexcept = default;
 122     //! Indexing operator to make the class work as the raw array.
 123     ValueType& operator[](int i) { return x_[i]; }
 124     //! Indexing operator to make the class work as the raw array.
 125     ValueType operator[](int i) const { return x_[i]; }
 126     //! Allow inplace addition for BasicVector
 127     BasicVector<ValueType>& operator+=(const BasicVector<ValueType>& right)
 128     {
 129         return *this = *this + right;
 130     }
 131     //! Allow inplace subtraction for BasicVector
 132     BasicVector<ValueType>& operator-=(const BasicVector<ValueType>& right)
 133     {
 134         return *this = *this - right;
 135     }
 136     //! Allow vector addition
 137     BasicVector<ValueType> operator+(const BasicVector<ValueType>& right) const
 138     {
 139         return { x_[0] + right[0], x_[1] + right[1], x_[2] + right[2] };
 140     }
 141     //! Allow vector subtraction
 142     BasicVector<ValueType> operator-(const BasicVector<ValueType>& right) const
 143     {
 144         return { x_[0] - right[0], x_[1] - right[1], x_[2] - right[2] };
 145     }
 146     //! Allow vector scalar division
 147     BasicVector<ValueType> operator/(const ValueType& right) const
 148     {
 149         assert((right != 0 && "Cannot divide by zero"));
 150
 151         return *this * (1 / right);
 152     }
 153     //! Scale vector by a scalar
 154     BasicVector<ValueType>& operator*=(const ValueType& right)
 155     {
 156         x_[0] *= right;
 157         x_[1] *= right;
 158         x_[2] *= right;
 159
 160         return *this;
 161     }
 162     //! Divide vector by a scalar
 163     BasicVector<ValueType>& operator/=(const ValueType& right)
 164     {
 165         assert((right != 0 && "Cannot divide by zero"));
 166
 167         return *this *= 1 / right;
 168     }
 169     //! Return dot product
 170     ValueType dot(const BasicVector<ValueType>& right) const
 171     {
 172         return x_[0] * right[0] + x_[1] * right[1] + x_[2] * right[2];
 173     }
 174
 175     //! Allow vector vector multiplication (cross product)
 176     BasicVector<ValueType> cross(const BasicVector<ValueType>& right) const
 177     {
 178         return { x_[YY] * right.x_[ZZ] - x_[ZZ] * right.x_[YY],
 179                  x_[ZZ] * right.x_[XX] - x_[XX] * right.x_[ZZ],
 180                  x_[XX] * right.x_[YY] - x_[YY] * right.x_[XX] };
 181     }
 182
 183     //! Return normalized to unit vector
 184     BasicVector<ValueType> unitVector() const
 185     {
 186         const ValueType vectorNorm = norm();
 187         assert((vectorNorm != 0 && "unitVector() should not be called with a zero vector"));
 188
 189         return *this / vectorNorm;
 190     }
 191
 192     //! Length^2 of vector
 193     ValueType norm2() const { return dot(*this); }
 194
 195     //! Norm or length of vector
 196     ValueType norm() const { return std::sqrt(norm2()); }
 197
 198     //! cast to RVec
 199     BasicVector<real> toRVec() const { return { real(x_[0]), real(x_[1]), real(x_[2]) }; }
 200
 201     //! cast to IVec
 202     BasicVector<int> toIVec() const
 203     {
 204         return { static_cast<int>(x_[0]), static_cast<int>(x_[1]), static_cast<int>(x_[2]) };
 205     }
 206
 207     //! cast to DVec
 208     BasicVector<double> toDVec() const { return { double(x_[0]), double(x_[1]), double(x_[2]) }; }
 209
 210     //! Converts to a raw C array where implicit conversion does not work.
 211     RawArray& as_vec() { return x_; }
 212     //! Converts to a raw C array where implicit conversion does not work.
 213     const RawArray& as_vec() const { return x_; }
 214     //! Makes BasicVector usable in contexts where a raw C array is expected.
 215     operator RawArray&() { return x_; }
 216     //! Makes BasicVector usable in contexts where a raw C array is expected.
 217     operator const RawArray&() const { return x_; }
 218
 219 private:
 220     RawArray x_;
 221 };
 222
 223 //! Allow vector scalar multiplication
 224 template<typename ValueType>
 225 BasicVector<ValueType> operator*(const BasicVector<ValueType>& basicVector, const ValueType& scalar)
 226 {
 227     return { basicVector[0] * scalar, basicVector[1] * scalar, basicVector[2] * scalar };
 228 }
 229
 230 //! Allow scalar vector multiplication
 231 template<typename ValueType>
 232 BasicVector<ValueType> operator*(const ValueType& scalar, const BasicVector<ValueType>& basicVector)
 233 {
 234     return { scalar * basicVector[0], scalar * basicVector[1], scalar * basicVector[2] };
 235 }
 236
 237 /*! \brief
 238  * unitv for gmx::BasicVector
 239  */
 240 template<typename VectorType>
 241 static inline VectorType unitVector(const VectorType& v)
 242 {
 243     return v.unitVector();
 244 }
 245
 246 /*! \brief
 247  * norm for gmx::BasicVector
 248  */
 249 template<typename ValueType>
 250 static inline ValueType norm(BasicVector<ValueType> v)
 251 {
 252     return v.norm();
 253 }
 254
 255 /*! \brief
 256  * Square of the vector norm for gmx::BasicVector
 257  */
 258 template<typename ValueType>
 259 static inline ValueType norm2(BasicVector<ValueType> v)
 260 {
 261     return v.norm2();
 262 }
 263
 264 /*! \brief
 265  * cross product for gmx::BasicVector
 266  */
 267 template<typename VectorType>
 268 static inline VectorType cross(const VectorType& a, const VectorType& b)
 269 {
 270     return a.cross(b);
 271 }
 272
 273 /*! \brief
 274  * dot product for gmx::BasicVector
 275  */
 276 template<typename ValueType>
 277 static inline ValueType dot(BasicVector<ValueType> a, BasicVector<ValueType> b)
 278 {
 279     return a.dot(b);
 280 }
 281
 282 /*! \brief
 283  * Multiply two vectors element by element and return the result.
 284  */
 285 template<typename VectorType>
 286 static inline VectorType scaleByVector(const VectorType& a, const VectorType& b)
 287 {
 288     return { a[0] * b[0], a[1] * b[1], a[2] * b[2] };
 289 }
 290
 291 /*! \brief
 292  * Return the element-wise minimum of two vectors.
 293  */
 294 template<typename VectorType>
 295 static inline VectorType elementWiseMin(const VectorType& a, const VectorType& b)
 296 {
 297     return { std::min(a[0], b[0]), std::min(a[1], b[1]), std::min(a[2], b[2]) };
 298 }
 299
 300 /*! \brief
 301  * Return the element-wise maximum of two vectors.
 302  */
 303 template<typename VectorType>
 304 static inline VectorType elementWiseMax(const VectorType& a, const VectorType& b)
 305 {
 306     return { std::max(a[0], b[0]), std::max(a[1], b[1]), std::max(a[2], b[2]) };
 307 }
 308
 309 /*! \brief
 310  * Casts a gmx::BasicVector array into an equivalent raw C array.
 311  */
 312 template<typename ValueType>
 313 static inline typename BasicVector<ValueType>::RawArray* as_vec_array(BasicVector<ValueType>* x)
 314 {
 315     return reinterpret_cast<typename BasicVector<ValueType>::RawArray*>(x);
 316 }
 317
 318 /*! \brief
 319  * Casts a gmx::BasicVector array into an equivalent raw C array.
 320  */
 321 template<typename ValueType>
 322 static inline const typename BasicVector<ValueType>::RawArray* as_vec_array(const BasicVector<ValueType>* x)
 323 {
 324     return reinterpret_cast<const typename BasicVector<ValueType>::RawArray*>(x);
 325 }
 326
 327 //! Shorthand for C++ `rvec`-equivalent type.
 328 typedef BasicVector<real> RVec;
 329 //! Shorthand for C++ `dvec`-equivalent type.
 330 typedef BasicVector<double> DVec;
 331 //! Shorthand for C++ `ivec`-equivalent type.
 332 typedef BasicVector<int> IVec;
 333 //! Casts a gmx::RVec array into an `rvec` array.
 334 static inline rvec* as_rvec_array(RVec* x)
 335 {
 336     return as_vec_array(x);
 337 }
 338 //! Casts a gmx::RVec array into an `rvec` array.
 339 static inline const rvec* as_rvec_array(const RVec* x)
 340 {
 341     return as_vec_array(x);
 342 }
 343 //! Casts a gmx::DVec array into an `Dvec` array.
 344 static inline dvec* as_dvec_array(DVec* x)
 345 {
 346     return as_vec_array(x);
 347 }
 348 //! Casts a gmx::IVec array into an `ivec` array.
 349 static inline ivec* as_ivec_array(IVec* x)
 350 {
 351     return as_vec_array(x);
 352 }
 353
 354
 355 //! Casts a gmx::DVec array into an `dvec` array.
 356 static inline const dvec* as_dvec_array(const DVec* x)
 357 {
 358     return as_vec_array(x);
 359 }
 360 //! Casts a gmx::IVec array into an `ivec` array.
 361 static inline const ivec* as_ivec_array(const IVec* x)
 362 {
 363     return as_vec_array(x);
 364 }
 365
 366 //! Shorthand for C++ `ivec`-equivalent type.
 367 typedef BasicVector<int> IVec;
 368
 369 } // namespace gmx
 370
 371 #endif // include guard