api/legacy/include/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,2021, 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     //! Return whether all elements compare equal
 127     bool operator==(const BasicVector<ValueType>& right)
 128     {
 129         return x_[0] == right[0] && x_[1] == right[1] && x_[2] == right[2];
 130     }
 131     //! Return whether any elements compare unequal
 132     bool operator!=(const BasicVector<ValueType>& right)
 133     {
 134         return x_[0] != right[0] || x_[1] != right[1] || x_[2] != right[2];
 135     }
 136     //! Allow inplace addition for BasicVector
 137     BasicVector<ValueType>& operator+=(const BasicVector<ValueType>& right)
 138     {
 139         return *this = *this + right;
 140     }
 141     //! Allow inplace subtraction for BasicVector
 142     BasicVector<ValueType>& operator-=(const BasicVector<ValueType>& right)
 143     {
 144         return *this = *this - right;
 145     }
 146     //! Allow vector addition
 147     BasicVector<ValueType> operator+(const BasicVector<ValueType>& right) const
 148     {
 149         return { x_[0] + right[0], x_[1] + right[1], x_[2] + right[2] };
 150     }
 151     //! Allow vector subtraction
 152     BasicVector<ValueType> operator-(const BasicVector<ValueType>& right) const
 153     {
 154         return { x_[0] - right[0], x_[1] - right[1], x_[2] - right[2] };
 155     }
 156     //! Allow vector scalar division
 157     BasicVector<ValueType> operator/(const ValueType& right) const
 158     {
 159         assert((right != 0 && "Cannot divide by zero"));
 160
 161         return *this * (1 / right);
 162     }
 163     //! Scale vector by a scalar
 164     BasicVector<ValueType>& operator*=(const ValueType& right)
 165     {
 166         x_[0] *= right;
 167         x_[1] *= right;
 168         x_[2] *= right;
 169
 170         return *this;
 171     }
 172     //! Divide vector by a scalar
 173     BasicVector<ValueType>& operator/=(const ValueType& right)
 174     {
 175         assert((right != 0 && "Cannot divide by zero"));
 176
 177         return *this *= 1 / right;
 178     }
 179     //! Return dot product
 180     ValueType dot(const BasicVector<ValueType>& right) const
 181     {
 182         return x_[0] * right[0] + x_[1] * right[1] + x_[2] * right[2];
 183     }
 184
 185     //! Allow vector vector multiplication (cross product)
 186     BasicVector<ValueType> cross(const BasicVector<ValueType>& right) const
 187     {
 188         return { x_[YY] * right.x_[ZZ] - x_[ZZ] * right.x_[YY],
 189                  x_[ZZ] * right.x_[XX] - x_[XX] * right.x_[ZZ],
 190                  x_[XX] * right.x_[YY] - x_[YY] * right.x_[XX] };
 191     }
 192
 193     //! Return normalized to unit vector
 194     BasicVector<ValueType> unitVector() const
 195     {
 196         const ValueType vectorNorm = norm();
 197         assert((vectorNorm != 0 && "unitVector() should not be called with a zero vector"));
 198
 199         return *this / vectorNorm;
 200     }
 201
 202     //! Length^2 of vector
 203     ValueType norm2() const { return dot(*this); }
 204
 205     //! Norm or length of vector
 206     ValueType norm() const { return std::sqrt(norm2()); }
 207
 208     //! cast to RVec
 209     BasicVector<real> toRVec() const { return { real(x_[0]), real(x_[1]), real(x_[2]) }; }
 210
 211     //! cast to IVec
 212     BasicVector<int> toIVec() const
 213     {
 214         return { static_cast<int>(x_[0]), static_cast<int>(x_[1]), static_cast<int>(x_[2]) };
 215     }
 216
 217     //! cast to DVec
 218     BasicVector<double> toDVec() const { return { double(x_[0]), double(x_[1]), double(x_[2]) }; }
 219
 220     //! Converts to a raw C array where implicit conversion does not work.
 221     RawArray& as_vec() { return x_; }
 222     //! Converts to a raw C array where implicit conversion does not work.
 223     const RawArray& as_vec() const { return x_; }
 224     //! Makes BasicVector usable in contexts where a raw C array is expected.
 225     operator RawArray&() { return x_; }
 226     //! Makes BasicVector usable in contexts where a raw C array is expected.
 227     operator const RawArray&() const { return x_; }
 228
 229 private:
 230     RawArray x_;
 231 };
 232
 233 //! Allow vector scalar multiplication
 234 template<typename ValueType>
 235 BasicVector<ValueType> operator*(const BasicVector<ValueType>& basicVector, const ValueType& scalar)
 236 {
 237     return { basicVector[0] * scalar, basicVector[1] * scalar, basicVector[2] * scalar };
 238 }
 239
 240 //! Allow scalar vector multiplication
 241 template<typename ValueType>
 242 BasicVector<ValueType> operator*(const ValueType& scalar, const BasicVector<ValueType>& basicVector)
 243 {
 244     return { scalar * basicVector[0], scalar * basicVector[1], scalar * basicVector[2] };
 245 }
 246
 247 /*! \brief
 248  * unitv for gmx::BasicVector
 249  */
 250 template<typename VectorType>
 251 static inline VectorType unitVector(const VectorType& v)
 252 {
 253     return v.unitVector();
 254 }
 255
 256 /*! \brief
 257  * norm for gmx::BasicVector
 258  */
 259 template<typename ValueType>
 260 static inline ValueType norm(BasicVector<ValueType> v)
 261 {
 262     return v.norm();
 263 }
 264
 265 /*! \brief
 266  * Square of the vector norm for gmx::BasicVector
 267  */
 268 template<typename ValueType>
 269 static inline ValueType norm2(BasicVector<ValueType> v)
 270 {
 271     return v.norm2();
 272 }
 273
 274 /*! \brief
 275  * cross product for gmx::BasicVector
 276  */
 277 template<typename VectorType>
 278 static inline VectorType cross(const VectorType& a, const VectorType& b)
 279 {
 280     return a.cross(b);
 281 }
 282
 283 /*! \brief
 284  * dot product for gmx::BasicVector
 285  */
 286 template<typename ValueType>
 287 static inline ValueType dot(BasicVector<ValueType> a, BasicVector<ValueType> b)
 288 {
 289     return a.dot(b);
 290 }
 291
 292 /*! \brief
 293  * Multiply two vectors element by element and return the result.
 294  */
 295 template<typename VectorType>
 296 static inline VectorType scaleByVector(const VectorType& a, const VectorType& b)
 297 {
 298     return { a[0] * b[0], a[1] * b[1], a[2] * b[2] };
 299 }
 300
 301 /*! \brief
 302  * Return the element-wise minimum of two vectors.
 303  */
 304 template<typename VectorType>
 305 static inline VectorType elementWiseMin(const VectorType& a, const VectorType& b)
 306 {
 307     return { std::min(a[0], b[0]), std::min(a[1], b[1]), std::min(a[2], b[2]) };
 308 }
 309
 310 /*! \brief
 311  * Return the element-wise maximum of two vectors.
 312  */
 313 template<typename VectorType>
 314 static inline VectorType elementWiseMax(const VectorType& a, const VectorType& b)
 315 {
 316     return { std::max(a[0], b[0]), std::max(a[1], b[1]), std::max(a[2], b[2]) };
 317 }
 318
 319 /*! \brief
 320  * Casts a gmx::BasicVector array into an equivalent raw C array.
 321  */
 322 template<typename ValueType>
 323 static inline typename BasicVector<ValueType>::RawArray* as_vec_array(BasicVector<ValueType>* x)
 324 {
 325     return reinterpret_cast<typename BasicVector<ValueType>::RawArray*>(x);
 326 }
 327
 328 /*! \brief
 329  * Casts a gmx::BasicVector array into an equivalent raw C array.
 330  */
 331 template<typename ValueType>
 332 static inline const typename BasicVector<ValueType>::RawArray* as_vec_array(const BasicVector<ValueType>* x)
 333 {
 334     return reinterpret_cast<const typename BasicVector<ValueType>::RawArray*>(x);
 335 }
 336
 337 //! Shorthand for C++ `rvec`-equivalent type.
 338 typedef BasicVector<real> RVec;
 339 //! Shorthand for C++ `dvec`-equivalent type.
 340 typedef BasicVector<double> DVec;
 341 //! Shorthand for C++ `ivec`-equivalent type.
 342 typedef BasicVector<int> IVec;
 343 //! Casts a gmx::RVec array into an `rvec` array.
 344 static inline rvec* as_rvec_array(RVec* x)
 345 {
 346     return as_vec_array(x);
 347 }
 348 //! Casts a gmx::RVec array into an `rvec` array.
 349 static inline const rvec* as_rvec_array(const RVec* x)
 350 {
 351     return as_vec_array(x);
 352 }
 353 //! Casts a gmx::DVec array into an `Dvec` array.
 354 static inline dvec* as_dvec_array(DVec* x)
 355 {
 356     return as_vec_array(x);
 357 }
 358 //! Casts a gmx::IVec array into an `ivec` array.
 359 static inline ivec* as_ivec_array(IVec* x)
 360 {
 361     return as_vec_array(x);
 362 }
 363
 364
 365 //! Casts a gmx::DVec array into an `dvec` array.
 366 static inline const dvec* as_dvec_array(const DVec* x)
 367 {
 368     return as_vec_array(x);
 369 }
 370 //! Casts a gmx::IVec array into an `ivec` array.
 371 static inline const ivec* as_ivec_array(const IVec* x)
 372 {
 373     return as_vec_array(x);
 374 }
 375
 376 //! Shorthand for C++ `ivec`-equivalent type.
 377 typedef BasicVector<int> IVec;
 378
 379 } // namespace gmx
 380
 381 #endif // include guard