Line data    Source code

       1             : // SPDX-License-Identifier: BSD-3-Clause
       2             : // Copyright (c) 2023 Scipp contributors (https://github.com/scipp)
       3             : /// @file
       4             : /// @author Jan-Lukas Wynen
       5             : #pragma once
       6             : 
       7             : #include "scipp-core_export.h"
       8             : #include "scipp/common/index_composition.h"
       9             : #include "scipp/core/dimensions.h"
      10             : #include "scipp/core/strides.h"
      11             : 
      12             : /*
      13             :  * See page on Multi-Dimensional Indexing in developer documentation
      14             :  * for an explanation of the implementation.
      15             :  */
      16             : 
      17             : namespace scipp::core {
      18             : 
      19             : /// A flat index into a multi-dimensional view.
      20             : class SCIPP_CORE_EXPORT ViewIndex {
      21             : public:
      22             :   ViewIndex(const Dimensions &target_dimensions, const Strides &strides);
      23             : 
      24     2262396 :   constexpr void increment_outer() noexcept {
      25     2262396 :     for (scipp::index d = 0;
      26     4971445 :          (d < NDIM_OP_MAX - 1) && (m_coord[d] == m_shape[d]); ++d) {
      27     2709049 :       m_memory_index += m_delta[d + 1];
      28     2709049 :       ++m_coord[d + 1];
      29     2709049 :       m_coord[d] = 0;
      30             :     }
      31     2262396 :   }
      32    22686268 :   constexpr void increment() noexcept {
      33    22686268 :     m_memory_index += m_delta[0];
      34    22686268 :     ++m_coord[0];
      35    22686268 :     if (m_coord[0] == m_shape[0])
      36     2262396 :       increment_outer();
      37    22686268 :     ++m_view_index;
      38    22686268 :   }
      39             : 
      40     2416730 :   void set_index(const scipp::index index) noexcept {
      41     2416730 :     m_view_index = index;
      42     2416730 :     extract_indices(index, m_shape.begin(), m_shape.begin() + m_ndim,
      43             :                     m_coord.begin());
      44     4833460 :     m_memory_index = flat_index_from_strides(
      45     2416730 :         m_strides.begin(), m_strides.begin() + m_ndim, m_coord.begin());
      46     2416730 :   }
      47             : 
      48    22783829 :   [[nodiscard]] constexpr scipp::index get() const noexcept {
      49    22783829 :     return m_memory_index;
      50             :   }
      51     2199829 :   [[nodiscard]] constexpr scipp::index index() const noexcept {
      52     2199829 :     return m_view_index;
      53             :   }
      54             : 
      55    11664651 :   constexpr bool operator==(const ViewIndex &other) const noexcept {
      56    11664651 :     return m_view_index == other.m_view_index;
      57             :   }
      58             :   constexpr bool operator!=(const ViewIndex &other) const noexcept {
      59             :     return m_view_index != other.m_view_index;
      60             :   }
      61             : 
      62             : private:
      63             :   /// Index into memory.
      64             :   scipp::index m_memory_index{0};
      65             :   /// Index in iteration dimensions.
      66             :   scipp::index m_view_index{0};
      67             :   /// Steps in memory to advance one element.
      68             :   std::array<scipp::index, NDIM_OP_MAX> m_delta = {};
      69             :   /// Multi-dimensional index in iteration dimensions.
      70             :   std::array<scipp::index, NDIM_OP_MAX> m_coord = {};
      71             :   /// Shape in iteration dimensions.
      72             :   std::array<scipp::index, NDIM_OP_MAX> m_shape = {};
      73             :   /// Strides in memory.
      74             :   std::array<scipp::index, NDIM_OP_MAX> m_strides = {};
      75             :   /// Number of dimensions.
      76             :   int32_t m_ndim;
      77             : };
      78             : // NOTE:
      79             : // We investigated different containers for the m_delta, m_coord & m_extent
      80             : // arrays, and their impact on performance when iterating over a variable
      81             : // view.
      82             : // Using std::array or C-style arrays give good performance (7.5 Gb/s) as long
      83             : // as a range based loop is used:
      84             : //
      85             : //   for ( const auto x : view ) {
      86             : //
      87             : // If a loop which explicitly accesses the begin() and end() of the container
      88             : // is used, e.g.
      89             : //
      90             : //   for ( auto it = view.begin(); it != view.end(); ++it ) {
      91             : //
      92             : // then the results differ widely.
      93             : // - using std::array is 80x slower than above, at ~90 Mb/s
      94             : // - using C-style arrays is 20x slower than above, at ~330 Mb/s
      95             : //
      96             : // We can recover the maximum performance by storing the view.end() in a
      97             : // variable, e.g.
      98             : //
      99             : //   auto iend = view.end();
     100             : //   for ( auto it = view.begin(); it != iend; ++it ) {
     101             : //
     102             : // for both std::array and C-style arrays.
     103             : //
     104             : // Finally, when using C-style arrays, we get a compilation warning from L37
     105             : //
     106             : //   m_delta[d] -= m_delta[d2];
     107             : //
     108             : // with the GCC compiler:
     109             : //
     110             : //  warning: array subscript is above array bounds [-Warray-bounds]
     111             : //
     112             : // which disappears when switching to std::array. This warning is not given
     113             : // by the CLANG compiler, and is not fully understood as d2 is always less
     114             : // than d and should never overflow the array bounds.
     115             : // We decided to go with std::array as our final choice to avoid the warning,
     116             : // as the performance is identical to C-style arrays, as long as range based
     117             : // loops are used.
     118             : 
     119             : } // namespace scipp::core