Transpose.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines function transpose
6 
7 #pragma once
8 
9 #include <cstddef>
10 
11 #include "ErrorHandling/Assert.hpp"
12 #include "Utilities/Gsl.hpp"
14 
15 /// \ingroup NumericalAlgorithmsGroup
16 /// \brief Function to compute transposed data.
17 ///
18 /// Transpose the data pointed to by `data`, writing the result to the
19 /// location pointed to by `result`. See the ::transpose function for
20 /// a safer interface and for the meaning of the other arguments.
21 // This is not an overload of transpose because overload resolution
22 // would make non-intuitive choices of the return-by-pointer version
23 // below over this function.
24 template <typename T>
25 void raw_transpose(const gsl::not_null<T*> result, const T* const data,
26  const size_t chunk_size,
27  const size_t number_of_chunks) noexcept {
28  // The i outside loop order is faster, but that could be architecture
29  // dependent and so may need updating in the future. Changing this made the
30  // logical derivatives in 3D with 50 variables 20% faster.
31  for (size_t i = 0; i < chunk_size; ++i) {
32  for (size_t j = 0; j < number_of_chunks; ++j) {
33  // clang-tidy: pointer arithmetic
34  result.get()[j + number_of_chunks * i] = // NOLINT
35  data[i + chunk_size * j]; // NOLINT
36  }
37  }
38 }
39 
40 // @{
41 /// \ingroup NumericalAlgorithmsGroup
42 /// \brief Function to compute transposed data.
43 ///
44 /// The primary use of this function is to rearrange the memory layout so that
45 /// another function can operate on contiguous chunks of memory.
46 ///
47 /// \requires `result.size()` to be the product of `number_of_chunks` and
48 /// `chunk_size`, `u.size()` to be equal or greater than `result.size()`,
49 /// and that both `result` and `u` have a `data()` member function.
50 ///
51 /// \details The container `u` holds a contiguous array of data,
52 /// treated as a sequence of `number_of_chunks` contiguous sets of
53 /// entries of size `chunk_size`. The output `result` has its data
54 /// arranged such that the first `number_of_chunks` elements in
55 /// `result` will be the first element of each chunk of `u`. The
56 /// last `number_of_chunks` elements in `result` will be the last
57 /// (i.e. `chunk_size`-th) element of each chunk of `u`. If
58 /// `u.size()` is greater than `result.size()` the extra elements
59 /// of `u` are ignored.
60 ///
61 /// \note This is equivalent to treating the first part of `u` as a matrix and
62 /// filling `result` (or the returned object) with the transpose of that
63 /// matrix.
64 ///
65 /// If `u` represents a block of data indexed by
66 /// \f$(x, y, z, \ldots)\f$ with the first index varying fastest,
67 /// transpose serves to rotate the indices. If the extents are
68 /// \f$(X, Y, Z, \ldots)\f$, with product \f$N\f$, `transpose(u, X,
69 /// N/X)` reorders the data to be indexed \f$(y, z, \ldots, x)\f$,
70 /// `transpose(u, X*Y, N/X/Y)` reorders the data to be indexed
71 /// \f$(z, \ldots, x, y)\f$, etc.
72 ///
73 /// \example
74 /// \snippet Test_Transpose.cpp transpose_matrix
75 /// \snippet Test_Transpose.cpp transpose_by_not_null_example
76 /// \snippet Test_Transpose.cpp return_transpose_example
77 /// \snippet Test_Transpose.cpp partial_transpose_example
78 ///
79 /// \tparam U the type of data to be transposed
80 /// \tparam T the type of the transposed data
81 template <typename U, typename T>
82 void transpose(const gsl::not_null<T*> result, const U& u,
83  const size_t chunk_size,
84  const size_t number_of_chunks) noexcept {
85  ASSERT(chunk_size * number_of_chunks == result->size(),
86  "chunk_size = " << chunk_size << ", number_of_chunks = "
87  << number_of_chunks << ", size = " << result->size());
88  ASSERT(result->size() <= u.size(),
89  "result.size = " << result->size() << ", u.size = " << u.size());
90  raw_transpose(make_not_null(result->data()), u.data(), chunk_size,
91  number_of_chunks);
92 }
93 
94 template <typename U, typename T = U>
95 T transpose(const U& u, const size_t chunk_size,
96  const size_t number_of_chunks) noexcept {
97  T t = make_with_value<T>(u, 0.0);
98  transpose(make_not_null(&t), u, chunk_size, number_of_chunks);
99  return t;
100 }
101 // @}
void transpose(const gsl::not_null< T *> result, const U &u, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data.
Definition: Transpose.hpp:82
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Defines macro ASSERT.
void raw_transpose(const gsl::not_null< T *> result, const T *const data, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data.
Definition: Transpose.hpp:25
Defines functions and classes from the GSL.
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
Defines make_with_value.