Transpose.hpp
Go to the documentation of this file.
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.