Skip to content

Commit

Permalink
Merge pull request #1116 from LLNL/feature/yang39/array-view-stride
Browse files Browse the repository at this point in the history 
ArrayView: multi-dim subspans, non-row-major striding
  • Loading branch information
@publixsubfan
publixsubfan committed Jun 12, 2023
2 parents 20a033a + 8cdc947 commit 43df7b6
Show file tree
Hide file tree
Showing 9 changed files with 588 additions and 62 deletions.
5 changes: 5 additions & 0 deletions RELEASE-NOTES.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,10 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
- Quest: Adds ability to import volume fractions into `SamplingShaper` before processing `Klee` input
- Slam: adds a `slam::MappedVariableCardinality` policy to accelerate mapping flat indices
back to first-set indices when used in a `StaticRelation`
- Adds an `ArrayView(data, shape, strides)` constructor to support column-major and custom
striding layouts.
- Adds an `ArrayView::subspan()` overload for multi-dimensional subspans
- Adds an `axom::utilities::insertionSort()` method.

### Changed
- Fixed bug in `mint::mesh::UnstructuredMesh` constructors, affecting capacity.
Expand Down Expand Up @@ -103,6 +107,7 @@ The Axom project release numbers follow [Semantic Versioning](http://semver.org/
- Multimat: Ported field data/sparsity layout conversion methods to GPU.
- Multimat: `MultiMat::makeOtherRelation()` now runs on the GPU with an appropriately-set allocator ID.
- Multimat: `MultiMat::setCellMatRel(counts, indices)` now runs on the GPU, and accepts GPU-side data.
- Renames `ArrayView::spacing()` to `ArrayView::minStride()`.

### Fixed
- Fixed issues with CUDA build in CMake versions 3.14.5 and above. Now require CMake 3.18+
Expand Down
5 changes: 0 additions & 5 deletions src/axom/core/Array.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -344,11 +344,6 @@ class Array : public ArrayBase<T, DIM, Array<T, DIM, SPACE>>

/// @}

/*!
\brief Returns spacing between adjacent items.
*/
AXOM_HOST_DEVICE IndexType spacing() const { return 1; }

/// @}

/// \name Array methods to modify the data.
Expand Down
179 changes: 140 additions & 39 deletions src/axom/core/ArrayBase.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,6 @@ bool operator!=(const ArrayBase<T1, DIM, LArrayType>& lhs,
* IndexType size() const;
* T* data();
* const T* data() const;
* IndexType spacing() const;
* int getAllocatorID() const;
* \endcode
*
Expand Down Expand Up @@ -167,14 +166,31 @@ class ArrayBase
* \brief Parameterized constructor that sets up the array shape.
*
* \param [in] shape Array size in each direction.
* \param [in] min_stride The minimum stride between two consecutive
* elements in row-major order.
*/
AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, DIM>& shape)
AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, DIM>& shape,
int min_stride = 1)
: m_shape {shape}
{
m_strides[DIM - 1] = 1;
m_strides[DIM - 1] = min_stride;
updateStrides();
}

/*!
* \brief Parameterized constructor that sets up the array shape and stride.
*
* \param [in] shape Array size in each direction.
* \param [in] stride Array stride for each direction.
*/
AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, DIM>& shape,
const StackArray<IndexType, DIM>& stride)
: m_shape {shape}
, m_strides {stride}
{
validateShapeAndStride(shape, stride);
}

/*!
* \brief Copy constructor for arrays of different type
* Because the element type (T) and dimension (DIM) are still locked down,
Expand Down Expand Up @@ -288,20 +304,20 @@ class ArrayBase
*
* \param [in] idx the position of the value to return.
*
* \note equivalent to *(array.data() + idx * spacing()).
* \note equivalent to *(array.data() + idx * minStride()).
*
* \pre 0 <= idx < asDerived().size()
*/
AXOM_HOST_DEVICE T& flatIndex(const IndexType idx)
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * asDerived().minStride()];
}
/// \overload
AXOM_HOST_DEVICE RealConstT& flatIndex(const IndexType idx) const
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * asDerived().minStride()];
}
/// @}

Expand All @@ -312,6 +328,34 @@ class ArrayBase
std::swap(m_strides, other.m_strides);
}

/// \brief Returns the dimensions of the Array
AXOM_HOST_DEVICE const StackArray<IndexType, DIM>& shape() const
{
return m_shape;
}

/*!
* \brief Returns the memory strides of the Array.
*/
AXOM_HOST_DEVICE const StackArray<IndexType, DIM>& strides() const
{
return m_strides;
}

/*!
* \brief Returns the minimum stride between adjacent items.
*/
AXOM_HOST_DEVICE IndexType minStride() const
{
IndexType minStride = m_strides[0];
for(int dim = 1; dim < DIM; dim++)
{
minStride = axom::utilities::min(minStride, m_strides[dim]);
}
return minStride;
}

protected:
/// \brief Set the shape
AXOM_HOST_DEVICE void setShape(const StackArray<IndexType, DIM>& shape_)
{
Expand All @@ -326,23 +370,17 @@ class ArrayBase
updateStrides();
}

/// \brief Returns the dimensions of the Array
AXOM_HOST_DEVICE const StackArray<IndexType, DIM>& shape() const
{
return m_shape;
}

/*!
\brief Returns the logical strides of the Array.
Note: Memory stride is logical stride times spacing.
*/
AXOM_HOST_DEVICE const StackArray<IndexType, DIM>& strides() const
/// \brief Set the shape and stride
AXOM_HOST_DEVICE void setShapeAndStride(const StackArray<IndexType, DIM>& shape,
const StackArray<IndexType, DIM>& stride)
{
return m_strides;
#ifdef AXOM_DEBUG
validateShapeAndStride(shape, stride);
#endif
m_shape = shape;
m_strides = stride;
}

protected:
/*!
* \brief Returns the minimum "chunk size" that should be allocated
* For example, 2 would be the chunk size of a 2D array whose second dimension is of size 2.
Expand Down Expand Up @@ -398,17 +436,29 @@ class ArrayBase
return static_cast<const ArrayType&>(*this);
}

/*!
\brief Logical offset to get to the given multidimensional index.
To get from logical offset to memory offset, multiply by spacing().
*/
//// \brief Memory offset to get to the given multidimensional index.
AXOM_HOST_DEVICE IndexType offset(const StackArray<IndexType, DIM>& idx) const
{
return numerics::dot_product((const IndexType*)idx, m_strides.begin(), DIM);
}

/// Logical offset to a slice at the given lower-dimensional index.
/*!
* \brief Returns the size of the range of memory in which elements are
* located. This is equivalent to size() * minStride().
*
* offset() will return a value between [0, memorySize()).
*/
AXOM_HOST_DEVICE IndexType memorySize() const
{
IndexType maxSize = 0;
for(int dim = 0; dim < DIM; dim++)
{
maxSize = axom::utilities::max(maxSize, m_strides[dim] * m_shape[dim]);
}
return maxSize;
}

/// \brief Memory offset to a slice at the given lower-dimensional index.
template <int UDim>
AXOM_HOST_DEVICE IndexType offset(const StackArray<IndexType, UDim>& idx) const
{
Expand All @@ -423,8 +473,35 @@ class ArrayBase
/*! \brief Test if idx is within bounds */
AXOM_HOST_DEVICE inline bool inBounds(IndexType idx) const
{
return idx >= 0 && idx < asDerived().size();
return idx >= 0 && idx < memorySize();
}

/// \brief Checks a shape and stride array for correct bounds.
AXOM_HOST_DEVICE inline void validateShapeAndStride(
const StackArray<IndexType, DIM>& shape,
const StackArray<IndexType, DIM>& stride)
{
int sorted_dims[DIM];
for(int dim = 0; dim < DIM; dim++)
{
sorted_dims[dim] = dim;
}
// Sort the dimensions by stride.
axom::utilities::insertionSort(sorted_dims,
DIM,
[&](int dim_a, int dim_b) -> bool {
return stride[dim_a] < stride[dim_b];
});
// Work from the smallest-strided dimension to the largest-strided.
for(int dim = 0; dim < DIM - 1; dim++)
{
int minor_dim = sorted_dims[dim];
int major_dim = sorted_dims[dim + 1];
assert(stride[major_dim] >= stride[minor_dim] * shape[minor_dim]);
assert(stride[major_dim] % stride[minor_dim] == 0);
}
}

/// @}

/// \name Internal subarray slicing methods
Expand All @@ -450,7 +527,7 @@ class ArrayBase
{
const IndexType baseIdx = offset(idx);
assert(inBounds(baseIdx));
return asDerived().data()[baseIdx * asDerived().spacing()];
return asDerived().data()[baseIdx];
}

/// \overload
Expand All @@ -459,7 +536,7 @@ class ArrayBase
{
const IndexType baseIdx = offset(idx);
assert(inBounds(baseIdx));
return asDerived().data()[baseIdx * asDerived().spacing()];
return asDerived().data()[baseIdx];
}
/// @}

Expand Down Expand Up @@ -489,7 +566,14 @@ class ArrayBase<T, 1, ArrayType>

AXOM_HOST_DEVICE ArrayBase(IndexType = 0) { }

AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&) { }
AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&, int stride = 1)
: m_stride(stride)
{ }

AXOM_HOST_DEVICE ArrayBase(const StackArray<IndexType, 1>&,
const StackArray<IndexType, 1>& stride)
: m_stride(stride[0])
{ }

// Empty implementation because no member data
template <typename OtherArrayType>
Expand All @@ -509,6 +593,11 @@ class ArrayBase<T, 1, ArrayType>
return {{asDerived().size()}};
}

/*!
* \brief Returns the stride between adjacent items.
*/
AXOM_HOST_DEVICE IndexType minStride() const { return m_stride; }

/*!
* \brief Accessor, returns a reference to the given value.
* For multidimensional arrays, indexes into the (flat) raw data.
Expand All @@ -523,13 +612,13 @@ class ArrayBase<T, 1, ArrayType>
AXOM_HOST_DEVICE T& operator[](const IndexType idx)
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * m_stride];
}
/// \overload
AXOM_HOST_DEVICE RealConstT& operator[](const IndexType idx) const
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * m_stride];
}

/*!
Expand All @@ -538,20 +627,20 @@ class ArrayBase<T, 1, ArrayType>
*
* \param [in] idx the position of the value to return.
*
* \note equivalent to *(array.data() + idx * array.spacing()).
* \note equivalent to *(array.data() + idx * array.minStride()).
*
* \pre 0 <= idx < asDerived().size()
*/
AXOM_HOST_DEVICE T& flatIndex(const IndexType idx)
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * m_stride];
}
/// \overload
AXOM_HOST_DEVICE RealConstT& flatIndex(const IndexType idx) const
{
assert(inBounds(idx));
return asDerived().data()[idx * asDerived().spacing()];
return asDerived().data()[idx * m_stride];
}
/// @}

Expand All @@ -567,7 +656,7 @@ class ArrayBase<T, 1, ArrayType>
/*!
* \brief Returns the minimum "chunk size" that should be allocated
*/
IndexType blockSize() const { return 1; }
IndexType blockSize() const { return m_stride; }

/*!
* \brief Updates the internal dimensions and striding based on the insertion
Expand All @@ -591,12 +680,23 @@ class ArrayBase<T, 1, ArrayType>
/// \name Internal bounds-checking routines
/// @{

/*!
* \brief Returns the size of the range of memory in which elements are
* located. This is equivalent to size() * minStride().
*/
AXOM_HOST_DEVICE IndexType memorySize() const
{
return m_stride * shape()[0];
}

/*! \brief Test if idx is within bounds */
AXOM_HOST_DEVICE inline bool inBounds(IndexType idx) const
{
return idx >= 0 && idx < asDerived().size();
return idx >= 0 && idx < memorySize();
}
/// @}

int m_stride {1};
};

//------------------------------------------------------------------------------
Expand Down Expand Up @@ -1388,7 +1488,8 @@ class ArraySubslice
public:
AXOM_HOST_DEVICE ArraySubslice(BaseArray* array,
const StackArray<IndexType, NumIndices>& idxs)
: BaseClass(detail::takeLastElems<SliceDim>(array->shape()))
: BaseClass(detail::takeLastElems<SliceDim>(array->shape()),
detail::takeLastElems<SliceDim>(array->strides()))
, m_array(array)
, m_inds(idxs)
{ }
Expand Down Expand Up @@ -1416,7 +1517,7 @@ class ArraySubslice
IndexType offset = numerics::dot_product(m_inds.begin(),
m_array->strides().begin(),
NumIndices);
return m_array->data() + offset * m_array->spacing();
return m_array->data() + offset;
}

/// @}
Expand Down

0 comments on commit 43df7b6

Please sign in to comment.