ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
Vector.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#pragma once
21
22/**
23 * @file
24 *
25 * @brief Vector implementation and trait types
26 * for boost qvm interoperability.
27 */
28
29#include <boost/qvm/deduce_vec.hpp>
30#include <boost/qvm/vec_traits.hpp>
31
32#include "utils/Array.hpp"
33#include "utils/get.hpp"
34
35#include <algorithm>
36#include <cassert>
37#include <cmath>
38#include <cstddef>
39#include <functional>
40#include <initializer_list>
41#include <iterator>
42#include <numeric>
43#include <span>
44#include <type_traits>
45#include <vector>
46
47namespace Utils {
48
49template <typename T, std::size_t N> class Vector : public Array<T, N> {
50 using Base = Array<T, N>;
51
52public:
53 using Base::Base;
54 using Array<T, N>::at;
55 using Array<T, N>::operator[];
56 using Array<T, N>::front;
57 using Array<T, N>::back;
58 using Array<T, N>::data;
59 using Array<T, N>::begin;
60 using Array<T, N>::cbegin;
61 using Array<T, N>::end;
62 using Array<T, N>::cend;
63 using Array<T, N>::empty;
64 using Array<T, N>::size;
65 using Array<T, N>::max_size;
66 using Array<T, N>::fill;
67 using Array<T, N>::broadcast;
68 Vector() noexcept = default;
69 Vector(Vector const &) = default;
70 Vector &operator=(Vector const &) = default;
71
72 void swap(Vector &rhs) { std::swap_ranges(begin(), end(), rhs.begin()); }
73
74private:
75 constexpr void copy_init(T const *first, T const *last) noexcept {
76 auto it = begin();
77 while (first != last) {
78 *it++ = *first++;
79 }
80 }
81
82public:
83 template <class Range>
84 explicit constexpr Vector(Range const &rng)
85 : Vector(std::begin(rng), std::end(rng)) {}
86 explicit constexpr Vector(T const (&v)[N]) noexcept : Base() {
87 copy_init(std::begin(v), std::end(v));
88 }
89
90 constexpr Vector(std::initializer_list<T> v) : Base() {
91 if (N != v.size()) {
92 throw std::length_error(
93 "Construction of Vector from Container of wrong length.");
94 }
95
96 copy_init(v.begin(), v.end());
97 }
98
99 template <typename InputIterator>
100 Vector(InputIterator first, InputIterator last) : Base() {
101 if (std::distance(first, last) == N) {
102 std::copy_n(first, N, begin());
103 } else {
104 throw std::length_error(
105 "Construction of Vector from Container of wrong length.");
106 }
107 }
108
109 /** @brief Create a vector that has all entries set to the same value. */
110 DEVICE_QUALIFIER static constexpr Vector<T, N>
111 broadcast(typename Base::value_type const &value) noexcept {
112 Vector<T, N> ret{};
113 for (std::size_t i = 0u; i != N; ++i) {
114 ret[i] = value;
115 }
116 return ret;
117 }
118
119 std::vector<T> as_vector() const { return std::vector<T>(begin(), end()); }
120
121 operator std::vector<T>() const { return as_vector(); }
122
123 constexpr std::span<T, N> as_span() const {
124 return std::span<T, N>(const_cast<T *>(begin()), size());
125 }
126
127 constexpr operator std::span<T, N>() const { return as_span(); }
128
129 template <class U> explicit operator Vector<U, N>() const {
130 Vector<U, N> ret;
131
132 std::transform(begin(), end(), ret.begin(),
133 [](auto const &e) { return static_cast<U>(e); });
134
135 return ret;
136 }
137
138 T norm2() const { return (*this) * (*this); }
139 T norm() const { return std::sqrt(norm2()); }
140
141 /*
142 * @brief Normalize the vector.
143 *
144 * Normalize the vector by its length,
145 * if not zero, otherwise the vector is unchanged.
146 */
147
149 auto const l = norm();
150 if (l > T(0)) {
151 for (std::size_t i = 0u; i < N; ++i)
152 this->operator[](i) /= l;
153 }
154
155 return *this;
156 }
157
158 Vector normalized() const { return (*this) / (*this).norm(); }
159};
160
161template <class T> using Vector3 = Vector<T, 3>;
162
163template <std::size_t N> using VectorXd = Vector<double, N>;
169
170template <std::size_t N> using VectorXf = Vector<float, N>;
172
173template <std::size_t N> using VectorXi = Vector<int, N>;
175
176namespace detail {
177template <std::size_t N, typename T, typename U, typename Op>
178auto binary_op(Vector<T, N> const &a, Vector<U, N> const &b, Op op) {
179 using std::declval;
180
181 using R = decltype(op(declval<T>(), declval<U>()));
182 Vector<R, N> ret;
183
184 std::transform(std::begin(a), std::end(a), std::begin(b), std::begin(ret),
185 op);
186
187 return ret;
188}
189
190template <std::size_t N, typename T, typename Op>
191Vector<T, N> &binary_op_assign(Vector<T, N> &a, Vector<T, N> const &b, Op op) {
192 std::transform(std::begin(a), std::end(a), std::begin(b), std::begin(a), op);
193 return a;
194}
195
196template <std::size_t N, typename T, typename Op>
197constexpr bool all_of(Vector<T, N> const &a, Vector<T, N> const &b, Op op) {
198 for (unsigned int i = 0; i < N; i++) {
199 /* Short circuit */
200 if (!static_cast<bool>(op(a[i], b[i]))) {
201 return false;
202 }
203 }
204
205 return true;
206}
207} // namespace detail
208
209template <std::size_t N, typename T>
210constexpr bool operator<(Vector<T, N> const &a, Vector<T, N> const &b) {
211 return detail::all_of(a, b, std::less<T>());
212}
213
214template <std::size_t N, typename T>
215constexpr bool operator>(Vector<T, N> const &a, Vector<T, N> const &b) {
216 return detail::all_of(a, b, std::greater<T>());
217}
218
219template <std::size_t N, typename T>
220constexpr bool operator<=(Vector<T, N> const &a, Vector<T, N> const &b) {
221 return detail::all_of(a, b, std::less_equal<T>());
222}
223
224template <std::size_t N, typename T>
225constexpr bool operator>=(Vector<T, N> const &a, Vector<T, N> const &b) {
226 return detail::all_of(a, b, std::greater_equal<T>());
227}
228
229template <std::size_t N, typename T>
230constexpr bool operator==(Vector<T, N> const &a, Vector<T, N> const &b) {
231 return detail::all_of(a, b, std::equal_to<T>());
232}
233
234template <std::size_t N, typename T>
235constexpr bool operator!=(Vector<T, N> const &a, Vector<T, N> const &b) {
236 return not(a == b);
237}
238
239template <std::size_t N, typename T, typename U>
240auto operator+(Vector<T, N> const &a, Vector<U, N> const &b) {
241 return detail::binary_op(a, b, std::plus<>());
242}
243
244template <std::size_t N, typename T>
246 return detail::binary_op_assign(a, b, std::plus<T>());
247}
248
249template <std::size_t N, typename T, typename U>
250auto operator-(Vector<T, N> const &a, Vector<U, N> const &b) {
251 return detail::binary_op(a, b, std::minus<>());
252}
253
254template <std::size_t N, typename T>
256 Vector<T, N> ret;
257
258 std::transform(std::begin(a), std::end(a), std::begin(ret), std::negate<T>());
259
260 return ret;
261}
262
263template <std::size_t N, typename T>
265 return detail::binary_op_assign(a, b, std::minus<T>());
266}
267
268/* Scalar multiplication */
269template <std::size_t N, typename T, class U,
270 std::enable_if_t<std::is_arithmetic_v<U>, bool> = true>
271auto operator*(U const &a, Vector<T, N> const &b) {
272 using R = decltype(a * std::declval<T>());
273 Vector<R, N> ret;
274
275 std::transform(std::begin(b), std::end(b), std::begin(ret),
276 [a](T const &val) { return a * val; });
277
278 return ret;
279}
280
281template <std::size_t N, typename T, class U,
282 std::enable_if_t<std::is_arithmetic_v<U>, bool> = true>
283auto operator*(Vector<T, N> const &b, U const &a) {
284 using R = decltype(std::declval<T>() * a);
285 Vector<R, N> ret;
286
287 std::transform(std::begin(b), std::end(b), std::begin(ret),
288 [a](T const &val) { return a * val; });
289
290 return ret;
291}
292
293template <std::size_t N, typename T>
295 std::transform(std::begin(b), std::end(b), std::begin(b),
296 [a](T const &val) { return a * val; });
297 return b;
298}
299
300/* Scalar division */
301template <std::size_t N, typename T>
302Vector<T, N> operator/(Vector<T, N> const &a, T const &b) {
303 Vector<T, N> ret;
304
305 std::transform(std::begin(a), std::end(a), ret.begin(),
306 [b](T const &val) { return val / b; });
307 return ret;
308}
309
310template <std::size_t N, typename T>
311Vector<T, N> operator/(T const &a, Vector<T, N> const &b) {
312 Vector<T, N> ret;
313
314 std::transform(std::begin(b), std::end(b), ret.begin(),
315 [a](T const &val) { return a / val; });
316 return ret;
317}
318
319template <std::size_t N, typename T>
321 std::transform(std::begin(a), std::end(a), std::begin(a),
322 [b](T const &val) { return val / b; });
323 return a;
324}
325
326namespace detail {
327template <class T> struct is_vector : std::false_type {};
328template <class T, std::size_t N>
329struct is_vector<Vector<T, N>> : std::true_type {};
330} // namespace detail
331
332/* Scalar product */
333template <std::size_t N, typename T, class U,
334 class = std::enable_if_t<not(detail::is_vector<T>::value or
335 detail::is_vector<U>::value)>>
336auto operator*(Vector<T, N> const &a, Vector<U, N> const &b) {
337 using std::declval;
338 using R = decltype(declval<T>() * declval<U>());
339
340 return std::inner_product(std::begin(a), std::end(a), std::begin(b), R{});
341}
342
343template <
344 std::size_t N, typename T, class U,
345 class = std::enable_if_t<std::is_integral_v<T> and std::is_integral_v<U>>>
346auto operator%(Vector<T, N> const &a, Vector<U, N> const &b) {
347 using std::declval;
348 using R = decltype(declval<T>() % declval<U>());
349 Vector<R, N> ret;
350
351 std::transform(std::begin(a), std::end(a), std::begin(b), std::begin(ret),
352 [](T const &ai, U const &bi) { return ai % bi; });
353
354 return ret;
355}
356
357/* Componentwise square root */
358template <std::size_t N, typename T> Vector<T, N> sqrt(Vector<T, N> const &a) {
359 using std::sqrt;
360 Vector<T, N> ret;
361
362 std::transform(std::begin(a), std::end(a), ret.begin(),
363 [](T const &v) { return sqrt(v); });
364
365 return ret;
366}
367
368template <class T>
370 return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
371 a[0] * b[1] - a[1] * b[0]};
372}
373
374// Product of array elements.
375template <class T, std::size_t N> T product(Vector<T, N> const &v) {
376 return std::accumulate(v.cbegin(), v.cend(), T{1}, std::multiplies<T>());
377}
378
379template <class T, class U, std::size_t N>
381 using std::declval;
382 using R = decltype(declval<T>() * declval<U>());
383
384 Vector<R, N> ret;
385 std::transform(a.cbegin(), a.cend(), b.cbegin(), ret.begin(),
386 [](auto const &ai, auto const &bi) { return ai * bi; });
387
388 return ret;
389}
390
391// specializations for when one or both operands is a scalar depending on
392// compile time features (e.g. when PARTICLE_ANISOTROPY is not enabled)
393template <class T, class U, std::size_t N,
394 class = std::enable_if_t<not(detail::is_vector<T>::value)>>
395auto hadamard_product(T const &a, Vector<U, N> const &b) {
396 using std::declval;
397 using R = decltype(declval<T>() * declval<U>());
398
399 Vector<R, N> ret = a * b;
400
401 return ret;
402}
403
404template <class T, class U, std::size_t N,
405 class = std::enable_if_t<not(detail::is_vector<U>::value)>>
406auto hadamard_product(Vector<T, N> const &a, U const &b) {
407 using std::declval;
408 using R = decltype(declval<T>() * declval<U>());
409
410 Vector<R, N> ret = a * b;
411
412 return ret;
413}
414
415template <typename T, typename U,
416 class = std::enable_if_t<not(detail::is_vector<T>::value or
417 detail::is_vector<U>::value)>>
418auto hadamard_product(T const &a, U const &b) {
419 return a * b;
420}
421
422template <class T, class U, std::size_t N>
424 using std::declval;
425 using R = decltype(declval<T>() * declval<U>());
426
427 Vector<R, N> ret;
428 std::transform(a.cbegin(), a.cend(), b.cbegin(), ret.begin(),
429 [](auto const &ai, auto const &bi) { return ai / bi; });
430
431 return ret;
432}
433
434// specializations for when one or both operands is a scalar depending on
435// compile time features (e.g. when PARTICLE_ANISOTROPY is not enabled)
436template <class T, class U, std::size_t N,
437 class = std::enable_if_t<not(detail::is_vector<U>::value)>>
438auto hadamard_division(Vector<T, N> const &a, U const &b) {
439 using std::declval;
440 using R = decltype(declval<T>() * declval<U>());
441
442 Vector<R, N> ret = a / b;
443
444 return ret;
445}
446
447template <class T, class U, std::size_t N,
448 class = std::enable_if_t<not(detail::is_vector<T>::value)>>
449auto hadamard_division(T const &a, Vector<U, N> const &b) {
450 using std::declval;
451 using R = decltype(declval<T>() * declval<U>());
452
453 Vector<R, N> ret;
454 std::transform(std::begin(b), std::end(b), ret.begin(),
455 [a](T const &bi) { return a / bi; });
456 return ret;
457}
458
459template <typename T, typename U,
460 class = std::enable_if_t<not(detail::is_vector<T>::value or
461 detail::is_vector<U>::value)>>
462auto hadamard_division(T const &a, U const &b) {
463 return a / b;
464}
465
466template <typename T> Vector<T, 3> unit_vector(unsigned int i) {
467 if (i == 0u)
468 return {T{1}, T{0}, T{0}};
469 if (i == 1u)
470 return {T{0}, T{1}, T{0}};
471 if (i == 2u)
472 return {T{0}, T{0}, T{1}};
473 throw std::domain_error("coordinate out of range");
474}
475
476/**
477 * @brief Meta function to turns a Vector<1, T> into T.
478 */
479template <typename T> struct decay_to_scalar {};
480template <typename T, std::size_t N> struct decay_to_scalar<Vector<T, N>> {
482};
483
484template <typename T> struct decay_to_scalar<Vector<T, 1>> {
485 using type = T;
486};
487
488template <std::size_t I, class T, std::size_t N>
489typename std::tuple_element<I, Vector<T, N>>::type &
490get(Vector<T, N> &a) noexcept {
491 return a[I];
492}
493
494template <std::size_t I, class T, std::size_t N>
495const typename std::tuple_element<I, Vector<T, N>>::type &
496get(Vector<T, N> const &a) noexcept {
497 return a[I];
498}
499
500} // namespace Utils
501
502template <std::size_t I, class T, std::size_t N>
503struct std::tuple_element<I, Utils::Vector<T, N>> {
504 static_assert(I < N, "Utils::Vector index must be in range");
505 using type = typename std::enable_if_t<(I < N), T>;
506};
507
508template <class T, std::size_t N>
509struct std::tuple_size<Utils::Vector<T, N>>
510 : std::integral_constant<std::size_t, N> {};
511
512namespace boost {
513namespace qvm {
514
515template <class T, std::size_t N> struct vec_traits<::Utils::Vector<T, N>> {
516
517 static constexpr std::size_t dim = N;
518 using scalar_type = T;
519
520 template <std::size_t I>
521 static constexpr inline scalar_type &write_element(::Utils::Vector<T, N> &v) {
522 return v[I];
523 }
524
525 template <std::size_t I>
526 static constexpr inline scalar_type
528 return v[I];
529 }
530
531 static inline scalar_type read_element_idx(std::size_t i,
532 ::Utils::Vector<T, N> const &v) {
533 return v[i];
534 }
535 static inline scalar_type &write_element_idx(std::size_t i,
537 return v[i];
538 }
539};
540
541template <typename T> struct deduce_vec<Utils::Vector<T, 3>, 3> {
542 using type = typename Utils::Vector<T, 3>;
543};
544
545} // namespace qvm
546} // namespace boost
547
550UTILS_ARRAY_BOOST_CLASS(Utils::Vector, N, object_serializable)
Array implementation with CUDA support.
#define UTILS_ARRAY_BOOST_MPI_T(Container, N)
Mark array types as MPI data types.
Definition array.hpp:50
#define UTILS_ARRAY_BOOST_BIT_S(Container, N)
Mark array types as MPI bitwise serializable.
Definition array.hpp:63
#define UTILS_ARRAY_BOOST_CLASS(Container, N, ImplementationLevel)
Redefinition of BOOST_CLASS_IMPLEMENTATION for array types.
Definition array.hpp:78
#define UTILS_ARRAY_BOOST_TRACK(Container, N, TrackingLevel)
Redefinition of BOOST_CLASS_TRACKING for array types.
Definition array.hpp:97
T norm2() const
Definition Vector.hpp:138
Vector() noexcept=default
Vector & normalize()
Definition Vector.hpp:148
T norm() const
Definition Vector.hpp:139
constexpr Vector(T const (&v)[N]) noexcept
Definition Vector.hpp:86
constexpr Vector(Range const &rng)
Definition Vector.hpp:84
void swap(Vector &rhs)
Definition Vector.hpp:72
Vector normalized() const
Definition Vector.hpp:158
std::vector< T > as_vector() const
Definition Vector.hpp:119
Vector(InputIterator first, InputIterator last)
Definition Vector.hpp:100
constexpr Vector(std::initializer_list< T > v)
Definition Vector.hpp:90
constexpr std::span< T, N > as_span() const
Definition Vector.hpp:123
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value) noexcept
Create a vector that has all entries set to the same value.
Definition Vector.hpp:111
#define DEVICE_QUALIFIER
auto operator+(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:240
Vector< T, 3 > unit_vector(unsigned int i)
Definition Vector.hpp:466
T product(Vector< T, N > const &v)
Definition Vector.hpp:375
Vector< T, N > & operator-=(Vector< T, N > &a, Vector< T, N > const &b)
Definition Vector.hpp:264
Vector< T, N > & operator/=(Vector< T, N > &a, T const &b)
Definition Vector.hpp:320
Vector< T, N > & operator+=(Vector< T, N > &a, Vector< T, N > const &b)
Definition Vector.hpp:245
constexpr bool operator<(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:210
constexpr bool operator>=(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:225
constexpr bool operator>(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:215
Vector< T, 3 > vector_product(Vector< T, 3 > const &a, Vector< T, 3 > const &b)
Definition Vector.hpp:369
Vector< T, N > operator/(Vector< T, N > const &a, T const &b)
Definition Vector.hpp:302
Quaternion< T > operator*(const U &b, const Quaternion< T > &a)
Product quaternion and arithmetic type.
constexpr bool operator==(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:230
constexpr bool operator!=(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:235
Vector< T, N > sqrt(Vector< T, N > const &a)
Definition Vector.hpp:358
std::tuple_element< I, Array< T, N > >::type & get(Array< T, N > &a) noexcept
Definition Array.hpp:219
constexpr bool operator<=(Vector< T, N > const &a, Vector< T, N > const &b)
Definition Vector.hpp:220
auto operator%(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:346
auto hadamard_division(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:423
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:380
Vector< T, N > & operator*=(Vector< T, N > &b, T const &a)
Definition Vector.hpp:294
auto operator-(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:250
DEVICE_QUALIFIER constexpr reference at(size_type i)
Definition Array.hpp:97
DEVICE_QUALIFIER constexpr bool empty() const noexcept
Definition Array.hpp:164
DEVICE_QUALIFIER constexpr reference back()
Definition Array.hpp:126
DEVICE_QUALIFIER constexpr pointer data() noexcept
Definition Array.hpp:132
DEVICE_QUALIFIER constexpr size_type max_size() const noexcept
Definition Array.hpp:168
DEVICE_QUALIFIER constexpr iterator begin() noexcept
Definition Array.hpp:140
DEVICE_QUALIFIER constexpr const_iterator cbegin() const noexcept
Definition Array.hpp:148
DEVICE_QUALIFIER constexpr const_iterator cend() const noexcept
Definition Array.hpp:160
DEVICE_QUALIFIER constexpr size_type size() const noexcept
Definition Array.hpp:166
DEVICE_QUALIFIER constexpr reference front()
Definition Array.hpp:122
DEVICE_QUALIFIER constexpr iterator end() noexcept
Definition Array.hpp:152
DEVICE_QUALIFIER void fill(const value_type &value)
Definition Array.hpp:170
Meta function to turns a Vector<1, T> into T.
Definition Vector.hpp:479
static constexpr scalar_type read_element(::Utils::Vector< T, N > const &v)
Definition Vector.hpp:527
static scalar_type & write_element_idx(std::size_t i, ::Utils::Vector< T, N > &v)
Definition Vector.hpp:535
static constexpr scalar_type & write_element(::Utils::Vector< T, N > &v)
Definition Vector.hpp:521
static scalar_type read_element_idx(std::size_t i, ::Utils::Vector< T, N > const &v)
Definition Vector.hpp:531
typename std::enable_if_t<(I< N), T > type
Definition Vector.hpp:505