From 86b33908bec895ffedf9cb5668883f33db33722b Mon Sep 17 00:00:00 2001 From: Thomas Hahn Date: Wed, 26 Jun 2024 16:15:19 -0400 Subject: [PATCH] Redefine is_contiguous and is_strided_1d and introduce has_positive_strides in idx_map --- c++/nda/_impl_basic_array_view_common.hpp | 6 +++++ c++/nda/layout/idx_map.hpp | 26 +++++++++++------- c++/nda/layout_transforms.hpp | 3 ++- c++/nda/linalg/eigenelements.hpp | 3 ++- c++/nda/mpi/gather.hpp | 7 +++-- c++/nda/mpi/reduce.hpp | 7 +++-- c++/nda/mpi/scatter.hpp | 7 +++-- test/c++/nda_bugs_and_issues.cpp | 33 +++++++++++++++++++++++ 8 files changed, 75 insertions(+), 17 deletions(-) diff --git a/c++/nda/_impl_basic_array_view_common.hpp b/c++/nda/_impl_basic_array_view_common.hpp index 9346f3c9..b2a7c02c 100644 --- a/c++/nda/_impl_basic_array_view_common.hpp +++ b/c++/nda/_impl_basic_array_view_common.hpp @@ -82,6 +82,12 @@ */ [[nodiscard]] long is_contiguous() const noexcept { return lay.is_contiguous(); } +/** + * @brief Are all the strides of the memory layout of the view/array positive? + * @return True if the nda::idx_map has positive strides, false otherwise. + */ +[[nodiscard]] long has_positive_strides() const noexcept { return lay.has_positive_strides(); } + /** * @brief Is the view/array empty? * @return True if the view/array does not contain any elements. diff --git a/c++/nda/layout/idx_map.hpp b/c++/nda/layout/idx_map.hpp index 683a066d..d893b1cb 100644 --- a/c++/nda/layout/idx_map.hpp +++ b/c++/nda/layout/idx_map.hpp @@ -184,7 +184,7 @@ namespace nda { [[nodiscard]] std::array const &strides() const noexcept { return str; } /** - * @brief Get the value of the smallest stride. + * @brief Get the value of the smallest stride (positive or negative). * @return Stride of the fastest varying dimension. */ [[nodiscard]] long min_stride() const noexcept { return str[stride_order[Rank - 1]]; } @@ -193,31 +193,39 @@ namespace nda { * @brief Is the data contiguous in memory? * * @details The data is contiguous in memory if the size of the map is equal to the number of memory locations - * spanned by the map, i.e. the product of the largest stride times the length of the corresponding dimension. + * spanned by the map, i.e. the absolute value of the product of the largest stride times the length of the + * corresponding dimension. * * @return True if the data is contiguous in memory, false otherwise. */ [[nodiscard]] bool is_contiguous() const noexcept { auto s = size(); if (s == 0) return true; - auto const str_x_len = str * len; - return (*std::max_element(str_x_len.cbegin(), str_x_len.cend()) == s); + return (std::abs(str[stride_order[0]] * len[stride_order[0]]) == s); } + /** + * @brief Are all strides positive? + * @return True if all strides are positive, false otherwise. + */ + [[nodiscard]] bool has_positive_strides() const noexcept { return (*std::min_element(str.cbegin(), str.cend()) >= 0); } + /** * @brief Is the data strided in memory with a constant stride? * - * @details The data is strided in memory with a constant stride if the size of the map times the stride of the - * fastest dimension with an extent bigger than 1 is equal to the number of memory locations spanned by the map, - * i.e. the product of the largest stride times the length of the corresponding dimension. + * @details The data is strided in memory with a constant stride if the size of the map times the absolute value of + * the stride of the fastest dimension with an extent bigger than 1 is equal to the number of memory locations + * spanned by the map, i.e. the absolute value of the product of the largest stride times the length of the + * corresponding dimension. * * @return True if the data is strided in memory with a constant stride, false otherwise. */ [[nodiscard]] bool is_strided_1d() const noexcept { auto s = size(); if (s == 0) return true; - auto const str_x_len = str * len; - return (*std::max_element(str_x_len.cbegin(), str_x_len.cend()) == s * min_stride()); + int i = Rank - 1; + for (; len[stride_order[i]] == 1; --i); + return (std::abs(str[stride_order[0]] * len[stride_order[0]]) == s * std::abs(str[stride_order[i]])); } /** diff --git a/c++/nda/layout_transforms.hpp b/c++/nda/layout_transforms.hpp index 12332631..c7284f2b 100644 --- a/c++/nda/layout_transforms.hpp +++ b/c++/nda/layout_transforms.hpp @@ -109,7 +109,8 @@ namespace nda { // check size and contiguity of new shape EXPECTS_WITH_MESSAGE(a.size() == (std::accumulate(new_shape.cbegin(), new_shape.cend(), Int{1}, std::multiplies<>{})), "Error in nda::reshape: New shape has an incorrect number of elements"); - EXPECTS_WITH_MESSAGE(a.indexmap().is_contiguous(), "Error in nda::reshape: Only contiguous arrays/views are supported"); + EXPECTS_WITH_MESSAGE(a.is_contiguous(), "Error in nda::reshape: Only contiguous arrays/views are supported"); + EXPECTS_WITH_MESSAGE(a.has_positive_strides(), "Error in nda::reshape: Only arrays/views with positive strides are supported") // restrict supported layouts (why?) using A_t = std::remove_cvref_t; diff --git a/c++/nda/linalg/eigenelements.hpp b/c++/nda/linalg/eigenelements.hpp index c2b7393e..f0bfbc79 100644 --- a/c++/nda/linalg/eigenelements.hpp +++ b/c++/nda/linalg/eigenelements.hpp @@ -45,7 +45,8 @@ namespace nda::linalg { // runtime checks EXPECTS((not m.empty())); EXPECTS(is_matrix_square(m, true)); - EXPECTS(m.indexmap().is_contiguous()); + EXPECTS(m.is_contiguous()); + EXPECTS(m.has_positive_strides()); // set up the workspace int dim = m.extent(0); diff --git a/c++/nda/mpi/gather.hpp b/c++/nda/mpi/gather.hpp index 6b3b6c99..98ad79e9 100644 --- a/c++/nda/mpi/gather.hpp +++ b/c++/nda/mpi/gather.hpp @@ -97,7 +97,9 @@ struct mpi::lazy { template void invoke(T &&target) const { // NOLINT (temporary views are allowed here) // check if the arrays can be used in the MPI call - if (not target.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI gather for nda::Array: Target array needs to be contiguous"; + if (not target.is_contiguous() or not target.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI gather for nda::Array: Target array needs to be contiguous with positive strides"; + static_assert(std::decay_t::layout_t::stride_order_encoded == std::decay_t::layout_t::stride_order_encoded, "Error in MPI gather for nda::Array: Incompatible stride orders"); @@ -166,7 +168,8 @@ namespace nda { ArrayInitializer> auto mpi_gather(A &&a, mpi::communicator comm = {}, int root = 0, bool all = false) requires(is_regular_or_view_v) { - if (not a.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI gather for nda::Array: Array needs to be contiguous"; + if (not a.is_contiguous() or not a.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI gather for nda::Array: Array needs to be contiguous with positive strides"; return mpi::lazy{std::forward(a), comm, root, all}; } diff --git a/c++/nda/mpi/reduce.hpp b/c++/nda/mpi/reduce.hpp index 35cbe05b..cc2697b5 100644 --- a/c++/nda/mpi/reduce.hpp +++ b/c++/nda/mpi/reduce.hpp @@ -88,7 +88,9 @@ struct mpi::lazy { template void invoke(T &&target) const { // NOLINT (temporary views are allowed here) // check if the arrays can be used in the MPI call - if (not target.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI reduce for nda::Array: Target array needs to be contiguous"; + if (not target.is_contiguous() or not target.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI reduce for nda::Array: Target array needs to be contiguous with positive strides"; + static_assert(std::decay_t::layout_t::stride_order_encoded == std::decay_t::layout_t::stride_order_encoded, "Error in MPI reduce for nda::Array: Incompatible stride orders"); @@ -166,7 +168,8 @@ namespace nda { MPI_Op op = MPI_SUM) requires(is_regular_or_view_v) { - if (not a.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI reduce for nda::Array: Array needs to be contiguous"; + if (not a.is_contiguous() or not a.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI reduce for nda::Array: Array needs to be contiguous with positive strides"; return mpi::lazy{std::forward(a), comm, root, all, op}; } diff --git a/c++/nda/mpi/scatter.hpp b/c++/nda/mpi/scatter.hpp index c509f0b2..8f8e126c 100644 --- a/c++/nda/mpi/scatter.hpp +++ b/c++/nda/mpi/scatter.hpp @@ -91,7 +91,9 @@ struct mpi::lazy { */ template void invoke(T &&target) const { // NOLINT (temporary views are allowed here) - if (not target.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Target array needs to be contiguous"; + if (not target.is_contiguous() or not target.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Target array needs to be contiguous with positive strides"; + static_assert(std::decay_t::layout_t::stride_order_encoded == std::decay_t::layout_t::stride_order_encoded, "Error in MPI scatter for nda::Array: Incompatible stride orders"); @@ -156,7 +158,8 @@ namespace nda { ArrayInitializer> auto mpi_scatter(A &&a, mpi::communicator comm = {}, int root = 0, bool all = false) requires(is_regular_or_view_v) { - if (not a.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Array needs to be contiguous"; + if (not a.is_contiguous() or not a.has_positive_strides()) + NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Array needs to be contiguous with positive strides"; return mpi::lazy{std::forward(a), comm, root, all}; } diff --git a/test/c++/nda_bugs_and_issues.cpp b/test/c++/nda_bugs_and_issues.cpp index c82eb4ab..0352f040 100644 --- a/test/c++/nda_bugs_and_issues.cpp +++ b/test/c++/nda_bugs_and_issues.cpp @@ -136,3 +136,36 @@ TEST(NDA, AssignmentTo1DNegativeStridedViewsIssue) { B = 1; for (auto x : B) EXPECT_EQ(x, 1); } + +TEST(NDA, IsStrided1DAndIsContiguousIssue) { + // issue concerning the nda::idx_map::is_strided_1d and the nda::idx_map::is_contiguous function + nda::array A(10, 10); + EXPECT_TRUE(A.indexmap().is_contiguous()); + EXPECT_TRUE(A.indexmap().has_positive_strides()); + EXPECT_TRUE(A.indexmap().is_strided_1d()); + + auto A_v1 = A(nda::range::all, nda::range(9, -1, -1)); + EXPECT_TRUE(A_v1.indexmap().is_contiguous()); + EXPECT_FALSE(A_v1.indexmap().has_positive_strides()); + EXPECT_TRUE(A_v1.indexmap().is_strided_1d()); + + auto A_v2 = A(nda::range(9, -1, -1), nda::range::all); + EXPECT_TRUE(A_v2.indexmap().is_contiguous()); + EXPECT_FALSE(A_v2.indexmap().has_positive_strides()); + EXPECT_TRUE(A_v2.indexmap().is_strided_1d()); + + auto A_v3 = A(nda::range(9, -1, -1), nda::range(9, -1, -1)); + EXPECT_TRUE(A_v3.indexmap().is_contiguous()); + EXPECT_FALSE(A_v3.indexmap().has_positive_strides()); + EXPECT_TRUE(A_v3.indexmap().is_strided_1d()); + + nda::array B(10, 1); + EXPECT_TRUE(B.indexmap().is_contiguous()); + EXPECT_TRUE(B.indexmap().has_positive_strides()); + EXPECT_TRUE(B.indexmap().is_strided_1d()); + + auto B_v1 = B(nda::range(0, 10, 2), nda::range::all); + EXPECT_FALSE(B_v1.indexmap().is_contiguous()); + EXPECT_TRUE(B_v1.indexmap().has_positive_strides()); + EXPECT_TRUE(B_v1.indexmap().is_strided_1d()); +}