Redefine is_contiguous and is_strided_1d and introduce has_positive_s…

…trides in idx_map

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.

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;
+      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]]));
     }
 
     /**

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_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<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());
+}