Merge pull request #885 from SparrowLii/negative_strides

Fix memory continuity judgment when stride is negative

Author: bluss

blas-tests/tests/oper.rs

@@ -173,7 +173,7 @@ where
     S2: Data<Elem = A>,
 {
     let ((m, _), k) = (lhs.dim(), rhs.dim());
-    reference_mat_mul(lhs, &rhs.to_owned().into_shape((k, 1)).unwrap())
+    reference_mat_mul(lhs, &rhs.as_standard_layout().into_shape((k, 1)).unwrap())
         .into_shape(m)
         .unwrap()
 }
@@ -186,7 +186,7 @@ where
     S2: Data<Elem = A>,
 {
     let (m, (_, n)) = (lhs.dim(), rhs.dim());
-    reference_mat_mul(&lhs.to_owned().into_shape((1, m)).unwrap(), rhs)
+    reference_mat_mul(&lhs.as_standard_layout().into_shape((1, m)).unwrap(), rhs)
         .into_shape(n)
         .unwrap()
 }

src/dimension/dimension_trait.rs

@@ -287,17 +287,16 @@ pub trait Dimension:
             return true;
         }
         if dim.ndim() == 1 {
-            return false;
+            return strides[0] as isize == -1;
         }
         let order = strides._fastest_varying_stride_order();
         let strides = strides.slice();
 
-        // FIXME: Negative strides
         let dim_slice = dim.slice();
         let mut cstride = 1;
         for &i in order.slice() {
             // a dimension of length 1 can have unequal strides
-            if dim_slice[i] != 1 && strides[i] != cstride {
+            if dim_slice[i] != 1 && (strides[i] as isize).abs() as usize != cstride {
                 return false;
             }
             cstride *= dim_slice[i];
@@ -308,16 +307,17 @@ pub trait Dimension:
     /// Return the axis ordering corresponding to the fastest variation
     /// (in ascending order).
     ///
-    /// Assumes that no stride value appears twice. This cannot yield the correct
-    /// result the strides are not positive.
+    /// Assumes that no stride value appears twice.
     #[doc(hidden)]
     fn _fastest_varying_stride_order(&self) -> Self {
         let mut indices = self.clone();
         for (i, elt) in enumerate(indices.slice_mut()) {
             *elt = i;
         }
         let strides = self.slice();
-        indices.slice_mut().sort_by_key(|&i| strides[i]);
+        indices
+            .slice_mut()
+            .sort_by_key(|&i| (strides[i] as isize).abs());
         indices
     }
 
@@ -646,7 +646,7 @@ impl Dimension for Dim<[Ix; 2]> {
 
     #[inline]
     fn _fastest_varying_stride_order(&self) -> Self {
-        if get!(self, 0) as Ixs <= get!(self, 1) as Ixs {
+        if (get!(self, 0) as Ixs).abs() <= (get!(self, 1) as Ixs).abs() {
             Ix2(0, 1)
         } else {
             Ix2(1, 0)
@@ -806,7 +806,7 @@ impl Dimension for Dim<[Ix; 3]> {
         let mut order = Ix3(0, 1, 2);
         macro_rules! swap {
             ($stride:expr, $order:expr, $x:expr, $y:expr) => {
-                if $stride[$x] > $stride[$y] {
+                if ($stride[$x] as isize).abs() > ($stride[$y] as isize).abs() {
                     $stride.swap($x, $y);
                     $order.ixm().swap($x, $y);
                 }

src/dimension/mod.rs

@@ -46,15 +46,12 @@ pub fn stride_offset(n: Ix, stride: Ix) -> isize {
 /// There is overlap if, when iterating through the dimensions in order of
 /// increasing stride, the current stride is less than or equal to the maximum
 /// possible offset along the preceding axes. (Axes of length ≤1 are ignored.)
-///
-/// The current implementation assumes that strides of axes with length > 1 are
-/// nonnegative. Additionally, it does not check for overflow.
 pub fn dim_stride_overlap<D: Dimension>(dim: &D, strides: &D) -> bool {
     let order = strides._fastest_varying_stride_order();
     let mut sum_prev_offsets = 0;
     for &index in order.slice() {
         let d = dim[index];
-        let s = strides[index] as isize;
+        let s = (strides[index] as isize).abs();
         match d {
             0 => return false,
             1 => {}
@@ -210,11 +207,7 @@ where
 ///
 /// 2. The product of non-zero axis lengths must not exceed `isize::MAX`.
 ///
-/// 3. For axes with length > 1, the stride must be nonnegative. This is
-///    necessary to make sure the pointer cannot move backwards outside the
-///    slice. For axes with length ≤ 1, the stride can be anything.
-///
-/// 4. If the array will be empty (any axes are zero-length), the difference
+/// 3. If the array will be empty (any axes are zero-length), the difference
 ///    between the least address and greatest address accessible by moving
 ///    along all axes must be ≤ `data.len()`. (It's fine in this case to move
 ///    one byte past the end of the slice since the pointers will be offset but
@@ -225,13 +218,19 @@ where
 ///    `data.len()`. This and #3 ensure that all dereferenceable pointers point
 ///    to elements within the slice.
 ///
-/// 5. The strides must not allow any element to be referenced by two different
+/// 4. The strides must not allow any element to be referenced by two different
 ///    indices.
 ///
 /// Note that since slices cannot contain more than `isize::MAX` bytes,
 /// condition 4 is sufficient to guarantee that the absolute difference in
 /// units of `A` and in units of bytes between the least address and greatest
 /// address accessible by moving along all axes does not exceed `isize::MAX`.
+///
+/// Warning: This function is sufficient to check the invariants of ArrayBase only
+/// if the pointer to the first element of the array is chosen such that the element
+/// with the smallest memory address is at the start of data. (In other words, the
+/// pointer to the first element of the array must be computed using offset_from_ptr_to_memory
+/// so that negative strides are correctly handled.)
 pub(crate) fn can_index_slice<A, D: Dimension>(
     data: &[A],
     dim: &D,
@@ -248,7 +247,7 @@ fn can_index_slice_impl<D: Dimension>(
     dim: &D,
     strides: &D,
 ) -> Result<(), ShapeError> {
-    // Check condition 4.
+    // Check condition 3.
     let is_empty = dim.slice().iter().any(|&d| d == 0);
     if is_empty && max_offset > data_len {
         return Err(from_kind(ErrorKind::OutOfBounds));
@@ -257,15 +256,7 @@ fn can_index_slice_impl<D: Dimension>(
         return Err(from_kind(ErrorKind::OutOfBounds));
     }
 
-    // Check condition 3.
-    for (&d, &s) in izip!(dim.slice(), strides.slice()) {
-        let s = s as isize;
-        if d > 1 && s < 0 {
-            return Err(from_kind(ErrorKind::Unsupported));
-        }
-    }
-
-    // Check condition 5.
+    // Check condition 4.
     if !is_empty && dim_stride_overlap(dim, strides) {
         return Err(from_kind(ErrorKind::Unsupported));
     }
@@ -289,6 +280,19 @@ pub fn stride_offset_checked(dim: &[Ix], strides: &[Ix], index: &[Ix]) -> Option
     Some(offset)
 }
 
+/// Checks if strides are non-negative.
+pub fn strides_non_negative<D>(strides: &D) -> Result<(), ShapeError>
+where
+    D: Dimension,
+{
+    for &stride in strides.slice() {
+        if (stride as isize) < 0 {
+            return Err(from_kind(ErrorKind::Unsupported));
+        }
+    }
+    Ok(())
+}
+
 /// Implementation-specific extensions to `Dimension`
 pub trait DimensionExt {
     // note: many extensions go in the main trait if they need to be special-
@@ -394,6 +398,19 @@ fn to_abs_slice(axis_len: usize, slice: Slice) -> (usize, usize, isize) {
     (start, end, step)
 }
 
+/// This function computes the offset from the logically first element to the first element in
+/// memory of the array. The result is always <= 0.
+pub fn offset_from_ptr_to_memory<D: Dimension>(dim: &D, strides: &D) -> isize {
+    let offset = izip!(dim.slice(), strides.slice()).fold(0, |_offset, (d, s)| {
+        if (*s as isize) < 0 {
+            _offset + *s as isize * (*d as isize - 1)
+        } else {
+            _offset
+        }
+    });
+    offset
+}
+
 /// Modify dimension, stride and return data pointer offset
 ///
 /// **Panics** if stride is 0 or if any index is out of bounds.
@@ -693,13 +710,21 @@ mod test {
         let dim = (2, 3, 2).into_dimension();
         let strides = (5, 2, 1).into_dimension();
         assert!(super::dim_stride_overlap(&dim, &strides));
+        let strides = (-5isize as usize, 2, -1isize as usize).into_dimension();
+        assert!(super::dim_stride_overlap(&dim, &strides));
         let strides = (6, 2, 1).into_dimension();
         assert!(!super::dim_stride_overlap(&dim, &strides));
+        let strides = (6, -2isize as usize, 1).into_dimension();
+        assert!(!super::dim_stride_overlap(&dim, &strides));
         let strides = (6, 0, 1).into_dimension();
         assert!(super::dim_stride_overlap(&dim, &strides));
+        let strides = (-6isize as usize, 0, 1).into_dimension();
+        assert!(super::dim_stride_overlap(&dim, &strides));
         let dim = (2, 2).into_dimension();
         let strides = (3, 2).into_dimension();
         assert!(!super::dim_stride_overlap(&dim, &strides));
+        let strides = (3, -2isize as usize).into_dimension();
+        assert!(!super::dim_stride_overlap(&dim, &strides));
     }
 
     #[test]
@@ -736,7 +761,7 @@ mod test {
         can_index_slice::<i32, _>(&[1], &Ix1(2), &Ix1(1)).unwrap_err();
         can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(0)).unwrap_err();
         can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(1)).unwrap();
-        can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(-1isize as usize)).unwrap_err();
+        can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(-1isize as usize)).unwrap();
     }
 
     #[test]

src/impl_constructors.rs

@@ -19,6 +19,7 @@ use alloc::vec;
 use alloc::vec::Vec;
 
 use crate::dimension;
+use crate::dimension::offset_from_ptr_to_memory;
 use crate::error::{self, ShapeError};
 use crate::extension::nonnull::nonnull_from_vec_data;
 use crate::imp_prelude::*;
@@ -30,6 +31,7 @@ use crate::iterators::to_vec_mapped;
 use crate::StrideShape;
 #[cfg(feature = "std")]
 use crate::{geomspace, linspace, logspace};
+use rawpointer::PointerExt;
 
 
 /// # Constructor Methods for Owned Arrays
@@ -442,7 +444,8 @@ where
     ///
     /// 2. The product of non-zero axis lengths must not exceed `isize::MAX`.
     ///
-    /// 3. For axes with length > 1, the stride must be nonnegative.
+    /// 3. For axes with length > 1, the pointer cannot move outside the
+    ///    slice.
     ///
     /// 4. If the array will be empty (any axes are zero-length), the
     ///    difference between the least address and greatest address accessible
@@ -468,7 +471,7 @@ where
         // debug check for issues that indicates wrong use of this constructor
         debug_assert!(dimension::can_index_slice(&v, &dim, &strides).is_ok());
         ArrayBase {
-            ptr: nonnull_from_vec_data(&mut v),
+            ptr: nonnull_from_vec_data(&mut v).offset(-offset_from_ptr_to_memory(&dim, &strides)),
             data: DataOwned::new(v),
             strides,
             dim,
@@ -494,7 +497,7 @@ where
     ///
     /// This constructor is limited to elements where `A: Copy` (no destructors)
     /// to avoid users shooting themselves too hard in the foot.
-    /// 
+    ///
     /// (Also note that the constructors `from_shape_vec` and
     /// `from_shape_vec_unchecked` allow the user yet more control, in the sense
     /// that Arrays can be created from arbitrary vectors.)

src/impl_methods.rs

@@ -18,7 +18,8 @@ use crate::arraytraits;
 use crate::dimension;
 use crate::dimension::IntoDimension;
 use crate::dimension::{
-    abs_index, axes_of, do_slice, merge_axes, size_of_shape_checked, stride_offset, Axes,
+    abs_index, axes_of, do_slice, merge_axes, offset_from_ptr_to_memory, size_of_shape_checked,
+    stride_offset, Axes,
 };
 use crate::error::{self, ErrorKind, ShapeError};
 use crate::math_cell::MathCell;
@@ -169,12 +170,12 @@ where
 
     /// Return an uniquely owned copy of the array.
     ///
-    /// If the input array is contiguous and its strides are positive, then the
-    /// output array will have the same memory layout. Otherwise, the layout of
-    /// the output array is unspecified. If you need a particular layout, you
-    /// can allocate a new array with the desired memory layout and
-    /// [`.assign()`](#method.assign) the data. Alternatively, you can collect
-    /// an iterator, like this for a result in standard layout:
+    /// If the input array is contiguous, then the output array will have the same
+    /// memory layout. Otherwise, the layout of the output array is unspecified.
+    /// If you need a particular layout, you can allocate a new array with the
+    /// desired memory layout and [`.assign()`](#method.assign) the data.
+    /// Alternatively, you can collectan iterator, like this for a result in
+    /// standard layout:
     ///
     /// ```
     /// # use ndarray::prelude::*;
@@ -1296,9 +1297,6 @@ where
     }
 
     /// Return true if the array is known to be contiguous.
-    ///
-    /// Will detect c- and f-contig arrays correctly, but otherwise
-    /// There are some false negatives.
     pub(crate) fn is_contiguous(&self) -> bool {
         D::is_contiguous(&self.dim, &self.strides)
     }
@@ -1420,14 +1418,18 @@ where
     ///
     /// If this function returns `Some(_)`, then the elements in the slice
     /// have whatever order the elements have in memory.
-    ///
-    /// Implementation notes: Does not yet support negatively strided arrays.
     pub fn as_slice_memory_order(&self) -> Option<&[A]>
     where
         S: Data,
     {
         if self.is_contiguous() {
-            unsafe { Some(slice::from_raw_parts(self.ptr.as_ptr(), self.len())) }
+            let offset = offset_from_ptr_to_memory(&self.dim, &self.strides);
+            unsafe {
+                Some(slice::from_raw_parts(
+                    self.ptr.offset(offset).as_ptr(),
+                    self.len(),
+                ))
+            }
         } else {
             None
         }
@@ -1441,7 +1443,13 @@ where
     {
         if self.is_contiguous() {
             self.ensure_unique();
-            unsafe { Some(slice::from_raw_parts_mut(self.ptr.as_ptr(), self.len())) }
+            let offset = offset_from_ptr_to_memory(&self.dim, &self.strides);
+            unsafe {
+                Some(slice::from_raw_parts_mut(
+                    self.ptr.offset(offset).as_ptr(),
+                    self.len(),
+                ))
+            }
         } else {
             None
         }

src/impl_raw_views.rs

@@ -62,6 +62,11 @@ where
     ///     [`.offset()`] regardless of the starting point due to past offsets.
     ///
     /// * The product of non-zero axis lengths must not exceed `isize::MAX`.
+    /// 
+    /// * Strides must be non-negative.
+    ///
+    /// This function can use debug assertions to check some of these requirements,
+    /// but it's not a complete check.
     ///
     /// [`.offset()`]: https://doc.rust-lang.org/stable/std/primitive.pointer.html#method.offset
     pub unsafe fn from_shape_ptr<Sh>(shape: Sh, ptr: *const A) -> Self
@@ -73,6 +78,7 @@ where
         if cfg!(debug_assertions) {
             assert!(!ptr.is_null(), "The pointer must be non-null.");
             if let Strides::Custom(strides) = &shape.strides {
+                dimension::strides_non_negative(strides).unwrap();
                 dimension::max_abs_offset_check_overflow::<A, _>(&dim, &strides).unwrap();
             } else {
                 dimension::size_of_shape_checked(&dim).unwrap();
@@ -202,6 +208,11 @@ where
     ///     [`.offset()`] regardless of the starting point due to past offsets.
     ///
     /// * The product of non-zero axis lengths must not exceed `isize::MAX`.
+    /// 
+    /// * Strides must be non-negative.
+    ///
+    /// This function can use debug assertions to check some of these requirements,
+    /// but it's not a complete check.
     ///
     /// [`.offset()`]: https://doc.rust-lang.org/stable/std/primitive.pointer.html#method.offset
     pub unsafe fn from_shape_ptr<Sh>(shape: Sh, ptr: *mut A) -> Self
@@ -213,6 +224,7 @@ where
         if cfg!(debug_assertions) {
             assert!(!ptr.is_null(), "The pointer must be non-null.");
             if let Strides::Custom(strides) = &shape.strides {
+                dimension::strides_non_negative(strides).unwrap();
                 dimension::max_abs_offset_check_overflow::<A, _>(&dim, &strides).unwrap();
             } else {
                 dimension::size_of_shape_checked(&dim).unwrap();