diff --git a/c++/nda/_impl_basic_array_view_common.hpp b/c++/nda/_impl_basic_array_view_common.hpp index 12ee109d3..b2a7c02cc 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. @@ -516,7 +522,7 @@ void fill_with_scalar(Scalar const &scalar) noexcept { } else { const long stri = indexmap().min_stride(); const long Lstri = L * stri; - for (long i = 0; i < Lstri; i += stri) p[i] = scalar; + for (long i = 0; i != Lstri; i += stri) p[i] = scalar; } } else { // no compile-time memory layout guarantees diff --git a/c++/nda/layout/idx_map.hpp b/c++/nda/layout/idx_map.hpp index 683a066d9..1e297b28c 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; + while (len[stride_order[i]] == 1 and i > 0) --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 12332631e..c7284f2b3 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 c2b7393e6..f0bfbc79b 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 6b3b6c99c..98ad79e9f 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 35cbe05b4..cc2697b5a 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 c509f0b25..8f8e126c6 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_basic_array_and_view.cpp b/test/c++/nda_basic_array_and_view.cpp index a312a7397..43f8adb93 100644 --- a/test/c++/nda_basic_array_and_view.cpp +++ b/test/c++/nda_basic_array_and_view.cpp @@ -553,6 +553,28 @@ TEST_F(NDAArrayAndView, SliceAccessFull) { auto A_s2 = A_3d(nda::range::all, nda::range::all, nda::range::all); EXPECT_EQ(A_s2.shape(), shape_3d); EXPECT_EQ(A_s2, A_3d); + + // full slice with negative strides + auto A_s3 = A_3d(nda::range(1, -1, -1), nda::range(2, -1, -1), nda::range(3, -1, -1)); + EXPECT_EQ(A_s3.shape(), shape_3d); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 3; ++j) { + for (long k = 0; k < 4; ++k) EXPECT_EQ(A_s3(i, j, k), A_3d(1 - i, 2 - j, 3 - k)); + } + } + A_s3 = 42; + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 3; ++j) { + for (long k = 0; k < 4; ++k) EXPECT_EQ(A_3d(i, j, k), 42); + } + } + auto A_3d_copy = A_3d; + A_s3 = A_3d_copy; + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 3; ++j) { + for (long k = 0; k < 4; ++k) EXPECT_EQ(A_3d(i, j, k), A_3d_copy(1 - i, 2 - j, 3 - k)); + } + } } TEST_F(NDAArrayAndView, SliceAccess1D) { @@ -579,6 +601,32 @@ TEST_F(NDAArrayAndView, SliceAccess1D) { for (long i = 0; i < 2; ++i) EXPECT_EQ(C_s(i), A_3d(0, 1, 1 + 2 * i)); C_s = 42; EXPECT_NE(C, A_3d); + + // full 1D slice with negative strides + auto D = A_3d; + auto D_s = D(0, nda::range(2, -1, -1), 2); + static_assert(D_s.rank == 1); + EXPECT_EQ(D_s.size(), 3); + EXPECT_EQ(D_s.strides(), (std::array{-4})); + for (long i = 0; i < 3; ++i) EXPECT_EQ(D_s(i), A_3d(0, 2 - i, 2)); + D_s = 42; + EXPECT_NE(D, A_3d); + for (long i = 0; i < 3; ++i) EXPECT_EQ(D(0, i, 2), 42); + D_s = nda::array{1, 2, 3}; + for (long i = 0; i < 3; ++i) EXPECT_EQ(D(0, i, 2), 3 - i); + + // partial 1D slice with negative strides + auto E = A_3d; + auto E_s = E(0, 1, nda::range(3, 0, -2)); + static_assert(E_s.rank == 1); + EXPECT_EQ(E_s.size(), 2); + EXPECT_EQ(E_s.strides(), (std::array{-2})); + for (long i = 0; i < 2; ++i) EXPECT_EQ(E_s(i), A_3d(0, 1, 3 - 2 * i)); + E_s = 42; + EXPECT_NE(E, A_3d); + for (long i = 0; i < 2; ++i) EXPECT_EQ(E(0, 1, 1 + 2 * i), 42); + E_s = nda::array{1, 2}; + for (long i = 0; i < 2; ++i) EXPECT_EQ(E(0, 1, 1 + 2 * i), 2 - i); } TEST_F(NDAArrayAndView, SliceAccess2D) { @@ -634,6 +682,80 @@ TEST_F(NDAArrayAndView, SliceAccess2D) { static_assert(C_s3.rank == 2); EXPECT_EQ(C_s3.shape(), (std::array{2, 2})); EXPECT_EQ(C_s3.strides(), (std::array{12, 8})); + + // full 2D slice with negative strides + auto D = A_3d; + auto D_s = D(nda::range(1, -1, -1), 1, nda::range(3, -1, -1)); + static_assert(D_s.rank == 2); + EXPECT_EQ(D_s.shape(), (std::array{2, 4})); + EXPECT_EQ(D_s.strides(), (std::array{-12, -1})); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D_s(i, j), A_3d(1 - i, 1, 3 - j)); + } + D_s = 42; + EXPECT_NE(D, A_3d); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D(i, 1, j), 42); + } + D_s = A_3d(nda::range::all, 1, nda::range::all); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D(i, 1, j), A_3d(1 - i, 1, 3 - j)); + } + + D = A_3d; + auto D_s2 = D(1, nda::range(2, -1, -1), nda::range::all); + static_assert(D_s2.rank == 2); + EXPECT_EQ(D_s2.shape(), (std::array{3, 4})); + EXPECT_EQ(D_s2.strides(), (std::array{-4, 1})); + for (long i = 0; i < 3; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D_s2(i, j), A_3d(1, 2 - i, j)); + } + D_s2 = 42; + EXPECT_NE(D, A_3d); + for (long i = 0; i < 3; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D(1, i, j), 42); + } + D_s2 = A_3d(1, nda::range::all, nda::range::all); + for (long i = 0; i < 3; ++i) { + for (long j = 0; j < 4; ++j) EXPECT_EQ(D(1, i, j), A_3d(1, 2 - i, j)); + } + + // partial 2D slice with negative strides + auto E = A_3d; + auto E_s = E(0, nda::range(1, -1, -1), nda::range(3, 0, -2)); + static_assert(E_s.rank == 2); + EXPECT_EQ(E_s.shape(), (std::array{2, 2})); + EXPECT_EQ(E_s.strides(), (std::array{-4, -2})); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E_s(i, j), A_3d(0, 1 - i, 3 - 2 * j)); + } + E_s = 42; + EXPECT_NE(E, A_3d); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E(0, i, 1 + 2 * j), 42); + } + E_s = A_3d(0, nda::range(0, 2), nda::range(1, 4, 2)); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E(0, i, 1 + 2 * j), A_3d(0, 1 - i, 3 - 2 * j)); + } + + E = A_3d; + auto E_s2 = E(nda::range::all, nda::range(2, -1, -2), 3); + static_assert(E_s2.rank == 2); + EXPECT_EQ(E_s2.shape(), (std::array{2, 2})); + EXPECT_EQ(E_s2.strides(), (std::array{12, -8})); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E_s2(i, j), A_3d(i, 2 - 2 * j, 3)); + } + E_s2 = 42; + EXPECT_NE(E, A_3d); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E(i, 2 * j, 3), 42); + } + E_s2 = A_3d(nda::range::all, nda::range(0, 3, 2), 3); + for (long i = 0; i < 2; ++i) { + for (long j = 0; j < 2; ++j) EXPECT_EQ(E(i, 2 * j, 3), A_3d(i, 2 - 2 * j, 3)); + } } TEST_F(NDAArrayAndView, SliceAccess3D) { @@ -844,8 +966,6 @@ TEST_F(NDAArrayAndView, StrideOrderOfArrays) { #if defined(__has_feature) #if !__has_feature(address_sanitizer) -TEST_F(NDAArrayAndView, BadAlloc) { - EXPECT_THROW(nda::vector(long(1e16)), std::bad_alloc); -} +TEST_F(NDAArrayAndView, BadAlloc) { EXPECT_THROW(nda::vector(long(1e16)), std::bad_alloc); } #endif #endif diff --git a/test/c++/nda_bugs_and_issues.cpp b/test/c++/nda_bugs_and_issues.cpp index 77e95baf2..0352f040f 100644 --- a/test/c++/nda_bugs_and_issues.cpp +++ b/test/c++/nda_bugs_and_issues.cpp @@ -124,3 +124,48 @@ TEST(NDA, AmbiguousAssignmentOperatorIssue) { B = std::array{1, 2, 3}; for (int i = 0; i < 3; ++i) { EXPECT_EQ(B(i), (std::array{1, 2, 3})); } } + +TEST(NDA, AssignmentTo1DNegativeStridedViewsIssue) { + // issue concerning the assignment of scalar values to 1D strided views with negative strides + nda::array A(10); + for (int i = 0; auto &x : A) x = i++; + + auto B = A(nda::range(9, -1, -1)); + for (int i = 9; auto x : B) EXPECT_EQ(x, i--); + + 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()); +}