Bugfix and optimization for 1D Array/ArrayView

RELEASE-NOTES.md

@@ -46,10 +46,13 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
 
 ### Changed
 - Updates CMake code check targets to only use checked in files (via `git ls-files`, when available)
+- Core: Optimization for axom::Array indirection -- since the stride is always 1, we can remove the runtime multiplication
 
 ### Fixed
 - Primal: Fixes signs of `compute_moments` to match orientation convention in `primal::evaluate_area_integral`
 - Quest: Improves error handling/reporting when loading an invalid c2c contour
+- Core: ArrayView assigments/copies now copy the stride
+- Core: Array construction from strided ArrayView now correctly copies the strided elements
 
 ## [Version 0.14.0] - Release date 2026-03-31
 

src/axom/bump/views/Shapes.hpp

@@ -87,7 +87,7 @@ struct PointTraits
 
   AXOM_HOST_DEVICE constexpr static axom::StackArray<IndexType, 1> getFace(IndexType faceIndex)
   {
-#if !defined(AXOM_DEBUG)
+#ifdef NDEBUG
     AXOM_UNUSED_VAR(faceIndex);
 #endif
     assert(faceIndex == 0);
@@ -136,7 +136,7 @@ struct LineTraits
 
   AXOM_HOST_DEVICE constexpr static axom::StackArray<IndexType, 2> getFace(IndexType faceIndex)
   {
-#if !defined(AXOM_DEBUG)
+#ifdef NDEBUG
     AXOM_UNUSED_VAR(faceIndex);
 #endif
     assert(faceIndex == 0);
@@ -190,7 +190,7 @@ struct TriTraits
 
   AXOM_HOST_DEVICE constexpr static axom::StackArray<IndexType, 3> getFace(IndexType faceIndex)
   {
-#if !defined(AXOM_DEBUG)
+#ifdef NDEBUG
     AXOM_UNUSED_VAR(faceIndex);
 #endif
     assert(faceIndex == 0);
@@ -245,7 +245,7 @@ struct QuadTraits
 
   AXOM_HOST_DEVICE constexpr static axom::StackArray<IndexType, 4> getFace(IndexType faceIndex)
   {
-#if !defined(AXOM_DEBUG)
+#ifdef NDEBUG
     AXOM_UNUSED_VAR(faceIndex);
 #endif
     assert(faceIndex == 0);

src/axom/core/Array.hpp

@@ -983,17 +983,17 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE, StoragePolicy>>, pro
   void initialize(IndexType num_elements, IndexType capacity, bool should_default_construct = true);
 
   /*!
-   * \brief Helper function for initializing an Array instance with an existing
-   *  range of elements.
+   * \brief Helper function for initializing an Array instance with an existing range of elements.
    *
    * \param [in] data pointer to the existing array of elements
    * \param [in] num_elements the number of elements in the existing array
+   * \param [in] src_stride the inter-element stride between elements of the existing array
    * \param [in] data_space the memory space in which data has been allocated
-   * \param [in] user_provided_allocator true if the Array's allocator ID was
-   *  provided by the user
+   * \param [in] user_provided_allocator true if the Array's allocator ID was provided by the user
    */
   void initialize_from_other(const T* data,
                              IndexType num_elements,
+                             IndexType src_stride,
                              MemorySpace data_space,
                              bool user_provided_allocator);
 
@@ -1202,7 +1202,7 @@ Array<T, DIM, SPACE, StoragePolicy>::Array(std::initializer_list<T> elems, int a
   : m_allocator_id(allocator_id)
   , m_arrayOps(m_allocator_id, m_executeOnGPU)
 {
-  initialize_from_other(elems.begin(), elems.size(), MemorySpace::Dynamic, true);
+  initialize_from_other(elems.begin(), elems.size(), 1 /* stride */, MemorySpace::Dynamic, true);
 }
 
 //------------------------------------------------------------------------------
@@ -1275,6 +1275,7 @@ Array<T, DIM, SPACE, StoragePolicy>::Array(const ArrayBase<T, DIM, OtherArrayTyp
 {
   initialize_from_other(static_cast<const OtherArrayType&>(other).data(),
                         static_cast<const OtherArrayType&>(other).size(),
+                        other.minStride(),
                         axom::detail::getAllocatorSpace(m_allocator_id),
                         false);
 }
@@ -1289,6 +1290,7 @@ Array<T, DIM, SPACE, StoragePolicy>::Array(const ArrayBase<const T, DIM, OtherAr
 {
   initialize_from_other(static_cast<const OtherArrayType&>(other).data(),
                         static_cast<const OtherArrayType&>(other).size(),
+                        other.minStride(),
                         axom::detail::getAllocatorSpace(m_allocator_id),
                         false);
 }
@@ -1306,6 +1308,7 @@ Array<T, DIM, SPACE, StoragePolicy>::Array(const ArrayBase<T, DIM, OtherArrayTyp
 
   initialize_from_other(static_cast<const OtherArrayType&>(other).data(),
                         static_cast<const OtherArrayType&>(other).size(),
+                        other.minStride(),
                         axom::detail::getAllocatorSpace(src_allocator),
                         true);
 }
@@ -1323,6 +1326,7 @@ Array<T, DIM, SPACE, StoragePolicy>::Array(const ArrayBase<const T, DIM, OtherAr
 
   initialize_from_other(static_cast<const OtherArrayType&>(other).data(),
                         static_cast<const OtherArrayType&>(other).size(),
+                        other.minStride(),
                         axom::detail::getAllocatorSpace(src_allocator),
                         true);
 }
@@ -1388,7 +1392,7 @@ inline void Array<T, DIM, SPACE, StoragePolicy>::assign(InputIt first, InputIt l
   {
     tmp.push_back(*it);
   }
-  initialize_from_other(tmp.data(), tmp.size(), MemorySpace::Dynamic, true);
+  initialize_from_other(tmp.data(), tmp.size(), 1 /* stride */, MemorySpace::Dynamic, true);
 }
 
 //------------------------------------------------------------------------------
@@ -1715,6 +1719,7 @@ template <typename T, int DIM, MemorySpace SPACE, typename StoragePolicy>
 inline void Array<T, DIM, SPACE, StoragePolicy>::initialize_from_other(
   const T* other_data,
   IndexType num_elements,
+  IndexType src_stride,
   MemorySpace other_data_space,
   bool AXOM_DEBUG_PARAM(user_provided_allocator))
 {
@@ -1734,9 +1739,15 @@ inline void Array<T, DIM, SPACE, StoragePolicy>::initialize_from_other(
   m_executeOnGPU = axom::isDeviceAllocator(m_allocator_id);
   m_arrayOps = OpHelper {m_allocator_id, m_executeOnGPU};
   this->setCapacity(num_elements);
-  // Use fill_range to ensure that copy constructors are invoked for each
-  // element.
-  m_arrayOps.fill_range(m_data, 0, num_elements, other_data, other_data_space);
+  // Use strided copy when necessary, otherwise use efficient contiguous copy
+  if(src_stride == 1)
+  {
+    m_arrayOps.fill_range(m_data, 0, num_elements, other_data, other_data_space);
+  }
+  else
+  {
+    m_arrayOps.fill_range_strided(m_data, 0, num_elements, other_data, src_stride, other_data_space);
+  }
   this->updateNumElements(num_elements);
 }
 

src/axom/core/ArrayBase.hpp

@@ -7,15 +7,15 @@
 #ifndef AXOM_ARRAYBASE_HPP_
 #define AXOM_ARRAYBASE_HPP_
 
-#include "axom/config.hpp"                    // for compile-time defines
-#include "axom/core/Macros.hpp"               // for axom macros
-#include "axom/core/MDMapping.hpp"            // for index conversion
-#include "axom/core/memory_management.hpp"    // for memory allocation functions
-#include "axom/core/utilities/Utilities.hpp"  // for processAbort()
-#include "axom/core/Types.hpp"                // for IndexType definition
+#include "axom/config.hpp"
+#include "axom/core/Macros.hpp"
+#include "axom/core/MDMapping.hpp"
+#include "axom/core/memory_management.hpp"
+#include "axom/core/utilities/Utilities.hpp"
+#include "axom/core/Types.hpp"
 #include "axom/core/StackArray.hpp"
-#include "axom/core/numerics/matvecops.hpp"  // for dot_product
-#include "axom/core/execution/for_all.hpp"   // for for_all, *_EXEC
+#include "axom/core/numerics/matvecops.hpp"
+#include "axom/core/execution/for_all.hpp"
 
 // C/C++ includes
 #include <iostream>  // for std::cerr and std::ostream
@@ -562,6 +562,21 @@ class ArrayBase<T, 1, ArrayType>
 private:
   constexpr static bool is_array_view = detail::ArrayTraits<ArrayType>::is_view;
 
+  /*!
+   * \brief Empty stand-in for the stride of an owning 1D Array.
+   *
+   * Owning 1D arrays are always contiguous, so their stride is the compile-time constant 1.
+   * Encoding that in the type (rather than storing a runtime int that is invariantly 1)
+   * lets operator[] and flatIndex() compile down to data()[idx], with no runtime multiply on the
+   * address-generation path. ArrayViews support runtime spacing and continue to store their stride as an int.
+   */
+  struct UnitStrideTag
+  {
+    AXOM_HOST_DEVICE constexpr UnitStrideTag(int = 1) { }
+    AXOM_HOST_DEVICE constexpr operator int() const { return 1; }
+  };
+  using StrideStorage = typename std::conditional<is_array_view, int, UnitStrideTag>::type;
+
 public:
   /* If ArrayType is an ArrayView, we use shallow-const semantics, akin to
    * std::span; a const ArrayView will still allow for mutating the underlying
@@ -574,25 +589,42 @@ class ArrayBase<T, 1, ArrayType>
 
   AXOM_HOST_DEVICE ArrayBase(IndexType = 0) { }
 
-  AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&, int stride = 1) : m_stride(stride) { }
+  AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&, int stride = 1) : m_stride(stride)
+  {
+    assert(is_array_view || stride == 1);
+  }
 
   AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&, const StackArray<IndexType, 1>& stride)
     : m_stride(static_cast<int>(stride[0]))
-  { }
+  {
+    assert(is_array_view || stride[0] == 1);
+  }
 
   AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&, const MDMapping<1>& mapping)
-    : m_stride(mapping.strides()[0])
-  { }
+    : m_stride(static_cast<int>(mapping.strides()[0]))
+  {
+    assert(is_array_view || mapping.strides()[0] == 1);
+  }
 
-  // Empty implementation because no member data
+  /*!
+   * \brief Copy the stride from another 1D array-like object.
+   *
+   * When this array is a view, the source's spacing must be preserved so the
+   * view continues to address the same elements. When this array is owning,
+   * the stride is the compile-time unit stride regardless of the source
+   * (a deep copy from a strided view compacts into contiguous storage).
+   */
   template <typename OtherArrayType>
-  AXOM_HOST_DEVICE ArrayBase(const ArrayBase<typename std::remove_const<T>::type, 1, OtherArrayType>&)
+  AXOM_HOST_DEVICE ArrayBase(
+    const ArrayBase<typename std::remove_const<T>::type, 1, OtherArrayType>& other)
+    : m_stride(static_cast<int>(other.minStride()))
   { }
 
-  // Empty implementation because no member data
+  /// \overload
   template <typename OtherArrayType>
   AXOM_HOST_DEVICE ArrayBase(
-    const ArrayBase<const typename std::remove_const<T>::type, 1, OtherArrayType>&)
+    const ArrayBase<const typename std::remove_const<T>::type, 1, OtherArrayType>& other)
+    : m_stride(static_cast<int>(other.minStride()))
   { }
 
   /// \brief Returns the dimensions of the Array
@@ -606,7 +638,7 @@ class ArrayBase<T, 1, ArrayType>
   /*!
    * \brief Returns the stride between adjacent items.
    */
-  AXOM_HOST_DEVICE IndexType minStride() const { return m_stride; }
+  AXOM_HOST_DEVICE constexpr IndexType minStride() const { return m_stride; }
 
   /*!
    * \brief Accessor, returns a reference to the given value.
@@ -655,8 +687,8 @@ class ArrayBase<T, 1, ArrayType>
   /// @}
 
   /// \brief Swaps two ArrayBases
-  /// No member data, so this is a no-op
-  void swap(ArrayBase&) { }
+  /// Swaps the stride; this is a no-op for owning arrays (unit stride).
+  void swap(ArrayBase& other) { std::swap(m_stride, other.m_stride); }
 
   /// \brief Set the shape
   /// No member data, so this is a no-op
@@ -666,7 +698,7 @@ class ArrayBase<T, 1, ArrayType>
   /*!
    * \brief Returns the minimum "chunk size" that should be allocated
    */
-  IndexType blockSize() const { return m_stride; }
+  constexpr IndexType blockSize() const { return m_stride; }
 
   /*!
    * \brief Updates the internal dimensions and striding based on the insertion
@@ -700,7 +732,7 @@ class ArrayBase<T, 1, ArrayType>
   }
   /// @}
 
-  int m_stride {1};
+  StrideStorage m_stride {1};
 };
 
 //------------------------------------------------------------------------------
@@ -1073,6 +1105,67 @@ struct ArrayOps
     }
   }
 
+  /*!
+   * \brief Fills an uninitialized array with a strided range of objects of type T.
+   *
+   * Similar to fill_range, but handles source data with non-unit stride. When src_stride == 1,
+   * this behaves identically to fill_range (and uses the same fast path). When src_stride > 1,
+   * elements are copied from positions 0, src_stride, 2*src_stride, etc.
+   *
+   * \param [inout] array the array to fill
+   * \param [in] begin the index at which to begin placing elements
+   * \param [in] nelems the number of elements to copy
+   * \param [in] values pointer to the first element of the source data
+   * \param [in] src_stride spacing between consecutive source elements
+   * \param [in] valueSpace the memory space in which values resides
+   */
+  void fill_range_strided(T* array,
+                          IndexType begin,
+                          IndexType nelems,
+                          const T* values,
+                          IndexType src_stride,
+                          MemorySpace valueSpace)
+  {
+    if constexpr(std::is_trivially_copyable_v<T>)
+    {
+      if(src_stride == 1)
+      {
+        // Contiguous case - use efficient bulk copy
+        axom::copy(array + begin, values, sizeof(T) * nelems);
+      }
+      else
+      {
+        // Strided case - element-by-element copy
+        StagingBuffer dst_buf(space, array, begin, nelems);
+        DeviceStagingBuffer<T> src_buf(valueSpace, const_cast<T*>(values), 0, nelems * src_stride, true);
+
+        T* dst = dst_buf.getStagingBuffer();
+        const T* src = src_buf.getStagingBuffer();
+
+        for(IndexType i = 0; i < nelems; ++i)
+        {
+          dst[i] = src[i * src_stride];
+        }
+        // Staging buffers clean up automatically via destructors
+      }
+    }
+    else
+    {
+      // Non-trivially copyable - use placement new with stride
+      StagingBuffer dst_buf(space, array, begin, nelems);
+      DeviceStagingBuffer<T> src_buf(valueSpace, const_cast<T*>(values), 0, nelems * src_stride, true);
+
+      T* dst = dst_buf.getStagingBuffer();
+      const T* src = src_buf.getStagingBuffer();
+
+      for(IndexType i = 0; i < nelems; ++i)
+      {
+        new(&dst[i]) T(src[i * src_stride]);
+      }
+      // Staging buffers clean up automatically via destructors
+    }
+  }
+
   /*!
    * \brief Constructs a new element in uninitialized memory.
    *