Line data    Source code

       1           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 "Utilities/ErrorHandling/Assert.hpp"
      12             : #include "Utilities/Gsl.hpp"
      13             : #include "Utilities/MakeWithValue.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           1 : void raw_transpose(const gsl::not_null<T*> result, const T* const data,
      26             :                    const size_t chunk_size, const size_t number_of_chunks) {
      27             :   // The i outside loop order is faster, but that could be architecture
      28             :   // dependent and so may need updating in the future. Changing this made the
      29             :   // logical derivatives in 3D with 50 variables 20% faster.
      30             :   for (size_t i = 0; i < chunk_size; ++i) {
      31             :     for (size_t j = 0; j < number_of_chunks; ++j) {
      32             :       // clang-tidy: pointer arithmetic
      33             :       result.get()[j + number_of_chunks * i] =  // NOLINT
      34             :           data[i + chunk_size * j];             // NOLINT
      35             :     }
      36             :   }
      37             : }
      38             : 
      39             : /// @{
      40             : /// \ingroup NumericalAlgorithmsGroup
      41             : /// \brief Function to compute transposed data.
      42             : ///
      43             : /// The primary use of this function is to rearrange the memory layout so that
      44             : /// another function can operate on contiguous chunks of memory.
      45             : ///
      46             : /// \requires `result.size()` to be the product of `number_of_chunks` and
      47             : /// `chunk_size`, `u.size()` to be equal or greater than `result.size()`,
      48             : /// and that both `result` and `u` have a `data()` member function.
      49             : ///
      50             : /// \details The container `u` holds a contiguous array of data,
      51             : /// treated as a sequence of `number_of_chunks` contiguous sets of
      52             : /// entries of size `chunk_size`.  The output `result` has its data
      53             : /// arranged such that the first `number_of_chunks` elements in
      54             : /// `result` will be the first element of each chunk of `u`.  The
      55             : /// last `number_of_chunks` elements in `result` will be the last
      56             : /// (i.e.  `chunk_size`-th) element of each chunk of `u`.  If
      57             : /// `u.size()` is greater than `result.size()` the extra elements
      58             : /// of `u` are ignored.
      59             : ///
      60             : /// \note This is equivalent to treating the first part of `u` as a matrix and
      61             : /// filling `result` (or the returned object) with the transpose of that
      62             : /// matrix.
      63             : ///
      64             : /// If `u` represents a block of data indexed by
      65             : /// \f$(x, y, z, \ldots)\f$ with the first index varying fastest,
      66             : /// transpose serves to rotate the indices.  If the extents are
      67             : /// \f$(X, Y, Z, \ldots)\f$, with product \f$N\f$, `transpose(u, X,
      68             : /// N/X)` reorders the data to be indexed \f$(y, z, \ldots, x)\f$,
      69             : /// `transpose(u, X*Y, N/X/Y)` reorders the data to be indexed
      70             : /// \f$(z, \ldots, x, y)\f$, etc.
      71             : ///
      72             : /// \example
      73             : /// \snippet Test_Transpose.cpp transpose_matrix
      74             : /// \snippet Test_Transpose.cpp transpose_by_not_null_example
      75             : /// \snippet Test_Transpose.cpp return_transpose_example
      76             : /// \snippet Test_Transpose.cpp partial_transpose_example
      77             : ///
      78             : /// \tparam U the type of data to be transposed
      79             : /// \tparam T the type of the transposed data
      80             : template <typename U, typename T>
      81           1 : void transpose(const gsl::not_null<T*> result, const U& u,
      82             :                const size_t chunk_size, const size_t number_of_chunks) {
      83             :   ASSERT(chunk_size * number_of_chunks == result->size(),
      84             :          "chunk_size = " << chunk_size << ", number_of_chunks = "
      85             :                          << number_of_chunks << ", size = " << result->size());
      86             :   ASSERT(result->size() <= u.size(),
      87             :          "result.size = " << result->size() << ", u.size = " << u.size());
      88             :   raw_transpose(make_not_null(result->data()), u.data(), chunk_size,
      89             :                 number_of_chunks);
      90             : }
      91             : 
      92             : template <typename U, typename T = U>
      93           1 : T transpose(const U& u, const size_t chunk_size,
      94             :             const size_t number_of_chunks) {
      95             :   T t = make_with_value<T>(u, 0.0);
      96             :   transpose(make_not_null(&t), u, chunk_size, number_of_chunks);
      97             :   return t;
      98             : }
      99             : /// @}