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             : namespace detail {
      16             : void transpose_impl(double* matrix_transpose, const double* matrix,
      17             :                     int32_t number_of_rows, int32_t number_of_columns);
      18             : }  // namespace detail
      19             : 
      20             : 
      21             : /// \ingroup NumericalAlgorithmsGroup
      22             : /// \brief Function to compute transposed data.
      23             : ///
      24             : /// Transpose the data pointed to by `data`, writing the result to the
      25             : /// location pointed to by `result`.  See the ::transpose function for
      26             : /// a safer interface and for the meaning of the other arguments.
      27             : // This is not an overload of transpose because overload resolution
      28             : // would make non-intuitive choices of the return-by-pointer version
      29             : // below over this function.
      30             : template <typename T>
      31           1 : void raw_transpose(const gsl::not_null<T*> result, const T* const data,
      32             :                    const size_t chunk_size, const size_t number_of_chunks) {
      33             :   if constexpr (std::is_same_v<double, T>) {
      34             :     detail::transpose_impl(result, data, static_cast<int32_t>(number_of_chunks),
      35             :                            static_cast<int32_t>(chunk_size));
      36             :   } else {
      37             :     // The i outside loop order is faster, but that could be architecture
      38             :     // dependent and so may need updating in the future. Changing this made the
      39             :     // logical derivatives in 3D with 50 variables 20% faster.
      40             :     for (size_t i = 0; i < chunk_size; ++i) {
      41             :       for (size_t j = 0; j < number_of_chunks; ++j) {
      42             :         // clang-tidy: pointer arithmetic
      43             :         result.get()[j + number_of_chunks * i] =  // NOLINT
      44             :             data[i + chunk_size * j];             // NOLINT
      45             :       }
      46             :     }
      47             :   }
      48             : }
      49             : 
      50             : /// @{
      51             : /// \ingroup NumericalAlgorithmsGroup
      52             : /// \brief Function to compute transposed data.
      53             : ///
      54             : /// The primary use of this function is to rearrange the memory layout so that
      55             : /// another function can operate on contiguous chunks of memory.
      56             : ///
      57             : /// \requires `result.size()` to be the product of `number_of_chunks` and
      58             : /// `chunk_size`, `u.size()` to be equal or greater than `result.size()`,
      59             : /// and that both `result` and `u` have a `data()` member function.
      60             : ///
      61             : /// \details The container `u` holds a contiguous array of data,
      62             : /// treated as a sequence of `number_of_chunks` contiguous sets of
      63             : /// entries of size `chunk_size`.  The output `result` has its data
      64             : /// arranged such that the first `number_of_chunks` elements in
      65             : /// `result` will be the first element of each chunk of `u`.  The
      66             : /// last `number_of_chunks` elements in `result` will be the last
      67             : /// (i.e.  `chunk_size`-th) element of each chunk of `u`.  If
      68             : /// `u.size()` is greater than `result.size()` the extra elements
      69             : /// of `u` are ignored.
      70             : ///
      71             : /// \note This is equivalent to treating the first part of `u` as a matrix and
      72             : /// filling `result` (or the returned object) with the transpose of that
      73             : /// matrix.
      74             : ///
      75             : /// If `u` represents a block of data indexed by
      76             : /// \f$(x, y, z, \ldots)\f$ with the first index varying fastest,
      77             : /// transpose serves to rotate the indices.  If the extents are
      78             : /// \f$(X, Y, Z, \ldots)\f$, with product \f$N\f$, `transpose(u, X,
      79             : /// N/X)` reorders the data to be indexed \f$(y, z, \ldots, x)\f$,
      80             : /// `transpose(u, X*Y, N/X/Y)` reorders the data to be indexed
      81             : /// \f$(z, \ldots, x, y)\f$, etc.
      82             : ///
      83             : /// \example
      84             : /// \snippet Test_Transpose.cpp transpose_matrix
      85             : /// \snippet Test_Transpose.cpp transpose_by_not_null_example
      86             : /// \snippet Test_Transpose.cpp return_transpose_example
      87             : /// \snippet Test_Transpose.cpp partial_transpose_example
      88             : ///
      89             : /// \tparam U the type of data to be transposed
      90             : /// \tparam T the type of the transposed data
      91             : template <typename U, typename T>
      92           1 : void transpose(const gsl::not_null<T*> result, const U& u,
      93             :                const size_t chunk_size, const size_t number_of_chunks) {
      94             :   ASSERT(chunk_size * number_of_chunks == result->size(),
      95             :          "chunk_size = " << chunk_size << ", number_of_chunks = "
      96             :                          << number_of_chunks << ", size = " << result->size());
      97             :   ASSERT(result->size() <= u.size(),
      98             :          "result.size = " << result->size() << ", u.size = " << u.size());
      99             :   raw_transpose(make_not_null(result->data()), u.data(), chunk_size,
     100             :                 number_of_chunks);
     101             : }
     102             : 
     103             : template <typename U, typename T = U>
     104           1 : T transpose(const U& u, const size_t chunk_size,
     105             :             const size_t number_of_chunks) {
     106             :   T t = make_with_value<T>(u, 0.0);
     107             :   transpose(make_not_null(&t), u, chunk_size, number_of_chunks);
     108             :   return t;
     109             : }
     110             : /// @}