Fixes for negative strides

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

c++/nda/layout/idx_map.hpp

@@ -184,7 +184,7 @@ namespace nda {
     [[nodiscard]] std::array<long, Rank> 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]]));
     }
 
     /**

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<A>;

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);

c++/nda/mpi/gather.hpp

@@ -97,7 +97,9 @@ struct mpi::lazy<mpi::tag::gather, A> {
   template <nda::Array T>
   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<A>::layout_t::stride_order_encoded == std::decay_t<T>::layout_t::stride_order_encoded,
                   "Error in MPI gather for nda::Array: Incompatible stride orders");
 
@@ -166,7 +168,8 @@ namespace nda {
   ArrayInitializer<std::remove_reference_t<A>> auto mpi_gather(A &&a, mpi::communicator comm = {}, int root = 0, bool all = false)
     requires(is_regular_or_view_v<A>)
   {
-    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<mpi::tag::gather, A>{std::forward<A>(a), comm, root, all};
   }
 

c++/nda/mpi/reduce.hpp

@@ -88,7 +88,9 @@ struct mpi::lazy<mpi::tag::reduce, A> {
   template <nda::Array T>
   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<A>::layout_t::stride_order_encoded == std::decay_t<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<A>)
   {
-    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<mpi::tag::reduce, A>{std::forward<A>(a), comm, root, all, op};
   }
 

c++/nda/mpi/scatter.hpp

@@ -91,7 +91,9 @@ struct mpi::lazy<mpi::tag::scatter, A> {
    */
   template <nda::Array T>
   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<A>::layout_t::stride_order_encoded == std::decay_t<T>::layout_t::stride_order_encoded,
                   "Error in MPI scatter for nda::Array: Incompatible stride orders");
 
@@ -156,7 +158,8 @@ namespace nda {
   ArrayInitializer<std::remove_reference_t<A>> auto mpi_scatter(A &&a, mpi::communicator comm = {}, int root = 0, bool all = false)
     requires(is_regular_or_view_v<A>)
   {
-    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<mpi::tag::scatter, A>{std::forward<A>(a), comm, root, all};
   }
 

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<long, 1>{-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<int, 1>{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<long, 1>{-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<int, 1>{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<long, 2>{2, 2}));
   EXPECT_EQ(C_s3.strides(), (std::array<long, 2>{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<long, 2>{2, 4}));
+  EXPECT_EQ(D_s.strides(), (std::array<long, 2>{-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<long, 2>{3, 4}));
+  EXPECT_EQ(D_s2.strides(), (std::array<long, 2>{-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<long, 2>{2, 2}));
+  EXPECT_EQ(E_s.strides(), (std::array<long, 2>{-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<long, 2>{2, 2}));
+  EXPECT_EQ(E_s2.strides(), (std::array<long, 2>{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<int>(long(1e16)), std::bad_alloc);
-}
+TEST_F(NDAArrayAndView, BadAlloc) { EXPECT_THROW(nda::vector<int>(long(1e16)), std::bad_alloc); }
 #endif
 #endif

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<int, 1> 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<int, 2> 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<int, 2> 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());
+}