Repository: rust-ndarray/ndarray
Pull request #1022 (2021-06-01 23:32:11)

Add try_into_owned_nocopy method
I've been working on improvements to the layout handling in `ndarray-linalg` to make the code easier to understand, fix bugs, and avoid unnecessary cloning and by-clone transposes. One of operations I need to implement is, "Convert an `ArrayBase` into a LAPACK-compatible `Array` (i.e. `strides = [1, s]`, with `s >= nrows`), while avoiding unnecessary cloning." I don't think it's possible to implement this today without adding more methods to `ArrayBase`. To address this, I propose adding a `try_into_owned_nocopy` method which is similar to `into_owned` but returns `Err(self)` if the data would need to be cloned, instead of actually performing the cloning. This way, if the data needs to be cloned, I can clone it into an array of the correct layout, like this:

```rust
match array.try_into_owned_nocopy() {
    Ok(o) if is_lapack_compatible(&o) => o,
    Ok(o) => clone_into_new_lapack_array(&o),
    Err(a) => clone_into_new_lapack_array(&a),
}
```

There are other ways this problem could be solved, but `try_into_owned_nocopy` is the best I came up with.

A couple of notes about the PR:

1. ~If we implement `DataOwned` for `CowRepr` in another PR, then I'd change the bound on `try_into_owned_nocopy` to be `S: DataOwned`, since there wouldn't be any reason to call it on any non-`DataOwned` types.~ Edit: I was wrong; the method is actually useful for `S: Data` after all, and there isn't a strong reason to restrict it to `S: DataOwned` anyway.

2. I'm not completely satisfied with the name `try_into_owned_nocopy`. Other suggestions would be welcome.


Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 42e5e658..a44568ae 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -75,5 +75,5 @@ jobs:
           toolchain: ${{ matrix.rust }}
           override: true
           components: clippy
-      - run: cargo clippy
+      - run: cargo clippy --features docs
 
diff --git a/Cargo.toml b/Cargo.toml
index 31bab8d3..75875d61 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,7 +1,7 @@
 [package]
 
 name = "ndarray"
-version = "0.15.2"
+version = "0.15.3"
 edition = "2018"
 authors = [
   "Ulrik Sverdrup \"bluss\"",
diff --git a/README.rst b/README.rst
index 50302488..a6ed9e90 100644
--- a/README.rst
+++ b/README.rst
@@ -15,7 +15,7 @@ __ https://docs.rs/ndarray/
    :alt: CI build status
 .. _build_status: https://github.com/rust-ndarray/ndarray/actions
 
-.. |crates| image:: http://meritbadge.herokuapp.com/ndarray
+.. |crates| image:: https://img.shields.io/crates/v/ndarray.svg
    :alt: ndarray at crates.io
 .. _crates: https://crates.io/crates/ndarray
 
diff --git a/RELEASES.md b/RELEASES.md
index da895d94..cfb5e7d3 100644
--- a/RELEASES.md
+++ b/RELEASES.md
@@ -1,3 +1,54 @@
+Version 0.15.3 (2021-06-05)
+===========================
+
+New features
+------------
+
+- New methods `.last/_mut()` for arrays and array views by [@jturner314]
+
+  https://github.com/rust-ndarray/ndarray/pull/1013
+
+Bug fixes
+---------
+
+- Fix `as_slice_memory_order_mut()` so that it never changes strides (the
+  memory layout) of the array when called.
+
+  This was a bug that impacted `ArcArray` (and for example not `Array` or `ArrayView/Mut`),
+  and multiple methods on `ArcArray` that use `as_slice_memory_order_mut` (for example `map_mut`).
+  Fix by [@jturner314].
+
+  https://github.com/rust-ndarray/ndarray/pull/1019
+
+API changes
+-----------
+
+- Array1 now implements `From<Box<[T]>>` by [@jturner314]
+
+  https://github.com/rust-ndarray/ndarray/pull/1016
+
+- ArcArray now implements `From<Array<...>>` by [@jturner314]
+
+  https://github.com/rust-ndarray/ndarray/pull/1021
+
+- CowArray now implements RawDataSubst by [@jturner314]
+
+  https://github.com/rust-ndarray/ndarray/pull/1020
+
+Other changes
+-------------
+
+- Mention unsharing in `.as_mut_ptr` docs by [@jturner314]
+
+  https://github.com/rust-ndarray/ndarray/pull/1017
+
+- Clarify and fix minor errors in push/append method docs by [@bluss] f21c668a
+
+- Fix several warnings in doc example code by [@bluss]
+
+  https://github.com/rust-ndarray/ndarray/pull/1009
+
+
 Version 0.15.2 (2021-05-17 🇳🇴)
 ================================
 
diff --git a/ndarray-rand/src/lib.rs b/ndarray-rand/src/lib.rs
index e293eb56..d9edabfd 100644
--- a/ndarray-rand/src/lib.rs
+++ b/ndarray-rand/src/lib.rs
@@ -8,15 +8,15 @@
 
 //! Constructors for randomized arrays: `rand` integration for `ndarray`.
 //!
-//! See [**`RandomExt`**](trait.RandomExt.html) for usage examples.
+//! See **[`RandomExt`]** for usage examples.
 //!
 //! ## Note
 //!
 //! `ndarray-rand` depends on [`rand` 0.8][rand].
 //!
 //! [`rand`][rand] and [`rand_distr`][rand_distr]
-//! are re-exported as sub-modules, [`ndarray_rand::rand`](rand/index.html)
-//! and [`ndarray_rand::rand_distr`](rand_distr/index.html) respectively.
+//! are re-exported as sub-modules, [`ndarray_rand::rand`](rand)
+//! and [`ndarray_rand::rand_distr`](rand_distr) respectively.
 //! You can use these submodules for guaranteed version compatibility or
 //! convenience.
 //!
@@ -60,7 +60,7 @@ pub mod rand_distr {
 /// Note that `SmallRng` is cheap to initialize and fast, but it may generate
 /// low-quality random numbers, and reproducibility is not guaranteed. See its
 /// documentation for information. You can select a different RNG with
-/// [`.random_using()`](#tymethod.random_using).
+/// [`.random_using()`](Self::random_using).
 pub trait RandomExt<S, A, D>
 where
     S: RawData<Elem = A>,
@@ -293,8 +293,8 @@ where
 /// if lanes from the original array should only be sampled once (*without replacement*) or
 /// multiple times (*with replacement*).
 ///
-/// [`sample_axis`]: trait.RandomExt.html#tymethod.sample_axis
-/// [`sample_axis_using`]: trait.RandomExt.html#tymethod.sample_axis_using
+/// [`sample_axis`]: RandomExt::sample_axis
+/// [`sample_axis_using`]: RandomExt::sample_axis_using
 #[derive(Debug, Clone)]
 pub enum SamplingStrategy {
     WithReplacement,
diff --git a/src/array_serde.rs b/src/array_serde.rs
index 41f7d60c..f1c9e121 100644
--- a/src/array_serde.rs
+++ b/src/array_serde.rs
@@ -276,17 +276,17 @@ where
 
         let _v = match v {
             Some(v) => v,
-            None => Err(de::Error::missing_field("v"))?,
+            None => return Err(de::Error::missing_field("v")),
         };
 
         let data = match data {
             Some(data) => data,
-            None => Err(de::Error::missing_field("data"))?,
+            None => return Err(de::Error::missing_field("data")),
         };
 
         let dim = match dim {
             Some(dim) => dim,
-            None => Err(de::Error::missing_field("dim"))?,
+            None => return Err(de::Error::missing_field("dim")),
         };
 
         if let Ok(array) = ArrayBase::from_shape_vec(dim, data) {
diff --git a/src/arraytraits.rs b/src/arraytraits.rs
index 421fc471..a7f22c1f 100644
--- a/src/arraytraits.rs
+++ b/src/arraytraits.rs
@@ -11,6 +11,7 @@ use std::iter::FromIterator;
 use std::iter::IntoIterator;
 use std::mem;
 use std::ops::{Index, IndexMut};
+use alloc::boxed::Box;
 use alloc::vec::Vec;
 
 use crate::imp_prelude::*;
@@ -148,6 +149,18 @@ where
 {
 }
 
+impl<A, S> From<Box<[A]>> for ArrayBase<S, Ix1>
+where
+    S: DataOwned<Elem = A>,
+{
+    /// Create a one-dimensional array from a boxed slice (no copying needed).
+    ///
+    /// **Panics** if the length is greater than `isize::MAX`.
+    fn from(b: Box<[A]>) -> Self {
+        Self::from_vec(b.into_vec())
+    }
+}
+
 impl<A, S> From<Vec<A>> for ArrayBase<S, Ix1>
 where
     S: DataOwned<Elem = A>,
@@ -354,6 +367,15 @@ where
     }
 }
 
+impl<A, D> From<Array<A, D>> for ArcArray<A, D>
+where
+    D: Dimension,
+{
+    fn from(arr: Array<A, D>) -> ArcArray<A, D> {
+        arr.into_shared()
+    }
+}
+
 /// Argument conversion into an array view
 ///
 /// The trait is parameterized over `A`, the element type, and `D`, the
diff --git a/src/data_repr.rs b/src/data_repr.rs
index 7047f201..6630f9dd 100644
--- a/src/data_repr.rs
+++ b/src/data_repr.rs
@@ -11,7 +11,7 @@ use rawpointer::PointerExt;
 /// Array's representation.
 ///
 /// *Don’t use this type directly—use the type alias
-/// [`Array`](type.Array.html) for the array type!*
+/// [`Array`](crate::Array) for the array type!*
 // Like a Vec, but with non-unique ownership semantics
 //
 // repr(C) to make it transmutable OwnedRepr<A> -> OwnedRepr<B> if
diff --git a/src/data_traits.rs b/src/data_traits.rs
index bd540b1b..b117792c 100644
--- a/src/data_traits.rs
+++ b/src/data_traits.rs
@@ -16,7 +16,9 @@ use std::ptr::NonNull;
 use alloc::sync::Arc;
 use alloc::vec::Vec;
 
-use crate::{ArrayBase, CowRepr, Dimension, OwnedArcRepr, OwnedRepr, RawViewRepr, ViewRepr};
+use crate::{
+    ArcArray, Array, ArrayBase, CowRepr, Dimension, OwnedArcRepr, OwnedRepr, RawViewRepr, ViewRepr,
+};
 
 /// Array representation trait.
 ///
@@ -27,6 +29,7 @@ use crate::{ArrayBase, CowRepr, Dimension, OwnedArcRepr, OwnedRepr, RawViewRepr,
 /// ***Note:*** `RawData` is not an extension interface at this point.
 /// Traits in Rust can serve many different roles. This trait is public because
 /// it is used as a bound on public methods.
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait RawData: Sized {
     /// The array element type.
     type Elem;
@@ -47,11 +50,15 @@ pub unsafe trait RawData: Sized {
 /// For an array with writable elements.
 ///
 /// ***Internal trait, see `RawData`.***
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait RawDataMut: RawData {
     /// If possible, ensures that the array has unique access to its data.
     ///
-    /// If `Self` provides safe mutable access to array elements, then it
-    /// **must** panic or ensure that the data is unique.
+    /// The implementer must ensure that if the input is contiguous, then the
+    /// output has the same strides as input.
+    ///
+    /// Additionally, if `Self` provides safe mutable access to array elements,
+    /// then this method **must** panic or ensure that the data is unique.
     #[doc(hidden)]
     fn try_ensure_unique<D>(_: &mut ArrayBase<Self, D>)
     where
@@ -71,6 +78,7 @@ pub unsafe trait RawDataMut: RawData {
 /// An array representation that can be cloned.
 ///
 /// ***Internal trait, see `RawData`.***
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait RawDataClone: RawData {
     #[doc(hidden)]
     /// Unsafe because, `ptr` must point inside the current storage.
@@ -93,20 +101,30 @@ pub unsafe trait RawDataClone: RawData {
 /// For an array with elements that can be accessed with safe code.
 ///
 /// ***Internal trait, see `RawData`.***
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait Data: RawData {
     /// Converts the array to a uniquely owned array, cloning elements if necessary.
     #[doc(hidden)]
     #[allow(clippy::wrong_self_convention)]
-    fn into_owned<D>(self_: ArrayBase<Self, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(self_: ArrayBase<Self, D>) -> Array<Self::Elem, D>
     where
         Self::Elem: Clone,
         D: Dimension;
 
+    /// Converts the array into `Array<A, D>` if this is possible without
+    /// cloning the array elements. Otherwise, returns `self_` unchanged.
+    #[doc(hidden)]
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension;
+
     /// Return a shared ownership (copy on write) array based on the existing one,
     /// cloning elements if necessary.
     #[doc(hidden)]
     #[allow(clippy::wrong_self_convention)]
-    fn to_shared<D>(self_: &ArrayBase<Self, D>) -> ArrayBase<OwnedArcRepr<Self::Elem>, D>
+    fn to_shared<D>(self_: &ArrayBase<Self, D>) -> ArcArray<Self::Elem, D>
     where
         Self::Elem: Clone,
         D: Dimension,
@@ -128,6 +146,7 @@ pub unsafe trait Data: RawData {
 // `RawDataMut::try_ensure_unique` implementation always panics or ensures that
 // the data is unique. You are also guaranteeing that `try_is_unique` always
 // returns `Some(_)`.
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait DataMut: Data + RawDataMut {
     /// Ensures that the array has unique access to its data.
     #[doc(hidden)]
@@ -230,14 +249,9 @@ where
             return;
         }
         if self_.dim.size() <= self_.data.0.len() / 2 {
-            // Create a new vec if the current view is less than half of
-            // backing data.
-            unsafe {
-                *self_ = ArrayBase::from_shape_vec_unchecked(
-                    self_.dim.clone(),
-                    self_.iter().cloned().collect(),
-                );
-            }
+            // Clone only the visible elements if the current view is less than
+            // half of backing data.
+            *self_ = self_.to_owned().into_shared();
             return;
         }
         let rcvec = &mut self_.data.0;
@@ -259,7 +273,7 @@ where
 }
 
 unsafe impl<A> Data for OwnedArcRepr<A> {
-    fn into_owned<D>(mut self_: ArrayBase<Self, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(mut self_: ArrayBase<Self, D>) -> Array<Self::Elem, D>
     where
         A: Clone,
         D: Dimension,
@@ -273,8 +287,29 @@ unsafe impl<A> Data for OwnedArcRepr<A> {
         }
     }
 
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension,
+    {
+        match Arc::try_unwrap(self_.data.0) {
+            Ok(owned_data) => unsafe {
+                // Safe because the data is equivalent.
+                Ok(ArrayBase::from_data_ptr(owned_data, self_.ptr)
+                    .with_strides_dim(self_.strides, self_.dim))
+            },
+            Err(arc_data) => unsafe {
+                // Safe because the data is equivalent; we're just
+                // reconstructing `self_`.
+                Err(ArrayBase::from_data_ptr(OwnedArcRepr(arc_data), self_.ptr)
+                    .with_strides_dim(self_.strides, self_.dim))
+            },
+        }
+    }
+
     #[allow(clippy::wrong_self_convention)]
-    fn to_shared<D>(self_: &ArrayBase<Self, D>) -> ArrayBase<OwnedArcRepr<Self::Elem>, D>
+    fn to_shared<D>(self_: &ArrayBase<Self, D>) -> ArcArray<Self::Elem, D>
     where
         Self::Elem: Clone,
         D: Dimension,
@@ -327,13 +362,23 @@ unsafe impl<A> RawDataMut for OwnedRepr<A> {
 
 unsafe impl<A> Data for OwnedRepr<A> {
     #[inline]
-    fn into_owned<D>(self_: ArrayBase<Self, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(self_: ArrayBase<Self, D>) -> Array<Self::Elem, D>
     where
         A: Clone,
         D: Dimension,
     {
         self_
     }
+
+    #[inline]
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension,
+    {
+        Ok(self_)
+    }
 }
 
 unsafe impl<A> DataMut for OwnedRepr<A> {}
@@ -363,7 +408,7 @@ where
         } else {
             0
         };
-        self.clone_from(&other);
+        self.clone_from(other);
         self.as_nonnull_mut().offset(our_off)
     }
 }
@@ -383,13 +428,22 @@ unsafe impl<'a, A> RawData for ViewRepr<&'a A> {
 }
 
 unsafe impl<'a, A> Data for ViewRepr<&'a A> {
-    fn into_owned<D>(self_: ArrayBase<Self, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(self_: ArrayBase<Self, D>) -> Array<Self::Elem, D>
     where
         Self::Elem: Clone,
         D: Dimension,
     {
         self_.to_owned()
     }
+
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension,
+    {
+        Err(self_)
+    }
 }
 
 unsafe impl<'a, A> RawDataClone for ViewRepr<&'a A> {
@@ -428,13 +482,22 @@ unsafe impl<'a, A> RawDataMut for ViewRepr<&'a mut A> {
 }
 
 unsafe impl<'a, A> Data for ViewRepr<&'a mut A> {
-    fn into_owned<D>(self_: ArrayBase<Self, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(self_: ArrayBase<Self, D>) -> Array<Self::Elem, D>
     where
         Self::Elem: Clone,
         D: Dimension,
     {
         self_.to_owned()
     }
+
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension,
+    {
+        Err(self_)
+    }
 }
 
 unsafe impl<'a, A> DataMut for ViewRepr<&'a mut A> {}
@@ -451,6 +514,7 @@ unsafe impl<'a, A> DataMut for ViewRepr<&'a mut A> {}
 // The array storage must be initially mutable - copy on write arrays may require copying for
 // unsharing storage before mutating it. The initially allocated storage must be mutable so
 // that it can be mutated directly - through .raw_view_mut_unchecked() - for initialization.
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait DataOwned: Data {
     /// Corresponding owned data with MaybeUninit elements
     type MaybeUninit: DataOwned<Elem = MaybeUninit<Self::Elem>>
@@ -469,6 +533,7 @@ pub unsafe trait DataOwned: Data {
 /// A representation that is a lightweight view.
 ///
 /// ***Internal trait, see `Data`.***
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait DataShared: Clone + Data + RawDataClone {}
 
 unsafe impl<A> DataShared for OwnedArcRepr<A> {}
@@ -590,7 +655,7 @@ where
 
 unsafe impl<'a, A> Data for CowRepr<'a, A> {
     #[inline]
-    fn into_owned<D>(self_: ArrayBase<CowRepr<'a, A>, D>) -> ArrayBase<OwnedRepr<Self::Elem>, D>
+    fn into_owned<D>(self_: ArrayBase<CowRepr<'a, A>, D>) -> Array<Self::Elem, D>
     where
         A: Clone,
         D: Dimension,
@@ -604,6 +669,22 @@ unsafe impl<'a, A> Data for CowRepr<'a, A> {
             },
         }
     }
+
+    fn try_into_owned_nocopy<D>(
+        self_: ArrayBase<Self, D>,
+    ) -> Result<Array<Self::Elem, D>, ArrayBase<Self, D>>
+    where
+        D: Dimension,
+    {
+        match self_.data {
+            CowRepr::View(_) => Err(self_),
+            CowRepr::Owned(data) => unsafe {
+                // safe because the data is equivalent so ptr, dims remain valid
+                Ok(ArrayBase::from_data_ptr(data, self_.ptr)
+                    .with_strides_dim(self_.strides, self_.dim))
+            },
+        }
+    }
 }
 
 unsafe impl<'a, A> DataMut for CowRepr<'a, A> where A: Clone {}
@@ -675,3 +756,13 @@ impl<'a, A: 'a, B: 'a> RawDataSubst<B> for ViewRepr<&'a mut A> {
     }
 }
 
+impl<'a, A: 'a, B: 'a> RawDataSubst<B> for CowRepr<'a, A> {
+    type Output = CowRepr<'a, B>;
+
+    unsafe fn data_subst(self) -> Self::Output {
+        match self {
+            CowRepr::View(view) => CowRepr::View(view.data_subst()),
+            CowRepr::Owned(owned) => CowRepr::Owned(owned.data_subst()),
+        }
+    }
+}
diff --git a/src/dimension/dim.rs b/src/dimension/dim.rs
index 97de6d84..19075f94 100644
--- a/src/dimension/dim.rs
+++ b/src/dimension/dim.rs
@@ -18,7 +18,7 @@ use crate::Ix;
 /// `Dim` describes the number of axes and the length of each axis
 /// in an array. It is also used as an index type.
 ///
-/// See also the [`Dimension` trait](trait.Dimension.html) for its methods and
+/// See also the [`Dimension`] trait for its methods and
 /// operations.
 ///
 /// # Examples
diff --git a/src/dimension/dimension_trait.rs b/src/dimension/dimension_trait.rs
index 1eace34a..a4d73b1b 100644
--- a/src/dimension/dimension_trait.rs
+++ b/src/dimension/dimension_trait.rs
@@ -379,6 +379,7 @@ pub trait Dimension:
 
 macro_rules! impl_insert_axis_array(
     ($n:expr) => (
+        #[inline]
         fn insert_axis(&self, axis: Axis) -> Self::Larger {
             debug_assert!(axis.index() <= $n);
             let mut out = [1; $n + 1];
@@ -422,7 +423,6 @@ impl Dimension for Dim<[Ix; 0]> {
     fn next_for(&self, _index: Self) -> Option<Self> {
         None
     }
-    #[inline]
     impl_insert_axis_array!(0);
     #[inline]
     fn try_remove_axis(&self, _ignore: Axis) -> Self::Smaller {
@@ -530,7 +530,6 @@ impl Dimension for Dim<[Ix; 1]> {
             None
         }
     }
-    #[inline]
     impl_insert_axis_array!(1);
     #[inline]
     fn try_remove_axis(&self, axis: Axis) -> Self::Smaller {
@@ -694,7 +693,6 @@ impl Dimension for Dim<[Ix; 2]> {
             None
         }
     }
-    #[inline]
     impl_insert_axis_array!(2);
     #[inline]
     fn try_remove_axis(&self, axis: Axis) -> Self::Smaller {
@@ -814,7 +812,6 @@ impl Dimension for Dim<[Ix; 3]> {
         }
         order
     }
-    #[inline]
     impl_insert_axis_array!(3);
     #[inline]
     fn try_remove_axis(&self, axis: Axis) -> Self::Smaller {
@@ -845,7 +842,6 @@ macro_rules! large_dim {
                 assert_eq!(ndim, $n);
                 Self::default()
             }
-            #[inline]
             $($insert_axis)*
             #[inline]
             fn try_remove_axis(&self, axis: Axis) -> Self::Smaller {
diff --git a/src/dimension/ndindex.rs b/src/dimension/ndindex.rs
index d9bac1d9..810a5b09 100644
--- a/src/dimension/ndindex.rs
+++ b/src/dimension/ndindex.rs
@@ -18,6 +18,7 @@ use crate::{
 /// a[[1, 1]] += 1;
 /// assert_eq!(a[(1, 1)], 4);
 /// ```
+#[allow(clippy::missing_safety_doc)] // TODO: Add doc
 pub unsafe trait NdIndex<E>: Debug {
     #[doc(hidden)]
     fn index_checked(&self, dim: &E, strides: &E) -> Option<isize>;
diff --git a/src/doc/ndarray_for_numpy_users/coord_transform.rs b/src/doc/ndarray_for_numpy_users/coord_transform.rs
index 09295be9..1529e874 100644
--- a/src/doc/ndarray_for_numpy_users/coord_transform.rs
+++ b/src/doc/ndarray_for_numpy_users/coord_transform.rs
@@ -51,40 +51,38 @@
 //! ```
 //! use ndarray::prelude::*;
 //!
-//! fn main() {
-//!     let nelems = 4;
-//!     let bunge = Array::ones((3, nelems));
-//!
-//!     let s1 = bunge.slice(s![0, ..]).mapv(f64::sin);
-//!     let c1 = bunge.slice(s![0, ..]).mapv(f64::cos);
-//!     let s2 = bunge.slice(s![1, ..]).mapv(f64::sin);
-//!     let c2 = bunge.slice(s![1, ..]).mapv(f64::cos);
-//!     let s3 = bunge.slice(s![2, ..]).mapv(f64::sin);
-//!     let c3 = bunge.slice(s![2, ..]).mapv(f64::cos);
-//!
-//!     let mut rmat = Array::zeros((3, 3, nelems).f());
-//!     for i in 0..nelems {
-//!         rmat[[0, 0, i]] = c1[i] * c3[i] - s1[i] * s3[i] * c2[i];
-//!         rmat[[0, 1, i]] = -c1[i] * s3[i] - s1[i] * c2[i] * c3[i];
-//!         rmat[[0, 2, i]] = s1[i] * s2[i];
-//!
-//!         rmat[[1, 0, i]] = s1[i] * c3[i] + c1[i] * c2[i] * s3[i];
-//!         rmat[[1, 1, i]] = -s1[i] * s3[i] + c1[i] * c2[i] * c3[i];
-//!         rmat[[1, 2, i]] = -c1[i] * s2[i];
-//!
-//!         rmat[[2, 0, i]] = s2[i] * s3[i];
-//!         rmat[[2, 1, i]] = s2[i] * c3[i];
-//!         rmat[[2, 2, i]] = c2[i];
-//!     }
-//!
-//!     let eye2d = Array::eye(3);
-//!
-//!     let mut rotated = Array::zeros((3, 3, nelems).f());
-//!     for i in 0..nelems {
-//!         rotated
-//!             .slice_mut(s![.., .., i])
-//!             .assign(&rmat.slice(s![.., .., i]).dot(&eye2d));
-//!     }
+//! let nelems = 4;
+//! let bunge = Array::ones((3, nelems));
+//!
+//! let s1 = bunge.slice(s![0, ..]).mapv(f64::sin);
+//! let c1 = bunge.slice(s![0, ..]).mapv(f64::cos);
+//! let s2 = bunge.slice(s![1, ..]).mapv(f64::sin);
+//! let c2 = bunge.slice(s![1, ..]).mapv(f64::cos);
+//! let s3 = bunge.slice(s![2, ..]).mapv(f64::sin);
+//! let c3 = bunge.slice(s![2, ..]).mapv(f64::cos);
+//!
+//! let mut rmat = Array::zeros((3, 3, nelems).f());
+//! for i in 0..nelems {
+//!     rmat[[0, 0, i]] = c1[i] * c3[i] - s1[i] * s3[i] * c2[i];
+//!     rmat[[0, 1, i]] = -c1[i] * s3[i] - s1[i] * c2[i] * c3[i];
+//!     rmat[[0, 2, i]] = s1[i] * s2[i];
+//!
+//!     rmat[[1, 0, i]] = s1[i] * c3[i] + c1[i] * c2[i] * s3[i];
+//!     rmat[[1, 1, i]] = -s1[i] * s3[i] + c1[i] * c2[i] * c3[i];
+//!     rmat[[1, 2, i]] = -c1[i] * s2[i];
+//!
+//!     rmat[[2, 0, i]] = s2[i] * s3[i];
+//!     rmat[[2, 1, i]] = s2[i] * c3[i];
+//!     rmat[[2, 2, i]] = c2[i];
+//! }
+//!
+//! let eye2d = Array::eye(3);
+//!
+//! let mut rotated = Array::zeros((3, 3, nelems).f());
+//! for i in 0..nelems {
+//!     rotated
+//!         .slice_mut(s![.., .., i])
+//!         .assign(&rmat.slice(s![.., .., i]).dot(&eye2d));
 //! }
 //! ```
 //!
@@ -96,37 +94,35 @@
 //! ```
 //! use ndarray::prelude::*;
 //!
-//! fn main() {
-//!     let nelems = 4;
-//!     let bunge = Array2::<f64>::ones((3, nelems));
-//!
-//!     let mut rmat = Array::zeros((3, 3, nelems).f());
-//!     azip!((mut rmat in rmat.axis_iter_mut(Axis(2)), bunge in bunge.axis_iter(Axis(1))) {
-//!         let s1 = bunge[0].sin();
-//!         let c1 = bunge[0].cos();
-//!         let s2 = bunge[1].sin();
-//!         let c2 = bunge[1].cos();
-//!         let s3 = bunge[2].sin();
-//!         let c3 = bunge[2].cos();
-//!
-//!         rmat[[0, 0]] = c1 * c3 - s1 * s3 * c2;
-//!         rmat[[0, 1]] = -c1 * s3 - s1 * c2 * c3;
-//!         rmat[[0, 2]] = s1 * s2;
-//!
-//!         rmat[[1, 0]] = s1 * c3 + c1 * c2 * s3;
-//!         rmat[[1, 1]] = -s1 * s3 + c1 * c2 * c3;
-//!         rmat[[1, 2]] = -c1 * s2;
-//!
-//!         rmat[[2, 0]] = s2 * s3;
-//!         rmat[[2, 1]] = s2 * c3;
-//!         rmat[[2, 2]] = c2;
-//!     });
-//!
-//!     let eye2d = Array2::<f64>::eye(3);
-//!
-//!     let mut rotated = Array3::<f64>::zeros((3, 3, nelems).f());
-//!     azip!((mut rotated in rotated.axis_iter_mut(Axis(2)), rmat in rmat.axis_iter(Axis(2))) {
-//!         rotated.assign(&rmat.dot(&eye2d));
-//!     });
-//! }
+//! let nelems = 4;
+//! let bunge = Array2::<f64>::ones((3, nelems));
+//!
+//! let mut rmat = Array::zeros((3, 3, nelems).f());
+//! azip!((mut rmat in rmat.axis_iter_mut(Axis(2)), bunge in bunge.axis_iter(Axis(1))) {
+//!     let s1 = bunge[0].sin();
+//!     let c1 = bunge[0].cos();
+//!     let s2 = bunge[1].sin();
+//!     let c2 = bunge[1].cos();
+//!     let s3 = bunge[2].sin();
+//!     let c3 = bunge[2].cos();
+//!
+//!     rmat[[0, 0]] = c1 * c3 - s1 * s3 * c2;
+//!     rmat[[0, 1]] = -c1 * s3 - s1 * c2 * c3;
+//!     rmat[[0, 2]] = s1 * s2;
+//!
+//!     rmat[[1, 0]] = s1 * c3 + c1 * c2 * s3;
+//!     rmat[[1, 1]] = -s1 * s3 + c1 * c2 * c3;
+//!     rmat[[1, 2]] = -c1 * s2;
+//!
+//!     rmat[[2, 0]] = s2 * s3;
+//!     rmat[[2, 1]] = s2 * c3;
+//!     rmat[[2, 2]] = c2;
+//! });
+//!
+//! let eye2d = Array2::<f64>::eye(3);
+//!
+//! let mut rotated = Array3::<f64>::zeros((3, 3, nelems).f());
+//! azip!((mut rotated in rotated.axis_iter_mut(Axis(2)), rmat in rmat.axis_iter(Axis(2))) {
+//!     rotated.assign(&rmat.dot(&eye2d));
+//! });
 //! ```
diff --git a/src/doc/ndarray_for_numpy_users/mod.rs b/src/doc/ndarray_for_numpy_users/mod.rs
index bedd1c8b..f9c05c61 100644
--- a/src/doc/ndarray_for_numpy_users/mod.rs
+++ b/src/doc/ndarray_for_numpy_users/mod.rs
@@ -23,7 +23,7 @@
 //!
 //! # Similarities
 //!
-//! `ndarray`'s array type ([`ArrayBase`][ArrayBase]), is very similar to
+//! `ndarray`'s array type ([`ArrayBase`]), is very similar to
 //! NumPy's array type (`numpy.ndarray`):
 //!
 //! * Arrays have a single element type.
@@ -70,12 +70,12 @@
 //! </td>
 //! <td>
 //!
-//! In `ndarray`, all arrays are instances of [`ArrayBase`][ArrayBase], but
-//! `ArrayBase` is generic over the ownership of the data. [`Array`][Array]
-//! owns its data; [`ArrayView`][ArrayView] is a view;
-//! [`ArrayViewMut`][ArrayViewMut] is a mutable view; [`CowArray`][CowArray]
+//! In `ndarray`, all arrays are instances of [`ArrayBase`], but
+//! `ArrayBase` is generic over the ownership of the data. [`Array`]
+//! owns its data; [`ArrayView`] is a view;
+//! [`ArrayViewMut`] is a mutable view; [`CowArray`]
 //! either owns its data or is a view (with copy-on-write mutation of the view
-//! variant); and [`ArcArray`][ArcArray] has a reference-counted pointer to its
+//! variant); and [`ArcArray`] has a reference-counted pointer to its
 //! data (with copy-on-write mutation). Arrays and views follow Rust's aliasing
 //! rules.
 //!
@@ -91,7 +91,7 @@
 //! <td>
 //!
 //! In `ndarray`, you can create fixed-dimension arrays, such as
-//! [`Array2`][Array2]. This takes advantage of the type system to help you
+//! [`Array2`]. This takes advantage of the type system to help you
 //! write correct code and also avoids small heap allocations for the shape and
 //! strides.
 //!
@@ -263,7 +263,7 @@
 //! Note that [`a.shape()`][.shape()], [`a.dim()`][.dim()], and
 //! [`a.raw_dim()`][.raw_dim()] all return the shape of the array, but as
 //! different types. `a.shape()` returns the shape as `&[Ix]`, (where
-//! [`Ix`][Ix] is `usize`) which is useful for general operations on the shape.
+//! [`Ix`] is `usize`) which is useful for general operations on the shape.
 //! `a.dim()` returns the shape as `D::Pattern`, which is useful for
 //! pattern-matching shapes. `a.raw_dim()` returns the shape as `D`, which is
 //! useful for creating other arrays of the same shape.
@@ -376,7 +376,7 @@
 //!
 //! </td><td>
 //!
-//! [`a * b`, `a + b`, etc.](../../struct.ArrayBase.html#arithmetic-operations)
+//! [`a * b`, `a + b`, etc.](ArrayBase#arithmetic-operations)
 //!
 //! </td><td>
 //!
@@ -540,17 +540,17 @@
 //! ## Iteration
 //!
 //! `ndarray` has lots of interesting iterators/producers that implement the
-//! [`NdProducer`][NdProducer] trait, which is a generalization of `Iterator`
+//! [`NdProducer`](crate::NdProducer) trait, which is a generalization of `Iterator`
 //! to multiple dimensions. This makes it possible to correctly and efficiently
 //! zip together slices/subviews of arrays in multiple dimensions with
-//! [`Zip`][Zip] or [`azip!()`][azip!]. The purpose of this is similar to
+//! [`Zip`] or [`azip!()`]. The purpose of this is similar to
 //! [`np.nditer`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.nditer.html),
-//! but [`Zip`][Zip] is implemented and used somewhat differently.
+//! but [`Zip`] is implemented and used somewhat differently.
 //!
 //! This table lists some of the iterators/producers which have a direct
 //! equivalent in NumPy. For a more complete introduction to producers and
 //! iterators, see [*Loops, Producers, and
-//! Iterators*](../../struct.ArrayBase.html#loops-producers-and-iterators).
+//! Iterators*](ArrayBase#loops-producers-and-iterators).
 //! Note that there are also variants of these iterators (with a `_mut` suffix)
 //! that yield `ArrayViewMut` instead of `ArrayView`.
 //!
@@ -570,88 +570,77 @@
 //! `a[:,4]` | [`a.column(4)`][.column()] or [`a.column_mut(4)`][.column_mut()] | view (or mutable view) of column 4 in a 2-D array
 //! `a.shape[0] == a.shape[1]` | [`a.is_square()`][.is_square()] | check if the array is square
 //!
-//! [.abs_diff_eq()]: ../../struct.ArrayBase.html#impl-AbsDiffEq<ArrayBase<S2%2C%20D>>
-//! [ArcArray]: ../../type.ArcArray.html
-//! [arr2()]: ../../fn.arr2.html
-//! [array!]: ../../macro.array.html
-//! [Array]: ../../type.Array.html
-//! [Array2]: ../../type.Array2.html
-//! [ArrayBase]: ../../struct.ArrayBase.html
-//! [ArrayView]: ../../type.ArrayView.html
-//! [ArrayViewMut]: ../../type.ArrayViewMut.html
-//! [.assign()]: ../../struct.ArrayBase.html#method.assign
-//! [.axis_iter()]: ../../struct.ArrayBase.html#method.axis_iter
-//! [azip!]: ../../macro.azip.html
-//! [.ncols()]: ../../struct.ArrayBase.html#method.ncols
-//! [.column()]: ../../struct.ArrayBase.html#method.column
-//! [.column_mut()]: ../../struct.ArrayBase.html#method.column_mut
-//! [concatenate!]: ../../macro.concatenate.html
-//! [concatenate()]: ../../fn.concatenate.html
-//! [CowArray]: ../../type.CowArray.html
-//! [::default()]: ../../struct.ArrayBase.html#method.default
-//! [.diag()]: ../../struct.ArrayBase.html#method.diag
-//! [.dim()]: ../../struct.ArrayBase.html#method.dim
-//! [::eye()]: ../../struct.ArrayBase.html#method.eye
-//! [.fill()]: ../../struct.ArrayBase.html#method.fill
-//! [.fold()]: ../../struct.ArrayBase.html#method.fold
-//! [.fold_axis()]: ../../struct.ArrayBase.html#method.fold_axis
-//! [::from_elem()]: ../../struct.ArrayBase.html#method.from_elem
-//! [::from_iter()]: ../../struct.ArrayBase.html#method.from_iter
-//! [::from_diag()]: ../../struct.ArrayBase.html#method.from_diag
-//! [::from_shape_fn()]: ../../struct.ArrayBase.html#method.from_shape_fn
-//! [::from_shape_vec()]: ../../struct.ArrayBase.html#method.from_shape_vec
-//! [::from_shape_vec_unchecked()]: ../../struct.ArrayBase.html#method.from_shape_vec_unchecked
-//! [::from_vec()]: ../../struct.ArrayBase.html#method.from_vec
-//! [.index()]: ../../struct.ArrayBase.html#impl-Index<I>
-//! [.indexed_iter()]: ../../struct.ArrayBase.html#method.indexed_iter
-//! [.insert_axis()]: ../../struct.ArrayBase.html#method.insert_axis
-//! [.is_empty()]: ../../struct.ArrayBase.html#method.is_empty
-//! [.is_square()]: ../../struct.ArrayBase.html#method.is_square
-//! [.iter()]: ../../struct.ArrayBase.html#method.iter
-//! [Ix]: ../../type.Ix.html
-//! [.len()]: ../../struct.ArrayBase.html#method.len
-//! [.len_of()]: ../../struct.ArrayBase.html#method.len_of
-//! [::linspace()]: ../../struct.ArrayBase.html#method.linspace
-//! [::logspace()]: ../../struct.ArrayBase.html#method.logspace
-//! [::geomspace()]: ../../struct.ArrayBase.html#method.geomspace
-//! [.map()]: ../../struct.ArrayBase.html#method.map
-//! [.map_axis()]: ../../struct.ArrayBase.html#method.map_axis
-//! [.map_inplace()]: ../../struct.ArrayBase.html#method.map_inplace
-//! [.mapv()]: ../../struct.ArrayBase.html#method.mapv
-//! [.mapv_inplace()]: ../../struct.ArrayBase.html#method.mapv_inplace
-//! [.mapv_into()]: ../../struct.ArrayBase.html#method.mapv_into
-//! [matrix-* dot]: ../../struct.ArrayBase.html#method.dot-1
-//! [.mean()]: ../../struct.ArrayBase.html#method.mean
-//! [.mean_axis()]: ../../struct.ArrayBase.html#method.mean_axis
-//! [.ndim()]: ../../struct.ArrayBase.html#method.ndim
-//! [NdProducer]: ../../trait.NdProducer.html
-//! [::ones()]: ../../struct.ArrayBase.html#method.ones
-//! [.outer_iter()]: ../../struct.ArrayBase.html#method.outer_iter
-//! [::range()]: ../../struct.ArrayBase.html#method.range
-//! [.raw_dim()]: ../../struct.ArrayBase.html#method.raw_dim
-//! [.reversed_axes()]: ../../struct.ArrayBase.html#method.reversed_axes
-//! [.row()]: ../../struct.ArrayBase.html#method.row
-//! [.row_mut()]: ../../struct.ArrayBase.html#method.row_mut
-//! [.nrows()]: ../../struct.ArrayBase.html#method.nrows
-//! [s!]: ../../macro.s.html
-//! [.sum()]: ../../struct.ArrayBase.html#method.sum
-//! [.slice()]: ../../struct.ArrayBase.html#method.slice
-//! [.slice_axis()]: ../../struct.ArrayBase.html#method.slice_axis
-//! [.slice_collapse()]: ../../struct.ArrayBase.html#method.slice_collapse
-//! [.slice_move()]: ../../struct.ArrayBase.html#method.slice_move
-//! [.slice_mut()]: ../../struct.ArrayBase.html#method.slice_mut
-//! [.shape()]: ../../struct.ArrayBase.html#method.shape
-//! [stack!]: ../../macro.stack.html
-//! [stack()]: ../../fn.stack.html
-//! [.strides()]: ../../struct.ArrayBase.html#method.strides
-//! [.index_axis()]: ../../struct.ArrayBase.html#method.index_axis
-//! [.sum_axis()]: ../../struct.ArrayBase.html#method.sum_axis
-//! [.t()]: ../../struct.ArrayBase.html#method.t
-//! [vec-* dot]: ../../struct.ArrayBase.html#method.dot
-//! [.for_each()]: ../../struct.ArrayBase.html#method.for_each
-//! [::zeros()]: ../../struct.ArrayBase.html#method.zeros
-//! [Zip]: ../../struct.Zip.html
+//! [.abs_diff_eq()]: ArrayBase#impl-AbsDiffEq<ArrayBase<S2%2C%20D>>
+//! [.assign()]: ArrayBase::assign
+//! [.axis_iter()]: ArrayBase::axis_iter
+//! [.ncols()]: ArrayBase::ncols
+//! [.column()]: ArrayBase::column
+//! [.column_mut()]: ArrayBase::column_mut
+//! [concatenate()]: crate::concatenate()
+//! [::default()]: ArrayBase::default
+//! [.diag()]: ArrayBase::diag
+//! [.dim()]: ArrayBase::dim
+//! [::eye()]: ArrayBase::eye
+//! [.fill()]: ArrayBase::fill
+//! [.fold()]: ArrayBase::fold
+//! [.fold_axis()]: ArrayBase::fold_axis
+//! [::from_elem()]: ArrayBase::from_elem
+//! [::from_iter()]: ArrayBase::from_iter
+//! [::from_diag()]: ArrayBase::from_diag
+//! [::from_shape_fn()]: ArrayBase::from_shape_fn
+//! [::from_shape_vec()]: ArrayBase::from_shape_vec
+//! [::from_shape_vec_unchecked()]: ArrayBase::from_shape_vec_unchecked
+//! [::from_vec()]: ArrayBase::from_vec
+//! [.index()]: ArrayBase#impl-Index<I>
+//! [.indexed_iter()]: ArrayBase::indexed_iter
+//! [.insert_axis()]: ArrayBase::insert_axis
+//! [.is_empty()]: ArrayBase::is_empty
+//! [.is_square()]: ArrayBase::is_square
+//! [.iter()]: ArrayBase::iter
+//! [.len()]: ArrayBase::len
+//! [.len_of()]: ArrayBase::len_of
+//! [::linspace()]: ArrayBase::linspace
+//! [::logspace()]: ArrayBase::logspace
+//! [::geomspace()]: ArrayBase::geomspace
+//! [.map()]: ArrayBase::map
+//! [.map_axis()]: ArrayBase::map_axis
+//! [.map_inplace()]: ArrayBase::map_inplace
+//! [.mapv()]: ArrayBase::mapv
+//! [.mapv_inplace()]: ArrayBase::mapv_inplace
+//! [.mapv_into()]: ArrayBase::mapv_into
+//! [matrix-* dot]: ArrayBase::dot-1
+//! [.mean()]: ArrayBase::mean
+//! [.mean_axis()]: ArrayBase::mean_axis
+//! [.ndim()]: ArrayBase::ndim
+//! [::ones()]: ArrayBase::ones
+//! [.outer_iter()]: ArrayBase::outer_iter
+//! [::range()]: ArrayBase::range
+//! [.raw_dim()]: ArrayBase::raw_dim
+//! [.reversed_axes()]: ArrayBase::reversed_axes
+//! [.row()]: ArrayBase::row
+//! [.row_mut()]: ArrayBase::row_mut
+//! [.nrows()]: ArrayBase::nrows
+//! [.sum()]: ArrayBase::sum
+//! [.slice()]: ArrayBase::slice
+//! [.slice_axis()]: ArrayBase::slice_axis
+//! [.slice_collapse()]: ArrayBase::slice_collapse
+//! [.slice_move()]: ArrayBase::slice_move
+//! [.slice_mut()]: ArrayBase::slice_mut
+//! [.shape()]: ArrayBase::shape
+//! [stack()]: crate::stack()
+//! [.strides()]: ArrayBase::strides
+//! [.index_axis()]: ArrayBase::index_axis
+//! [.sum_axis()]: ArrayBase::sum_axis
+//! [.t()]: ArrayBase::t
+//! [vec-* dot]: ArrayBase::dot
+//! [.for_each()]: ArrayBase::for_each
+//! [::zeros()]: ArrayBase::zeros
+//! [`Zip`]: crate::Zip
 
 pub mod coord_transform;
 pub mod rk_step;
 pub mod simple_math;
+
+// This is to avoid putting `crate::` everywhere
+#[allow(unused_imports)]
+use crate::imp_prelude::*;
diff --git a/src/doc/ndarray_for_numpy_users/rk_step.rs b/src/doc/ndarray_for_numpy_users/rk_step.rs
index 92935f3e..0448e070 100644
--- a/src/doc/ndarray_for_numpy_users/rk_step.rs
+++ b/src/doc/ndarray_for_numpy_users/rk_step.rs
@@ -120,10 +120,10 @@
 //! * Don't return a newly allocated `f_new` array. If the caller wants this
 //!   information, they can get it from the last row of `k`.
 //!
-//! * Use [`c.mul_add(h, t)`][f64.mul_add()] instead of `t + c * h`. This is
+//! * Use [`c.mul_add(h, t)`](f64::mul_add) instead of `t + c * h`. This is
 //!   faster and reduces the floating-point error. It might also be beneficial
-//!   to use [`.scaled_add()`][.scaled_add()] or a combination of
-//!   [`azip!()`][azip!] and [`.mul_add()`][f64.mul_add()] on the arrays in
+//!   to use [`.scaled_add()`] or a combination of
+//!   [`azip!()`] and [`.mul_add()`](f64::mul_add) on the arrays in
 //!   some places, but that's not demonstrated in the example below.
 //!
 //! ```
@@ -168,9 +168,7 @@
 //! # fn main() { let _ = rk_step::<fn(_, ArrayView1<'_, f64>, ArrayViewMut1<'_, f64>)>; }
 //! ```
 //!
-//! [f64.mul_add()]: https://doc.rust-lang.org/std/primitive.f64.html#method.mul_add
-//! [.scaled_add()]: ../../../struct.ArrayBase.html#method.scaled_add
-//! [azip!]: ../../../macro.azip.html
+//! [`.scaled_add()`]: crate::ArrayBase::scaled_add
 //!
 //! ### SciPy license
 //!
diff --git a/src/free_functions.rs b/src/free_functions.rs
index 95eb7dc8..2b30e0bd 100644
--- a/src/free_functions.rs
+++ b/src/free_functions.rs
@@ -14,7 +14,7 @@ use alloc::vec::Vec;
 use crate::imp_prelude::*;
 use crate::{dimension, ArcArray1, ArcArray2};
 
-/// Create an [**`Array`**](type.Array.html) with one, two or
+/// Create an **[`Array`]** with one, two or
 /// three dimensions.
 ///
 /// ```
@@ -169,6 +169,7 @@ pub fn aview_mut2<A, V: FixedInitializer<Elem = A>>(xs: &mut [V]) -> ArrayViewMu
 }
 
 /// Fixed-size array used for array initialization
+#[allow(clippy::missing_safety_doc)] // Should not be implemented downstream and to be deprecated.
 pub unsafe trait FixedInitializer {
     type Elem;
     fn as_init_slice(&self) -> &[Self::Elem];
diff --git a/src/impl_constructors.rs b/src/impl_constructors.rs
index 46d0c76c..df5e295c 100644
--- a/src/impl_constructors.rs
+++ b/src/impl_constructors.rs
@@ -227,9 +227,32 @@ where
     {
         let n = diag.len();
         let mut arr = Self::zeros((n, n));
-        arr.diag_mut().assign(&diag);
+        arr.diag_mut().assign(diag);
         arr
     }
+
+    /// Create a square 2D matrix of the specified size, with the specified
+    /// element along the diagonal and zeros elsewhere.
+    ///
+    /// **Panics** if `n * n` would overflow `isize`.
+    ///
+    /// ```rust
+    /// use ndarray::{array, Array2};
+    ///
+    /// let array = Array2::from_diag_elem(2, 5.);
+    /// assert_eq!(array, array![[5., 0.], [0., 5.]]);
+    /// ```
+    pub fn from_diag_elem(n: usize, elem: A) -> Self
+    where
+        S: DataMut,
+        A: Clone + Zero,
+    {
+        let mut eye = Self::zeros((n, n));
+        for a_ii in eye.diag_mut() {
+            *a_ii = elem.clone();
+        }
+        eye
+    }
 }
 
 #[cfg(not(debug_assertions))]
@@ -269,7 +292,7 @@ macro_rules! size_of_shape_checked_unwrap {
 /// column major (“f” order) memory layout instead of the default row major.
 /// For example `Array::zeros((5, 6).f())` makes a column major 5 × 6 array.
 ///
-/// Use [`IxDyn`](type.IxDyn.html) for the shape to create an array with dynamic
+/// Use [`type@IxDyn`] for the shape to create an array with dynamic
 /// number of axes.
 ///
 /// Finally, the few constructors that take a completely general
diff --git a/src/impl_cow.rs b/src/impl_cow.rs
index b880fd62..5cd00077 100644
--- a/src/impl_cow.rs
+++ b/src/impl_cow.rs
@@ -11,8 +11,6 @@ use crate::imp_prelude::*;
 /// Methods specific to `CowArray`.
 ///
 /// ***See also all methods for [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
 impl<'a, A, D> CowArray<'a, A, D>
 where
     D: Dimension,
diff --git a/src/impl_methods.rs b/src/impl_methods.rs
index 28098f68..745a8e60 100644
--- a/src/impl_methods.rs
+++ b/src/impl_methods.rs
@@ -179,7 +179,7 @@ where
     /// 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.
+    /// desired memory layout and [`.assign()`](Self::assign) the data.
     /// Alternatively, you can collectan iterator, like this for a result in
     /// standard layout:
     ///
@@ -217,7 +217,7 @@ where
                 )
             }
         } else {
-            self.map(|x| x.clone())
+            self.map(A::clone)
         }
     }
 
@@ -241,6 +241,34 @@ where
         S::into_owned(self)
     }
 
+    /// Converts the array into `Array<A, D>` if this is possible without
+    /// cloning the array elements. Otherwise, returns `self` unchanged.
+    ///
+    /// ```
+    /// use ndarray::{array, rcarr2, ArcArray2, Array2};
+    ///
+    /// // Reference-counted, clone-on-write `ArcArray`.
+    /// let a: ArcArray2<_> = rcarr2(&[[1., 2.], [3., 4.]]);
+    /// {
+    ///     // Another reference to the same data.
+    ///     let b: ArcArray2<_> = a.clone();
+    ///     // Since there are two references to the same data, `.into_owned()`
+    ///     // would require cloning the data, so `.try_into_owned_nocopy()`
+    ///     // returns `Err`.
+    ///     assert!(b.try_into_owned_nocopy().is_err());
+    /// }
+    /// // Here, since the second reference has been dropped, the `ArcArray`
+    /// // can be converted into an `Array` without cloning the data.
+    /// let unique: Array2<_> = a.try_into_owned_nocopy().unwrap();
+    /// assert_eq!(unique, array![[1., 2.], [3., 4.]]);
+    /// ```
+    pub fn try_into_owned_nocopy(self) -> Result<Array<A, D>, Self>
+    where
+        S: Data,
+    {
+        S::try_into_owned_nocopy(self)
+    }
+
     /// Turn the array into a shared ownership (copy on write) array,
     /// without any copying.
     pub fn into_shared(self) -> ArcArray<A, D>
@@ -256,6 +284,19 @@ where
 
     /// Returns a reference to the first element of the array, or `None` if it
     /// is empty.
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use ndarray::Array3;
+    ///
+    /// let mut a = Array3::<f64>::zeros([3, 4, 2]);
+    /// a[[0, 0, 0]] = 42.;
+    /// assert_eq!(a.first(), Some(&42.));
+    ///
+    /// let b = Array3::<f64>::zeros([3, 0, 5]);
+    /// assert_eq!(b.first(), None);
+    /// ```
     pub fn first(&self) -> Option<&A>
     where
         S: Data,
@@ -269,6 +310,19 @@ where
 
     /// Returns a mutable reference to the first element of the array, or
     /// `None` if it is empty.
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use ndarray::Array3;
+    ///
+    /// let mut a = Array3::<f64>::zeros([3, 4, 2]);
+    /// *a.first_mut().unwrap() = 42.;
+    /// assert_eq!(a[[0, 0, 0]], 42.);
+    ///
+    /// let mut b = Array3::<f64>::zeros([3, 0, 5]);
+    /// assert_eq!(b.first_mut(), None);
+    /// ```
     pub fn first_mut(&mut self) -> Option<&mut A>
     where
         S: DataMut,
@@ -280,6 +334,66 @@ where
         }
     }
 
+    /// Returns a reference to the last element of the array, or `None` if it
+    /// is empty.
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use ndarray::Array3;
+    ///
+    /// let mut a = Array3::<f64>::zeros([3, 4, 2]);
+    /// a[[2, 3, 1]] = 42.;
+    /// assert_eq!(a.last(), Some(&42.));
+    ///
+    /// let b = Array3::<f64>::zeros([3, 0, 5]);
+    /// assert_eq!(b.last(), None);
+    /// ```
+    pub fn last(&self) -> Option<&A>
+    where
+        S: Data,
+    {
+        if self.is_empty() {
+            None
+        } else {
+            let mut index = self.raw_dim();
+            for ax in 0..index.ndim() {
+                index[ax] -= 1;
+            }
+            Some(unsafe { self.uget(index) })
+        }
+    }
+
+    /// Returns a mutable reference to the last element of the array, or `None`
+    /// if it is empty.
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use ndarray::Array3;
+    ///
+    /// let mut a = Array3::<f64>::zeros([3, 4, 2]);
+    /// *a.last_mut().unwrap() = 42.;
+    /// assert_eq!(a[[2, 3, 1]], 42.);
+    ///
+    /// let mut b = Array3::<f64>::zeros([3, 0, 5]);
+    /// assert_eq!(b.last_mut(), None);
+    /// ```
+    pub fn last_mut(&mut self) -> Option<&mut A>
+    where
+        S: DataMut,
+    {
+        if self.is_empty() {
+            None
+        } else {
+            let mut index = self.raw_dim();
+            for ax in 0..index.ndim() {
+                index[ax] -= 1;
+            }
+            Some(unsafe { self.uget_mut(index) })
+        }
+    }
+
     /// Return an iterator of references to the elements of the array.
     ///
     /// Elements are visited in the *logical order* of the array, which
@@ -314,7 +428,7 @@ where
     ///
     /// Iterator element type is `(D::Pattern, &A)`.
     ///
-    /// See also [`Zip::indexed`](struct.Zip.html)
+    /// See also [`Zip::indexed`]
     pub fn indexed_iter(&self) -> IndexedIter<'_, A, D>
     where
         S: Data,
@@ -458,9 +572,9 @@ where
     /// collapsed, as in [`.collapse_axis()`], rather than removed, as in
     /// [`.slice_move()`] or [`.index_axis_move()`].
     ///
-    /// [`.collapse_axis()`]: #method.collapse_axis
-    /// [`.slice_move()`]: #method.slice_move
-    /// [`.index_axis_move()`]: #method.index_axis_move
+    /// [`.collapse_axis()`]: Self::collapse_axis
+    /// [`.slice_move()`]: Self::slice_move
+    /// [`.index_axis_move()`]: Self::index_axis_move
     ///
     /// See [*Slicing*](#slicing) for full documentation.
     /// See also [`s!`], [`SliceArg`], and [`SliceInfo`](crate::SliceInfo).
@@ -830,7 +944,7 @@ where
 
     /// Collapses the array to `index` along the axis and removes the axis.
     ///
-    /// See [`.index_axis()`](#method.index_axis) and [*Subviews*](#subviews) for full documentation.
+    /// See [`.index_axis()`](Self::index_axis) and [*Subviews*](#subviews) for full documentation.
     ///
     /// **Panics** if `axis` or `index` is out of bounds.
     pub fn index_axis_move(mut self, axis: Axis, index: usize) -> ArrayBase<S, D::Smaller>
@@ -1313,6 +1427,50 @@ where
         Windows::new(self.view(), window_size)
     }
 
+    /// Returns a producer which traverses over all windows of a given length along an axis.
+    ///
+    /// The windows are all distinct, possibly-overlapping views. The shape of each window
+    /// is the shape of `self`, with the length of `axis` replaced with `window_size`.
+    ///
+    /// **Panics** if `axis` is out-of-bounds or if `window_size` is zero.
+    ///
+    /// ```
+    /// use ndarray::{Array3, Axis, s};
+    ///
+    /// let arr = Array3::from_shape_fn([4, 5, 2], |(i, j, k)| i * 100 + j * 10 + k);
+    /// let correct = vec![
+    ///     arr.slice(s![.., 0..3, ..]),
+    ///     arr.slice(s![.., 1..4, ..]),
+    ///     arr.slice(s![.., 2..5, ..]),
+    /// ];
+    /// for (window, correct) in arr.axis_windows(Axis(1), 3).into_iter().zip(&correct) {
+    ///     assert_eq!(window, correct);
+    ///     assert_eq!(window.shape(), &[4, 3, 2]);
+    /// }
+    /// ```
+    pub fn axis_windows(&self, axis: Axis, window_size: usize) -> Windows<'_, A, D>
+    where
+        S: Data,
+    {
+        let axis_index = axis.index();
+
+        ndassert!(
+            axis_index < self.ndim(),
+            concat!(
+                "Window axis {} does not match array dimension {} ",
+                "(with array of shape {:?})"
+            ),
+            axis_index,
+            self.ndim(),
+            self.shape()
+        );
+
+        let mut size = self.raw_dim();
+        size[axis_index] = window_size;
+
+        Windows::new(self.view(), size)
+    }
+
     // Return (length, stride) for diagonal
     fn diag_params(&self) -> (Ix, Ixs) {
         /* empty shape has len 1 */
@@ -1446,6 +1604,15 @@ where
     }
 
     /// Return a mutable pointer to the first element in the array.
+    ///
+    /// This method attempts to unshare the data. If `S: DataMut`, then the
+    /// data is guaranteed to be uniquely held on return.
+    ///
+    /// # Warning
+    ///
+    /// When accessing elements through this pointer, make sure to use strides
+    /// obtained *after* calling this method, since the process of unsharing
+    /// the data may change the strides.
     #[inline(always)]
     pub fn as_mut_ptr(&mut self) -> *mut A
     where
@@ -1462,6 +1629,9 @@ where
     }
 
     /// Return a raw mutable view of the array.
+    ///
+    /// This method attempts to unshare the data. If `S: DataMut`, then the
+    /// data is guaranteed to be uniquely held on return.
     #[inline]
     pub fn raw_view_mut(&mut self) -> RawArrayViewMut<A, D>
     where
@@ -1536,6 +1706,10 @@ where
 
     /// Return the array’s data as a slice if it is contiguous,
     /// return `None` otherwise.
+    ///
+    /// In the contiguous case, in order to return a unique reference, this
+    /// method unshares the data if necessary, but it preserves the existing
+    /// strides.
     pub fn as_slice_memory_order_mut(&mut self) -> Option<&mut [A]>
     where
         S: DataMut,
@@ -2190,7 +2364,7 @@ where
             if let Some(self_s) = self.as_slice_memory_order_mut() {
                 if let Some(rhs_s) = rhs.as_slice_memory_order() {
                     for (s, r) in self_s.iter_mut().zip(rhs_s) {
-                        f(s, &r);
+                        f(s, r);
                     }
                     return;
                 }
@@ -2371,6 +2545,48 @@ where
         self
     }
 
+    /// Consume the array, call `f` by **v**alue on each element, and return an
+    /// owned array with the new values. Works for **any** `F: FnMut(A)->B`.
+    ///
+    /// If `A` and `B` are the same type then the map is performed by delegating
+    /// to [`mapv_into`] and then converting into an owned array. This avoids
+    /// unnecessary memory allocations in [`mapv`].
+    ///
+    /// If `A` and `B` are different types then a new array is allocated and the
+    /// map is performed as in [`mapv`].
+    ///
+    /// Elements are visited in arbitrary order.
+    ///
+    /// [`mapv_into`]: ArrayBase::mapv_into
+    /// [`mapv`]: ArrayBase::mapv
+    pub fn mapv_into_any<B, F>(self, mut f: F) -> Array<B, D>
+    where
+        S: DataMut,
+        F: FnMut(A) -> B,
+        A: Clone + 'static,
+        B: 'static,
+    {
+        if core::any::TypeId::of::<A>() == core::any::TypeId::of::<B>() {
+            // A and B are the same type.
+            // Wrap f in a closure of type FnMut(A) -> A .
+            let f = |a| {
+                let b = f(a);
+                // Safe because A and B are the same type.
+                unsafe { unlimited_transmute::<B, A>(b) }
+            };
+            // Delegate to mapv_into() using the wrapped closure.
+            // Convert output to a uniquely owned array of type Array<A, D>.
+            let output = self.mapv_into(f).into_owned();
+            // Change the return type from Array<A, D> to Array<B, D>.
+            // Again, safe because A and B are the same type.
+            unsafe { unlimited_transmute::<Array<A, D>, Array<B, D>>(output) }
+        } else {
+            // A and B are not the same type.
+            // Fallback to mapv().
+            self.mapv(f)
+        }
+    }
+
     /// Modify the array in place by calling `f` by mutable reference on each element.
     ///
     /// Elements are visited in arbitrary order.
diff --git a/src/impl_owned_array.rs b/src/impl_owned_array.rs
index efe67a04..e5893274 100644
--- a/src/impl_owned_array.rs
+++ b/src/impl_owned_array.rs
@@ -17,8 +17,6 @@ use crate::Zip;
 /// Methods specific to `Array0`.
 ///
 /// ***See also all methods for [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
 impl<A> Array<A, Ix0> {
     /// Returns the single element in the array without cloning it.
     ///
@@ -56,8 +54,6 @@ impl<A> Array<A, Ix0> {
 /// Methods specific to `Array`.
 ///
 /// ***See also all methods for [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
 impl<A, D> Array<A, D>
 where
     D: Dimension,
@@ -75,8 +71,6 @@ where
 /// Methods specific to `Array2`.
 ///
 /// ***See also all methods for [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
 impl<A> Array<A, Ix2> {
     /// Append a row to an array
     ///
diff --git a/src/impl_raw_views.rs b/src/impl_raw_views.rs
index 2ac5c08c..90cfa637 100644
--- a/src/impl_raw_views.rs
+++ b/src/impl_raw_views.rs
@@ -75,7 +75,7 @@ where
             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();
+                dimension::max_abs_offset_check_overflow::<A, _>(&dim, strides).unwrap();
             } else {
                 dimension::size_of_shape_checked(&dim).unwrap();
             }
@@ -217,7 +217,7 @@ where
             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();
+                dimension::max_abs_offset_check_overflow::<A, _>(&dim, strides).unwrap();
             } else {
                 dimension::size_of_shape_checked(&dim).unwrap();
             }
diff --git a/src/impl_special_element_types.rs b/src/impl_special_element_types.rs
index bf4384e4..5d5f1849 100644
--- a/src/impl_special_element_types.rs
+++ b/src/impl_special_element_types.rs
@@ -15,8 +15,6 @@ use crate::RawDataSubst;
 /// Methods specific to arrays with `MaybeUninit` elements.
 ///
 /// ***See also all methods for [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
 impl<A, S, D> ArrayBase<S, D>
 where
     S: RawDataSubst<A, Elem=MaybeUninit<A>>,
diff --git a/src/impl_views/conversions.rs b/src/impl_views/conversions.rs
index d0f91b22..7498add4 100644
--- a/src/impl_views/conversions.rs
+++ b/src/impl_views/conversions.rs
@@ -53,9 +53,6 @@ where
 /// Methods specific to `ArrayView0`.
 ///
 /// ***See also all methods for [`ArrayView`] and [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
-/// [`ArrayView`]: struct.ArrayView.html
 impl<'a, A> ArrayView<'a, A, Ix0> {
     /// Consume the view and return a reference to the single element in the array.
     ///
@@ -82,9 +79,6 @@ impl<'a, A> ArrayView<'a, A, Ix0> {
 /// Methods specific to `ArrayViewMut0`.
 ///
 /// ***See also all methods for [`ArrayViewMut`] and [`ArrayBase`]***
-///
-/// [`ArrayBase`]: struct.ArrayBase.html
-/// [`ArrayViewMut`]: struct.ArrayViewMut.html
 impl<'a, A> ArrayViewMut<'a, A, Ix0> {
     /// Consume the mutable view and return a mutable reference to the single element in the array.
     ///
diff --git a/src/impl_views/indexing.rs b/src/impl_views/indexing.rs
index d8ed956f..c73e870e 100644
--- a/src/impl_views/indexing.rs
+++ b/src/impl_views/indexing.rs
@@ -59,7 +59,7 @@ pub trait IndexLonger<I> {
     /// See also [the `get` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.get
+    /// [1]: ArrayBase::get
     ///
     /// **Panics** if index is out of bounds.
     fn index(self, index: I) -> Self::Output;
@@ -73,8 +73,8 @@ pub trait IndexLonger<I> {
     /// See also [the `get` method][1] (and [`get_mut`][2]) which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.get
-    /// [2]: struct.ArrayBase.html#method.get_mut
+    /// [1]: ArrayBase::get
+    /// [2]: ArrayBase::get_mut
     ///
     /// **Panics** if index is out of bounds.
     fn get(self, index: I) -> Option<Self::Output>;
@@ -87,7 +87,7 @@ pub trait IndexLonger<I> {
     /// See also [the `uget` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.uget
+    /// [1]: ArrayBase::uget
     ///
     /// **Note:** only unchecked for non-debug builds of ndarray.
     ///
@@ -113,7 +113,7 @@ where
     /// See also [the `get` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.get
+    /// [1]: ArrayBase::get
     ///
     /// **Panics** if index is out of bounds.
     fn index(self, index: I) -> &'a A {
@@ -133,7 +133,7 @@ where
     /// See also [the `uget` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.uget
+    /// [1]: ArrayBase::uget
     ///
     /// **Note:** only unchecked for non-debug builds of ndarray.
     unsafe fn uget(self, index: I) -> &'a A {
@@ -158,7 +158,7 @@ where
     /// See also [the `get_mut` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.get_mut
+    /// [1]: ArrayBase::get_mut
     ///
     /// **Panics** if index is out of bounds.
     fn index(mut self, index: I) -> &'a mut A {
@@ -177,7 +177,7 @@ where
     /// See also [the `get_mut` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.get_mut
+    /// [1]: ArrayBase::get_mut
     ///
     fn get(mut self, index: I) -> Option<&'a mut A> {
         debug_bounds_check!(self, index);
@@ -195,7 +195,7 @@ where
     /// See also [the `uget_mut` method][1] which works for all arrays and array
     /// views.
     ///
-    /// [1]: struct.ArrayBase.html#method.uget_mut
+    /// [1]: ArrayBase::uget_mut
     ///
     /// **Note:** only unchecked for non-debug builds of ndarray.
     unsafe fn uget(mut self, index: I) -> &'a mut A {
diff --git a/src/impl_views/splitting.rs b/src/impl_views/splitting.rs
index a36ae4dd..f2c0dc41 100644
--- a/src/impl_views/splitting.rs
+++ b/src/impl_views/splitting.rs
@@ -121,7 +121,7 @@ where
     /// [`MultiSliceArg`], [`s!`], [`SliceArg`](crate::SliceArg), and
     /// [`SliceInfo`](crate::SliceInfo).
     ///
-    /// [`.multi_slice_mut()`]: struct.ArrayBase.html#method.multi_slice_mut
+    /// [`.multi_slice_mut()`]: ArrayBase::multi_slice_mut
     ///
     /// **Panics** if any of the following occur:
     ///
diff --git a/src/iterators/chunks.rs b/src/iterators/chunks.rs
index e41c1bf2..ba0e789f 100644
--- a/src/iterators/chunks.rs
+++ b/src/iterators/chunks.rs
@@ -28,7 +28,7 @@ type BaseProducerMut<'a, A, D> = ArrayViewMut<'a, A, D>;
 
 /// Exact chunks producer and iterable.
 ///
-/// See [`.exact_chunks()`](../struct.ArrayBase.html#method.exact_chunks) for more
+/// See [`.exact_chunks()`](ArrayBase::exact_chunks) for more
 /// information.
 //#[derive(Debug)]
 pub struct ExactChunks<'a, A, D> {
@@ -88,7 +88,7 @@ where
 
 /// Exact chunks iterator.
 ///
-/// See [`.exact_chunks()`](../struct.ArrayBase.html#method.exact_chunks) for more
+/// See [`.exact_chunks()`](ArrayBase::exact_chunks) for more
 /// information.
 pub struct ExactChunksIter<'a, A, D> {
     iter: ElementsBase<'a, A, D>,
@@ -118,7 +118,7 @@ impl_ndproducer! {
 
 /// Exact chunks producer and iterable.
 ///
-/// See [`.exact_chunks_mut()`](../struct.ArrayBase.html#method.exact_chunks_mut)
+/// See [`.exact_chunks_mut()`](ArrayBase::exact_chunks_mut)
 /// for more information.
 //#[derive(Debug)]
 pub struct ExactChunksMut<'a, A, D> {
@@ -222,7 +222,7 @@ impl_iterator! {
 
 /// Exact chunks iterator.
 ///
-/// See [`.exact_chunks_mut()`](../struct.ArrayBase.html#method.exact_chunks_mut)
+/// See [`.exact_chunks_mut()`](ArrayBase::exact_chunks_mut)
 /// for more information.
 pub struct ExactChunksIterMut<'a, A, D> {
     iter: ElementsBaseMut<'a, A, D>,
diff --git a/src/iterators/iter.rs b/src/iterators/iter.rs
index 7352e5e1..7b6d80cd 100644
--- a/src/iterators/iter.rs
+++ b/src/iterators/iter.rs
@@ -4,7 +4,7 @@
 //! implementation structs.
 //!
 //!
-//! See also [`NdProducer`](../trait.NdProducer.html).
+//! See also [`NdProducer`](crate::NdProducer).
 
 pub use crate::dimension::Axes;
 pub use crate::indexes::{Indices, IndicesIter};
diff --git a/src/iterators/lanes.rs b/src/iterators/lanes.rs
index 2163c58a..cda783a8 100644
--- a/src/iterators/lanes.rs
+++ b/src/iterators/lanes.rs
@@ -23,7 +23,7 @@ impl_ndproducer! {
     }
 }
 
-/// See [`.lanes()`](../struct.ArrayBase.html#method.lanes)
+/// See [`.lanes()`](ArrayBase::lanes)
 /// for more information.
 pub struct Lanes<'a, A, D> {
     base: ArrayView<'a, A, D>,
@@ -92,7 +92,7 @@ where
     }
 }
 
-/// See [`.lanes_mut()`](../struct.ArrayBase.html#method.lanes_mut)
+/// See [`.lanes_mut()`](ArrayBase::lanes_mut)
 /// for more information.
 pub struct LanesMut<'a, A, D> {
     base: ArrayViewMut<'a, A, D>,
diff --git a/src/iterators/mod.rs b/src/iterators/mod.rs
index bb618e5b..b2e4b9c6 100644
--- a/src/iterators/mod.rs
+++ b/src/iterators/mod.rs
@@ -315,7 +315,7 @@ pub enum ElementsRepr<S, C> {
 ///
 /// Iterator element type is `&'a A`.
 ///
-/// See [`.iter()`](../struct.ArrayBase.html#method.iter) for more information.
+/// See [`.iter()`](ArrayBase::iter) for more information.
 pub struct Iter<'a, A, D> {
     inner: ElementsRepr<SliceIter<'a, A>, ElementsBase<'a, A, D>>,
 }
@@ -330,7 +330,7 @@ pub struct ElementsBase<'a, A, D> {
 ///
 /// Iterator element type is `&'a mut A`.
 ///
-/// See [`.iter_mut()`](../struct.ArrayBase.html#method.iter_mut) for more information.
+/// See [`.iter_mut()`](ArrayBase::iter_mut) for more information.
 pub struct IterMut<'a, A, D> {
     inner: ElementsRepr<SliceIterMut<'a, A>, ElementsBaseMut<'a, A, D>>,
 }
@@ -354,12 +354,12 @@ impl<'a, A, D: Dimension> ElementsBaseMut<'a, A, D> {
 
 /// An iterator over the indexes and elements of an array.
 ///
-/// See [`.indexed_iter()`](../struct.ArrayBase.html#method.indexed_iter) for more information.
+/// See [`.indexed_iter()`](ArrayBase::indexed_iter) for more information.
 #[derive(Clone)]
 pub struct IndexedIter<'a, A, D>(ElementsBase<'a, A, D>);
 /// An iterator over the indexes and elements of an array (mutable).
 ///
-/// See [`.indexed_iter_mut()`](../struct.ArrayBase.html#method.indexed_iter_mut) for more information.
+/// See [`.indexed_iter_mut()`](ArrayBase::indexed_iter_mut) for more information.
 pub struct IndexedIterMut<'a, A, D>(ElementsBaseMut<'a, A, D>);
 
 impl<'a, A, D> IndexedIter<'a, A, D>
@@ -681,7 +681,7 @@ where
 /// An iterator that traverses over all axes but one, and yields a view for
 /// each lane along that axis.
 ///
-/// See [`.lanes()`](../struct.ArrayBase.html#method.lanes) for more information.
+/// See [`.lanes()`](ArrayBase::lanes) for more information.
 pub struct LanesIter<'a, A, D> {
     inner_len: Ix,
     inner_stride: Ixs,
@@ -732,7 +732,7 @@ where
 /// An iterator that traverses over all dimensions but the innermost,
 /// and yields each inner row (mutable).
 ///
-/// See [`.lanes_mut()`](../struct.ArrayBase.html#method.lanes_mut)
+/// See [`.lanes_mut()`](ArrayBase::lanes_mut)
 /// for more information.
 pub struct LanesIterMut<'a, A, D> {
     inner_len: Ix,
@@ -927,8 +927,8 @@ where
 ///
 /// Iterator element type is `ArrayView<'a, A, D>`.
 ///
-/// See [`.outer_iter()`](../struct.ArrayBase.html#method.outer_iter)
-/// or [`.axis_iter()`](../struct.ArrayBase.html#method.axis_iter)
+/// See [`.outer_iter()`](ArrayBase::outer_iter)
+/// or [`.axis_iter()`](ArrayBase::axis_iter)
 /// for more information.
 #[derive(Debug)]
 pub struct AxisIter<'a, A, D> {
@@ -1024,8 +1024,8 @@ where
 ///
 /// Iterator element type is `ArrayViewMut<'a, A, D>`.
 ///
-/// See [`.outer_iter_mut()`](../struct.ArrayBase.html#method.outer_iter_mut)
-/// or [`.axis_iter_mut()`](../struct.ArrayBase.html#method.axis_iter_mut)
+/// See [`.outer_iter_mut()`](ArrayBase::outer_iter_mut)
+/// or [`.axis_iter_mut()`](ArrayBase::axis_iter_mut)
 /// for more information.
 pub struct AxisIterMut<'a, A, D> {
     iter: AxisIterCore<A, D>,
@@ -1222,7 +1222,7 @@ impl<'a, A, D: Dimension> NdProducer for AxisIterMut<'a, A, D> {
 ///
 /// Iterator element type is `ArrayView<'a, A, D>`.
 ///
-/// See [`.axis_chunks_iter()`](../struct.ArrayBase.html#method.axis_chunks_iter) for more information.
+/// See [`.axis_chunks_iter()`](ArrayBase::axis_chunks_iter) for more information.
 pub struct AxisChunksIter<'a, A, D> {
     iter: AxisIterCore<A, D>,
     /// Index of the partial chunk (the chunk smaller than the specified chunk
@@ -1401,7 +1401,7 @@ macro_rules! chunk_iter_impl {
 ///
 /// Iterator element type is `ArrayViewMut<'a, A, D>`.
 ///
-/// See [`.axis_chunks_iter_mut()`](../struct.ArrayBase.html#method.axis_chunks_iter_mut)
+/// See [`.axis_chunks_iter_mut()`](ArrayBase::axis_chunks_iter_mut)
 /// for more information.
 pub struct AxisChunksIterMut<'a, A, D> {
     iter: AxisIterCore<A, D>,
@@ -1445,6 +1445,7 @@ send_sync_read_write!(ElementsBaseMut);
 ///
 /// The iterator must produce exactly the number of elements it reported or
 /// diverge before reaching the end.
+#[allow(clippy::missing_safety_doc)] // not nameable downstream
 pub unsafe trait TrustedIterator {}
 
 use crate::indexes::IndicesIterF;
diff --git a/src/iterators/windows.rs b/src/iterators/windows.rs
index 4538f7ab..c47bfece 100644
--- a/src/iterators/windows.rs
+++ b/src/iterators/windows.rs
@@ -6,7 +6,7 @@ use crate::NdProducer;
 
 /// Window producer and iterable
 ///
-/// See [`.windows()`](../struct.ArrayBase.html#method.windows) for more
+/// See [`.windows()`](ArrayBase::windows) for more
 /// information.
 pub struct Windows<'a, A, D> {
     base: ArrayView<'a, A, D>,
@@ -86,7 +86,7 @@ where
 
 /// Window iterator.
 ///
-/// See [`.windows()`](../struct.ArrayBase.html#method.windows) for more
+/// See [`.windows()`](ArrayBase::windows) for more
 /// information.
 pub struct WindowsIter<'a, A, D> {
     iter: ElementsBase<'a, A, D>,
diff --git a/src/itertools.rs b/src/itertools.rs
index 96732f90..4fa8145a 100644
--- a/src/itertools.rs
+++ b/src/itertools.rs
@@ -64,11 +64,9 @@ where
 /// The special cases of one and two arguments produce the equivalent of
 /// `$a.into_iter()` and `$a.into_iter().zip($b)` respectively.
 ///
-/// Prefer this macro `izip!()` over [`multizip`] for the performance benefits
+/// Prefer this macro `izip!()` over `multizip` for the performance benefits
 /// of using the standard library `.zip()`.
 ///
-/// [`multizip`]: fn.multizip.html
-///
 /// ```
 /// #[macro_use] extern crate itertools;
 /// # fn main() {
diff --git a/src/lib.rs b/src/lib.rs
index 8cf1a495..f773226f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -15,6 +15,8 @@
     clippy::manual_map, // is not an error
     clippy::while_let_on_iterator, // is not an error
     clippy::from_iter_instead_of_collect, // using from_iter is good style
+    clippy::if_then_panic, // is not an error
+    clippy::redundant_closure, // false positives clippy #7812
 )]
 #![doc(test(attr(deny(warnings))))]
 #![doc(test(attr(allow(unused_variables))))]
@@ -24,16 +26,16 @@
 //! The `ndarray` crate provides an *n*-dimensional container for general elements
 //! and for numerics.
 //!
-//! In *n*-dimensional we include for example 1-dimensional rows or columns,
+//! In *n*-dimensional we include, for example, 1-dimensional rows or columns,
 //! 2-dimensional matrices, and higher dimensional arrays. If the array has *n*
 //! dimensions, then an element in the array is accessed by using that many indices.
 //! Each dimension is also called an *axis*.
 //!
-//! - **[`ArrayBase`](struct.ArrayBase.html)**:
+//! - **[`ArrayBase`]**:
 //!   The *n*-dimensional array type itself.<br>
 //!   It is used to implement both the owned arrays and the views; see its docs
 //!   for an overview of all array features.<br>
-//! - The main specific array type is **[`Array`](type.Array.html)**, which owns
+//! - The main specific array type is **[`Array`]**, which owns
 //! its elements.
 //!
 //! ## Highlights
@@ -45,8 +47,8 @@
 //! - Higher order operations and arithmetic are performant
 //! - Array views can be used to slice and mutate any `[T]` data using
 //!   `ArrayView::from` and `ArrayViewMut::from`.
-//! - [`Zip`](struct.Zip.html) for lock step function application across two or more arrays or other
-//!   item producers ([`NdProducer`](trait.NdProducer.html) trait).
+//! - [`Zip`] for lock step function application across two or more arrays or other
+//!   item producers ([`NdProducer`] trait).
 //!
 //! ## Crate Status
 //!
@@ -99,13 +101,13 @@
 //!
 //! ## Documentation
 //!
-//! * The docs for [`ArrayBase`](struct.ArrayBase.html) provide an overview of
+//! * The docs for [`ArrayBase`] provide an overview of
 //!   the *n*-dimensional array type. Other good pages to look at are the
-//!   documentation for the [`s![]`](macro.s.html) and
-//!   [`azip!()`](macro.azip.html) macros.
+//!   documentation for the [`s![]`](s!) and
+//!   [`azip!()`](azip!) macros.
 //!
 //! * If you have experience with NumPy, you may also be interested in
-//!   [`ndarray_for_numpy_users`](doc/ndarray_for_numpy_users/index.html).
+//!   [`ndarray_for_numpy_users`](doc::ndarray_for_numpy_users).
 //!
 //! ## The ndarray ecosystem
 //!
@@ -258,12 +260,6 @@ pub type Ixs = isize;
 /// [`ArrayViewMut`] refer to `ArrayBase` with different types for the data
 /// container.
 ///
-/// [`Array`]: type.Array.html
-/// [`ArcArray`]: type.ArcArray.html
-/// [`ArrayView`]: type.ArrayView.html
-/// [`ArrayViewMut`]: type.ArrayViewMut.html
-/// [`CowArray`]: type.CowArray.html
-///
 /// ## Contents
 ///
 /// + [Array](#array)
@@ -286,7 +282,7 @@ pub type Ixs = isize;
 ///
 /// ## `Array`
 ///
-/// [`Array`](type.Array.html) is an owned array that owns the underlying array
+/// [`Array`] is an owned array that owns the underlying array
 /// elements directly (just like a `Vec`) and it is the default way to create and
 /// store n-dimensional data. `Array<A, D>` has two type parameters: `A` for
 /// the element type, and `D` for the dimensionality. A particular
@@ -305,17 +301,16 @@ pub type Ixs = isize;
 ///
 /// ## `ArcArray`
 ///
-/// [`ArcArray`](type.ArcArray.html) is an owned array with reference counted
+/// [`ArcArray`] is an owned array with reference counted
 /// data (shared ownership).
 /// Sharing requires that it uses copy-on-write for mutable operations.
 /// Calling a method for mutating elements on `ArcArray`, for example
-/// [`view_mut()`](#method.view_mut) or [`get_mut()`](#method.get_mut),
+/// [`view_mut()`](Self::view_mut) or [`get_mut()`](Self::get_mut),
 /// will break sharing and require a clone of the data (if it is not uniquely held).
 ///
 /// ## `CowArray`
 ///
-/// [`CowArray`](type.CowArray.html) is analogous to
-/// [`std::borrow::Cow`](https://doc.rust-lang.org/std/borrow/enum.Cow.html).
+/// [`CowArray`] is analogous to [`std::borrow::Cow`].
 /// It can represent either an immutable view or a uniquely owned array. If a
 /// `CowArray` instance is the immutable view variant, then calling a method
 /// for mutating elements in the array will cause it to be converted into the
@@ -380,15 +375,15 @@ pub type Ixs = isize;
 ///
 /// Important traits and types for dimension and indexing:
 ///
-/// - A [`Dim`](struct.Dim.html) value represents a dimensionality or index.
-/// - Trait [`Dimension`](trait.Dimension.html) is implemented by all
+/// - A [`struct@Dim`] value represents a dimensionality or index.
+/// - Trait [`Dimension`] is implemented by all
 /// dimensionalities. It defines many operations for dimensions and indices.
-/// - Trait [`IntoDimension`](trait.IntoDimension.html) is used to convert into a
+/// - Trait [`IntoDimension`] is used to convert into a
 /// `Dim` value.
-/// - Trait [`ShapeBuilder`](trait.ShapeBuilder.html) is an extension of
+/// - Trait [`ShapeBuilder`] is an extension of
 /// `IntoDimension` and is used when constructing an array. A shape describes
 /// not just the extent of each axis but also their strides.
-/// - Trait [`NdIndex`](trait.NdIndex.html) is an extension of `Dimension` and is
+/// - Trait [`NdIndex`] is an extension of `Dimension` and is
 /// for values that can be used with indexing syntax.
 ///
 ///
@@ -403,10 +398,10 @@ pub type Ixs = isize;
 ///
 /// ## Loops, Producers and Iterators
 ///
-/// Using [`Zip`](struct.Zip.html) is the most general way to apply a procedure
+/// Using [`Zip`] is the most general way to apply a procedure
 /// across one or several arrays or *producers*.
 ///
-/// [`NdProducer`](trait.NdProducer.html) is like an iterable but for
+/// [`NdProducer`] is like an iterable but for
 /// multidimensional data. All producers have dimensions and axes, like an
 /// array view, and they can be split and used with parallelization using `Zip`.
 ///
@@ -478,12 +473,12 @@ pub type Ixs = isize;
 /// [`.columns()`][gc], [`.columns_mut()`][gcm],
 /// [`.lanes(axis)`][l], [`.lanes_mut(axis)`][lm].
 ///
-/// [gr]: #method.rows
-/// [grm]: #method.rows_mut
-/// [gc]: #method.columns
-/// [gcm]: #method.columns_mut
-/// [l]: #method.lanes
-/// [lm]: #method.lanes_mut
+/// [gr]: Self::rows
+/// [grm]: Self::rows_mut
+/// [gc]: Self::columns
+/// [gcm]: Self::columns_mut
+/// [l]: Self::lanes
+/// [lm]: Self::lanes_mut
 ///
 /// Yes, for 2D arrays `.rows()` and `.outer_iter()` have about the same
 /// effect:
@@ -497,7 +492,7 @@ pub type Ixs = isize;
 /// the array. Slicing methods include [`.slice()`], [`.slice_mut()`],
 /// [`.slice_move()`], and [`.slice_collapse()`].
 ///
-/// The slicing argument can be passed using the macro [`s![]`](macro.s!.html),
+/// The slicing argument can be passed using the macro [`s![]`](s!),
 /// which will be used in all examples. (The explicit form is an instance of
 /// [`SliceInfo`] or another type which implements [`SliceArg`]; see their docs
 /// for more information.)
@@ -509,11 +504,10 @@ pub type Ixs = isize;
 /// [`.slice_collapse()`] panics on `NewAxis` elements and behaves like
 /// [`.collapse_axis()`] by preserving the number of dimensions.
 ///
-/// [`.slice()`]: #method.slice
-/// [`.slice_mut()`]: #method.slice_mut
-/// [`.slice_move()`]: #method.slice_move
-/// [`.slice_collapse()`]: #method.slice_collapse
-/// [`NewAxis`]: struct.NewAxis.html
+/// [`.slice()`]: Self::slice
+/// [`.slice_mut()`]: Self::slice_mut
+/// [`.slice_move()`]: Self::slice_move
+/// [`.slice_collapse()`]: Self::slice_collapse
 ///
 /// When slicing arrays with generic dimensionality, creating an instance of
 /// [`SliceInfo`] to pass to the multi-axis slicing methods like [`.slice()`]
@@ -522,18 +516,18 @@ pub type Ixs = isize;
 /// or to create a view and then slice individual axes of the view using
 /// methods such as [`.slice_axis_inplace()`] and [`.collapse_axis()`].
 ///
-/// [`.slice_each_axis()`]: #method.slice_each_axis
-/// [`.slice_each_axis_mut()`]: #method.slice_each_axis_mut
-/// [`.slice_each_axis_inplace()`]: #method.slice_each_axis_inplace
-/// [`.slice_axis_inplace()`]: #method.slice_axis_inplace
-/// [`.collapse_axis()`]: #method.collapse_axis
+/// [`.slice_each_axis()`]: Self::slice_each_axis
+/// [`.slice_each_axis_mut()`]: Self::slice_each_axis_mut
+/// [`.slice_each_axis_inplace()`]: Self::slice_each_axis_inplace
+/// [`.slice_axis_inplace()`]: Self::slice_axis_inplace
+/// [`.collapse_axis()`]: Self::collapse_axis
 ///
 /// It's possible to take multiple simultaneous *mutable* slices with
 /// [`.multi_slice_mut()`] or (for [`ArrayViewMut`] only)
 /// [`.multi_slice_move()`].
 ///
-/// [`.multi_slice_mut()`]: #method.multi_slice_mut
-/// [`.multi_slice_move()`]: type.ArrayViewMut.html#method.multi_slice_move
+/// [`.multi_slice_mut()`]: Self::multi_slice_mut
+/// [`.multi_slice_move()`]: ArrayViewMut#method.multi_slice_move
 ///
 /// ```
 /// use ndarray::{arr2, arr3, s, ArrayBase, DataMut, Dimension, NewAxis, Slice};
@@ -631,16 +625,16 @@ pub type Ixs = isize;
 /// Methods for selecting an individual subview take two arguments: `axis` and
 /// `index`.
 ///
-/// [`.axis_iter()`]: #method.axis_iter
-/// [`.axis_iter_mut()`]: #method.axis_iter_mut
-/// [`.fold_axis()`]: #method.fold_axis
-/// [`.index_axis()`]: #method.index_axis
-/// [`.index_axis_inplace()`]: #method.index_axis_inplace
-/// [`.index_axis_mut()`]: #method.index_axis_mut
-/// [`.index_axis_move()`]: #method.index_axis_move
-/// [`.collapse_axis()`]: #method.collapse_axis
-/// [`.outer_iter()`]: #method.outer_iter
-/// [`.outer_iter_mut()`]: #method.outer_iter_mut
+/// [`.axis_iter()`]: Self::axis_iter
+/// [`.axis_iter_mut()`]: Self::axis_iter_mut
+/// [`.fold_axis()`]: Self::fold_axis
+/// [`.index_axis()`]: Self::index_axis
+/// [`.index_axis_inplace()`]: Self::index_axis_inplace
+/// [`.index_axis_mut()`]: Self::index_axis_mut
+/// [`.index_axis_move()`]: Self::index_axis_move
+/// [`.collapse_axis()`]: Self::collapse_axis
+/// [`.outer_iter()`]: Self::outer_iter
+/// [`.outer_iter_mut()`]: Self::outer_iter_mut
 ///
 /// ```
 ///
@@ -723,7 +717,7 @@ pub type Ixs = isize;
 ///
 /// ### Binary Operators with Array and Scalar
 ///
-/// The trait [`ScalarOperand`](trait.ScalarOperand.html) marks types that can be used in arithmetic
+/// The trait [`ScalarOperand`] marks types that can be used in arithmetic
 /// with arrays directly. For a scalar `K` the following combinations of operands
 /// are supported (scalar can be on either the left or right side, but
 /// `ScalarOperand` docs has the detailed condtions).
@@ -746,7 +740,7 @@ pub type Ixs = isize;
 /// Arrays support limited *broadcasting*, where arithmetic operations with
 /// array operands of different sizes can be carried out by repeating the
 /// elements of the smaller dimension array. See
-/// [`.broadcast()`](#method.broadcast) for a more detailed
+/// [`.broadcast()`](Self::broadcast) for a more detailed
 /// description.
 ///
 /// ```
@@ -892,12 +886,12 @@ pub type Ixs = isize;
 /// </td>
 /// <td>
 ///
-/// [`CowArray::from(a)`](type.CowArray.html#impl-From<ArrayBase<OwnedRepr<A>%2C%20D>>)
+/// [`CowArray::from(a)`](CowArray#impl-From<ArrayBase<OwnedRepr<A>%2C%20D>>)
 ///
 /// </td>
 /// <td>
 ///
-/// [`CowArray::from(a.into_owned())`](type.CowArray.html#impl-From<ArrayBase<OwnedRepr<A>%2C%20D>>)
+/// [`CowArray::from(a.into_owned())`](CowArray#impl-From<ArrayBase<OwnedRepr<A>%2C%20D>>)
 ///
 /// </td>
 /// <td>
@@ -907,12 +901,12 @@ pub type Ixs = isize;
 /// </td>
 /// <td>
 ///
-/// [`CowArray::from(a)`](type.CowArray.html#impl-From<ArrayBase<ViewRepr<%26%27a%20A>%2C%20D>>)
+/// [`CowArray::from(a)`](CowArray#impl-From<ArrayBase<ViewRepr<%26%27a%20A>%2C%20D>>)
 ///
 /// </td>
 /// <td>
 ///
-/// [`CowArray::from(a.view())`](type.CowArray.html#impl-From<ArrayBase<ViewRepr<%26%27a%20A>%2C%20D>>)
+/// [`CowArray::from(a.view())`](CowArray#impl-From<ArrayBase<ViewRepr<%26%27a%20A>%2C%20D>>)
 ///
 /// </td>
 /// </tr>
@@ -1041,19 +1035,19 @@ pub type Ixs = isize;
 ///
 /// Input | Output | Methods
 /// ------|--------|--------
-/// `Vec<A>` | `ArrayBase<S: DataOwned, Ix1>` | [`::from_vec()`](#method.from_vec)
-/// `Vec<A>` | `ArrayBase<S: DataOwned, D>` | [`::from_shape_vec()`](#method.from_shape_vec)
-/// `&[A]` | `ArrayView1<A>` | [`::from()`](type.ArrayView.html#method.from)
-/// `&[A]` | `ArrayView<A, D>` | [`::from_shape()`](type.ArrayView.html#method.from_shape)
-/// `&mut [A]` | `ArrayViewMut1<A>` | [`::from()`](type.ArrayViewMut.html#method.from)
-/// `&mut [A]` | `ArrayViewMut<A, D>` | [`::from_shape()`](type.ArrayViewMut.html#method.from_shape)
-/// `&ArrayBase<S, Ix1>` | `Vec<A>` | [`.to_vec()`](#method.to_vec)
-/// `Array<A, D>` | `Vec<A>` | [`.into_raw_vec()`](type.Array.html#method.into_raw_vec)<sup>[1](#into_raw_vec)</sup>
-/// `&ArrayBase<S, D>` | `&[A]` | [`.as_slice()`](#method.as_slice)<sup>[2](#req_contig_std)</sup>, [`.as_slice_memory_order()`](#method.as_slice_memory_order)<sup>[3](#req_contig)</sup>
-/// `&mut ArrayBase<S: DataMut, D>` | `&mut [A]` | [`.as_slice_mut()`](#method.as_slice_mut)<sup>[2](#req_contig_std)</sup>, [`.as_slice_memory_order_mut()`](#method.as_slice_memory_order_mut)<sup>[3](#req_contig)</sup>
-/// `ArrayView<A, D>` | `&[A]` | [`.to_slice()`](type.ArrayView.html#method.to_slice)<sup>[2](#req_contig_std)</sup>
-/// `ArrayViewMut<A, D>` | `&mut [A]` | [`.into_slice()`](type.ArrayViewMut.html#method.into_slice)<sup>[2](#req_contig_std)</sup>
-/// `Array0<A>` | `A` | [`.into_scalar()`](type.Array.html#method.into_scalar)
+/// `Vec<A>` | `ArrayBase<S: DataOwned, Ix1>` | [`::from_vec()`](Self::from_vec)
+/// `Vec<A>` | `ArrayBase<S: DataOwned, D>` | [`::from_shape_vec()`](Self::from_shape_vec)
+/// `&[A]` | `ArrayView1<A>` | [`::from()`](ArrayView#method.from)
+/// `&[A]` | `ArrayView<A, D>` | [`::from_shape()`](ArrayView#method.from_shape)
+/// `&mut [A]` | `ArrayViewMut1<A>` | [`::from()`](ArrayViewMut#method.from)
+/// `&mut [A]` | `ArrayViewMut<A, D>` | [`::from_shape()`](ArrayViewMut#method.from_shape)
+/// `&ArrayBase<S, Ix1>` | `Vec<A>` | [`.to_vec()`](Self::to_vec)
+/// `Array<A, D>` | `Vec<A>` | [`.into_raw_vec()`](Array#method.into_raw_vec)<sup>[1](#into_raw_vec)</sup>
+/// `&ArrayBase<S, D>` | `&[A]` | [`.as_slice()`](Self::as_slice)<sup>[2](#req_contig_std)</sup>, [`.as_slice_memory_order()`](Self::as_slice_memory_order)<sup>[3](#req_contig)</sup>
+/// `&mut ArrayBase<S: DataMut, D>` | `&mut [A]` | [`.as_slice_mut()`](Self::as_slice_mut)<sup>[2](#req_contig_std)</sup>, [`.as_slice_memory_order_mut()`](Self::as_slice_memory_order_mut)<sup>[3](#req_contig)</sup>
+/// `ArrayView<A, D>` | `&[A]` | [`.to_slice()`](ArrayView#method.to_slice)<sup>[2](#req_contig_std)</sup>
+/// `ArrayViewMut<A, D>` | `&mut [A]` | [`.into_slice()`](ArrayViewMut#method.into_slice)<sup>[2](#req_contig_std)</sup>
+/// `Array0<A>` | `A` | [`.into_scalar()`](Array#method.into_scalar)
 ///
 /// <sup><a name="into_raw_vec">1</a></sup>Returns the data in memory order.
 ///
@@ -1066,16 +1060,16 @@ pub type Ixs = isize;
 /// conversions to/from `Vec`s/slices. See
 /// [below](#constructor-methods-for-owned-arrays) for more constructors.
 ///
-/// [ArrayView::reborrow()]: type.ArrayView.html#method.reborrow
-/// [ArrayViewMut::reborrow()]: type.ArrayViewMut.html#method.reborrow
-/// [.into_dimensionality()]: #method.into_dimensionality
-/// [.into_dyn()]: #method.into_dyn
-/// [.into_owned()]: #method.into_owned
-/// [.into_shared()]: #method.into_shared
-/// [.to_owned()]: #method.to_owned
-/// [.map()]: #method.map
-/// [.view()]: #method.view
-/// [.view_mut()]: #method.view_mut
+/// [ArrayView::reborrow()]: ArrayView#method.reborrow
+/// [ArrayViewMut::reborrow()]: ArrayViewMut#method.reborrow
+/// [.into_dimensionality()]: Self::into_dimensionality
+/// [.into_dyn()]: Self::into_dyn
+/// [.into_owned()]: Self::into_owned
+/// [.into_shared()]: Self::into_shared
+/// [.to_owned()]: Self::to_owned
+/// [.map()]: Self::map
+/// [.view()]: Self::view
+/// [.view_mut()]: Self::view_mut
 ///
 /// ### Conversions from Nested `Vec`s/`Array`s
 ///
@@ -1117,7 +1111,7 @@ pub type Ixs = isize;
 /// If you don't know ahead-of-time the shape of the final array, then the
 /// cleanest solution is generally to append the data to a flat `Vec`, and then
 /// convert it to an `Array` at the end with
-/// [`::from_shape_vec()`](#method.from_shape_vec). You just have to be careful
+/// [`::from_shape_vec()`](Self::from_shape_vec). You just have to be careful
 /// that the layout of the data (the order of the elements in the flat `Vec`)
 /// is correct.
 ///
@@ -1140,10 +1134,9 @@ pub type Ixs = isize;
 ///
 /// If neither of these options works for you, and you really need to convert
 /// nested `Vec`/`Array` instances to an `Array`, the cleanest solution is
-/// generally to use
-/// [`Iterator::flatten()`](https://doc.rust-lang.org/std/iter/trait.Iterator.html#method.flatten)
+/// generally to use [`Iterator::flatten()`]
 /// to get a flat `Vec`, and then convert the `Vec` to an `Array` with
-/// [`::from_shape_vec()`](#method.from_shape_vec), like this:
+/// [`::from_shape_vec()`](Self::from_shape_vec), like this:
 ///
 /// ```rust
 /// use ndarray::{array, Array2, Array3};
@@ -1306,20 +1299,20 @@ where
 /// It can act as both an owner as the data as well as a shared reference (view
 /// like).
 /// Calling a method for mutating elements on `ArcArray`, for example
-/// [`view_mut()`](struct.ArrayBase.html#method.view_mut) or
-/// [`get_mut()`](struct.ArrayBase.html#method.get_mut), will break sharing and
+/// [`view_mut()`](ArrayBase::view_mut) or
+/// [`get_mut()`](ArrayBase::get_mut), will break sharing and
 /// require a clone of the data (if it is not uniquely held).
 ///
 /// `ArcArray` uses atomic reference counting like `Arc`, so it is `Send` and
 /// `Sync` (when allowed by the element type of the array too).
 ///
-/// [**`ArrayBase`**](struct.ArrayBase.html) is used to implement both the owned
+/// **[`ArrayBase`]** is used to implement both the owned
 /// arrays and the views; see its docs for an overview of all array features.
 ///
 /// See also:
 ///
-/// + [Constructor Methods for Owned Arrays](struct.ArrayBase.html#constructor-methods-for-owned-arrays)
-/// + [Methods For All Array Types](struct.ArrayBase.html#methods-for-all-array-types)
+/// + [Constructor Methods for Owned Arrays](ArrayBase#constructor-methods-for-owned-arrays)
+/// + [Methods For All Array Types](ArrayBase#methods-for-all-array-types)
 pub type ArcArray<A, D> = ArrayBase<OwnedArcRepr<A>, D>;
 
 /// An array that owns its data uniquely.
@@ -1330,18 +1323,18 @@ pub type ArcArray<A, D> = ArrayBase<OwnedArcRepr<A>, D>;
 /// The `Array<A, D>` is parameterized by `A` for the element type and `D` for
 /// the dimensionality.
 ///
-/// [**`ArrayBase`**](struct.ArrayBase.html) is used to implement both the owned
+/// **[`ArrayBase`]** is used to implement both the owned
 /// arrays and the views; see its docs for an overview of all array features.
 ///
 /// See also:
 ///
-/// + [Constructor Methods for Owned Arrays](struct.ArrayBase.html#constructor-methods-for-owned-arrays)
-/// + [Methods For All Array Types](struct.ArrayBase.html#methods-for-all-array-types)
+/// + [Constructor Methods for Owned Arrays](ArrayBase#constructor-methods-for-owned-arrays)
+/// + [Methods For All Array Types](ArrayBase#methods-for-all-array-types)
 /// + Dimensionality-specific type alises
-/// [`Array1`](type.Array1.html),
-/// [`Array2`](type.Array2.html),
-/// [`Array3`](type.Array3.html), ...,
-/// [`ArrayD`](type.ArrayD.html),
+/// [`Array1`],
+/// [`Array2`],
+/// [`Array3`], ...,
+/// [`ArrayD`],
 /// and so on.
 pub type Array<A, D> = ArrayBase<OwnedRepr<A>, D>;
 
@@ -1350,20 +1343,17 @@ pub type Array<A, D> = ArrayBase<OwnedRepr<A>, D>;
 /// An `CowArray` represents either a uniquely owned array or a view of an
 /// array. The `'a` corresponds to the lifetime of the view variant.
 ///
-/// This type is analogous to
-/// [`std::borrow::Cow`](https://doc.rust-lang.org/std/borrow/enum.Cow.html).
+/// This type is analogous to [`std::borrow::Cow`].
 /// If a `CowArray` instance is the immutable view variant, then calling a
 /// method for mutating elements in the array will cause it to be converted
 /// into the owned variant (by cloning all the elements) before the
 /// modification is performed.
 ///
-/// Array views have all the methods of an array (see [`ArrayBase`][ab]).
+/// Array views have all the methods of an array (see [`ArrayBase`]).
 ///
-/// See also [`ArcArray`](type.ArcArray.html), which also provides
+/// See also [`ArcArray`], which also provides
 /// copy-on-write behavior but has a reference-counted pointer to the data
 /// instead of either a view or a uniquely owned copy.
-///
-/// [ab]: struct.ArrayBase.html
 pub type CowArray<'a, A, D> = ArrayBase<CowRepr<'a, A>, D>;
 
 /// A read-only array view.
@@ -1374,11 +1364,9 @@ pub type CowArray<'a, A, D> = ArrayBase<CowRepr<'a, A>, D>;
 /// The `ArrayView<'a, A, D>` is parameterized by `'a` for the scope of the
 /// borrow, `A` for the element type and `D` for the dimensionality.
 ///
-/// Array views have all the methods of an array (see [`ArrayBase`][ab]).
-///
-/// See also [`ArrayViewMut`](type.ArrayViewMut.html).
+/// Array views have all the methods of an array (see [`ArrayBase`]).
 ///
-/// [ab]: struct.ArrayBase.html
+/// See also [`ArrayViewMut`].
 pub type ArrayView<'a, A, D> = ArrayBase<ViewRepr<&'a A>, D>;
 
 /// A read-write array view.
@@ -1389,11 +1377,9 @@ pub type ArrayView<'a, A, D> = ArrayBase<ViewRepr<&'a A>, D>;
 /// The `ArrayViewMut<'a, A, D>` is parameterized by `'a` for the scope of the
 /// borrow, `A` for the element type and `D` for the dimensionality.
 ///
-/// Array views have all the methods of an array (see [`ArrayBase`][ab]).
+/// Array views have all the methods of an array (see [`ArrayBase`]).
 ///
-/// See also [`ArrayView`](type.ArrayView.html).
-///
-/// [ab]: struct.ArrayBase.html
+/// See also [`ArrayView`].
 pub type ArrayViewMut<'a, A, D> = ArrayBase<ViewRepr<&'a mut A>, D>;
 
 /// A read-only array view without a lifetime.
@@ -1405,15 +1391,13 @@ pub type ArrayViewMut<'a, A, D> = ArrayBase<ViewRepr<&'a mut A>, D>;
 /// T` and `&T`, but `RawArrayView` has additional requirements that `*const T`
 /// does not, such as non-nullness.
 ///
-/// [`ArrayView`]: type.ArrayView.html
-///
 /// The `RawArrayView<A, D>` is parameterized by `A` for the element type and
 /// `D` for the dimensionality.
 ///
 /// Raw array views have all the methods of an array (see
-/// [`ArrayBase`](struct.ArrayBase.html)).
+/// [`ArrayBase`]).
 ///
-/// See also [`RawArrayViewMut`](type.RawArrayViewMut.html).
+/// See also [`RawArrayViewMut`].
 ///
 /// # Warning
 ///
@@ -1430,15 +1414,13 @@ pub type RawArrayView<A, D> = ArrayBase<RawViewRepr<*const A>, D>;
 /// relationship between `*mut T` and `&mut T`, but `RawArrayViewMut` has
 /// additional requirements that `*mut T` does not, such as non-nullness.
 ///
-/// [`ArrayViewMut`]: type.ArrayViewMut.html
-///
 /// The `RawArrayViewMut<A, D>` is parameterized by `A` for the element type
 /// and `D` for the dimensionality.
 ///
 /// Raw array views have all the methods of an array (see
-/// [`ArrayBase`](struct.ArrayBase.html)).
+/// [`ArrayBase`]).
 ///
-/// See also [`RawArrayView`](type.RawArrayView.html).
+/// See also [`RawArrayView`].
 ///
 /// # Warning
 ///
@@ -1451,7 +1433,7 @@ pub use data_repr::OwnedRepr;
 /// ArcArray's representation.
 ///
 /// *Don’t use this type directly—use the type alias
-/// [`ArcArray`](type.ArcArray.html) for the array type!*
+/// [`ArcArray`] for the array type!*
 #[derive(Debug)]
 pub struct OwnedArcRepr<A>(Arc<OwnedRepr<A>>);
 
@@ -1464,8 +1446,7 @@ impl<A> Clone for OwnedArcRepr<A> {
 /// Array pointer’s representation.
 ///
 /// *Don’t use this type directly—use the type aliases
-/// [`RawArrayView`](type.RawArrayView.html) /
-/// [`RawArrayViewMut`](type.RawArrayViewMut.html) for the array type!*
+/// [`RawArrayView`] / [`RawArrayViewMut`] for the array type!*
 #[derive(Copy, Clone)]
 // This is just a marker type, to carry the mutability and element type.
 pub struct RawViewRepr<A> {
@@ -1482,8 +1463,7 @@ impl<A> RawViewRepr<A> {
 /// Array view’s representation.
 ///
 /// *Don’t use this type directly—use the type aliases
-/// [`ArrayView`](type.ArrayView.html)
-/// / [`ArrayViewMut`](type.ArrayViewMut.html) for the array type!*
+/// [`ArrayView`] / [`ArrayViewMut`] for the array type!*
 #[derive(Copy, Clone)]
 // This is just a marker type, to carry the lifetime parameter.
 pub struct ViewRepr<A> {
@@ -1500,7 +1480,7 @@ impl<A> ViewRepr<A> {
 /// CowArray's representation.
 ///
 /// *Don't use this type directly—use the type alias
-/// [`CowArray`](type.CowArray.html) for the array type!*
+/// [`CowArray`] for the array type!*
 pub enum CowRepr<'a, A> {
     /// Borrowed data.
     View(ViewRepr<&'a A>),
diff --git a/src/linalg/impl_linalg.rs b/src/linalg/impl_linalg.rs
index 61b91eae..d851f2f0 100644
--- a/src/linalg/impl_linalg.rs
+++ b/src/linalg/impl_linalg.rs
@@ -8,6 +8,8 @@
 
 use crate::imp_prelude::*;
 use crate::numeric_util;
+#[cfg(feature = "blas")]
+use crate::dimension::offset_from_low_addr_ptr_to_logical_ptr;
 
 use crate::{LinalgScalar, Zip};
 
@@ -649,6 +651,12 @@ unsafe fn general_mat_vec_mul_impl<A, S1, S2>(
                             }
                         };
 
+                        // Low addr in memory pointers required for x, y
+                        let x_offset = offset_from_low_addr_ptr_to_logical_ptr(&x.dim, &x.strides);
+                        let x_ptr = x.ptr.as_ptr().sub(x_offset);
+                        let y_offset = offset_from_low_addr_ptr_to_logical_ptr(&y.dim, &y.strides);
+                        let y_ptr = y.ptr.as_ptr().sub(y_offset);
+
                         let x_stride = x.strides()[0] as blas_index;
                         let y_stride = y.strides()[0] as blas_index;
 
@@ -660,10 +668,10 @@ unsafe fn general_mat_vec_mul_impl<A, S1, S2>(
                             cast_as(&alpha),            // alpha
                             a.ptr.as_ptr() as *const _, // a
                             a_stride,                   // lda
-                            x.ptr.as_ptr() as *const _, // x
+                            x_ptr as *const _,          // x
                             x_stride,
-                            cast_as(&beta),           // beta
-                            y.ptr.as_ptr() as *mut _, // x
+                            cast_as(&beta),             // beta
+                            y_ptr as *mut _,            // y
                             y_stride,
                         );
                         return;
@@ -719,7 +727,10 @@ where
         return false;
     }
     let stride = a.strides()[0];
-    if stride > blas_index::max_value() as isize || stride < blas_index::min_value() as isize {
+    if stride == 0
+        || stride > blas_index::max_value() as isize
+        || stride < blas_index::min_value() as isize
+    {
         return false;
     }
     true
diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index da170c9e..afcd7e89 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -12,7 +12,6 @@ use num_traits::{self, FromPrimitive, Zero};
 use std::ops::{Add, Div, Mul};
 
 use crate::imp_prelude::*;
-use crate::itertools::enumerate;
 use crate::numeric_util;
 
 /// # Numerical Methods for Arrays
@@ -246,22 +245,16 @@ where
         A: Clone + Zero + Add<Output = A>,
         D: RemoveAxis,
     {
-        let n = self.len_of(axis);
-        let mut res = Array::zeros(self.raw_dim().remove_axis(axis));
-        let stride = self.strides()[axis.index()];
-        if self.ndim() == 2 && stride == 1 {
-            // contiguous along the axis we are summing
-            let ax = axis.index();
-            for (i, elt) in enumerate(&mut res) {
-                *elt = self.index_axis(Axis(1 - ax), i).sum();
-            }
+        let min_stride_axis = self.dim.min_stride_axis(&self.strides);
+        if axis == min_stride_axis {
+            crate::Zip::from(self.lanes(axis)).map_collect(|lane| lane.sum())
         } else {
-            for i in 0..n {
-                let view = self.index_axis(axis, i);
-                res = res + &view;
+            let mut res = Array::zeros(self.raw_dim().remove_axis(axis));
+            for subview in self.axis_iter(axis) {
+                res = res + &subview;
             }
+            res
         }
-        res
     }
 
     /// Return mean along `axis`.
diff --git a/src/parallel/impl_par_methods.rs b/src/parallel/impl_par_methods.rs
index 33d39334..381b86b5 100644
--- a/src/parallel/impl_par_methods.rs
+++ b/src/parallel/impl_par_methods.rs
@@ -177,6 +177,56 @@ macro_rules! zip_impl {
                 self.par_map_assign_into(into, f)
             }
 
+            /// Parallel version of `fold`.
+            ///
+            /// Splits the producer in multiple tasks which each accumulate a single value
+            /// using the `fold` closure. Those tasks are executed in parallel and their results
+            /// are then combined to a single value using the `reduce` closure.
+            ///
+            /// The `identity` closure provides the initial values for each of the tasks and
+            /// for the final reduction.
+            ///
+            /// This is a shorthand for calling `self.into_par_iter().fold(...).reduce(...)`.
+            ///
+            /// Note that it is often more efficient to parallelize not per-element but rather
+            /// based on larger chunks of an array like generalized rows and operating on each chunk
+            /// using a sequential variant of the accumulation.
+            /// For example, sum each row sequentially and in parallel, taking advatange of locality
+            /// and vectorization within each task, and then reduce their sums to the sum of the matrix.
+            ///
+            /// Also note that the splitting of the producer into multiple tasks is _not_ deterministic
+            /// which needs to be considered when the accuracy of such an operation is analyzed.
+            ///
+            /// ## Examples
+            ///
+            /// ```rust
+            /// use ndarray::{Array, Zip};
+            ///
+            /// let a = Array::<usize, _>::ones((128, 1024));
+            /// let b = Array::<usize, _>::ones(128);
+            ///
+            /// let weighted_sum = Zip::from(a.rows()).and(&b).par_fold(
+            ///     || 0,
+            ///     |sum, row, factor| sum + row.sum() * factor,
+            ///     |sum, other_sum| sum + other_sum,
+            /// );
+            ///
+            /// assert_eq!(weighted_sum, a.len());
+            /// ```
+            pub fn par_fold<ID, F, R, T>(self, identity: ID, fold: F, reduce: R) -> T
+            where
+                ID: Fn() -> T + Send + Sync + Clone,
+                F: Fn(T, $($p::Item),*) -> T + Send + Sync,
+                R: Fn(T, T) -> T + Send + Sync,
+                T: Send
+            {
+                self.into_par_iter()
+                    .fold(identity.clone(), move |accumulator, ($($p,)*)| {
+                        fold(accumulator, $($p),*)
+                    })
+                    .reduce(identity, reduce)
+            }
+
             );
         }
         )+
diff --git a/src/parallel/mod.rs b/src/parallel/mod.rs
index dd3b8934..48e97be4 100644
--- a/src/parallel/mod.rs
+++ b/src/parallel/mod.rs
@@ -32,6 +32,10 @@
 //! “unindexed”. Use ndarray’s [Zip] for lock step parallel iteration of
 //! multiple arrays or producers at a time.
 //!
+//! For the unindexed parallel iterators, an inherent method [`with_min_len`](Parallel::with_min_len)
+//! is provided to limit the number of elements each parallel task processes in way that is
+//! similar to Rayon's [`IndexedParallelIterator::with_min_len`](rayon::prelude::IndexedParallelIterator::with_min_len).
+//!
 //! # Examples
 //!
 //! ## Arrays and array views
@@ -42,16 +46,14 @@
 //! use ndarray::Array2;
 //! use ndarray::parallel::prelude::*;
 //!
-//! fn main() {
-//!     let mut a = Array2::<f64>::zeros((128, 128));
+//! let mut a = Array2::<f64>::zeros((128, 128));
 //!
-//!     // Parallel versions of regular array methods
-//!     a.par_map_inplace(|x| *x = x.exp());
-//!     a.par_mapv_inplace(f64::exp);
+//! // Parallel versions of regular array methods
+//! a.par_map_inplace(|x| *x = x.exp());
+//! a.par_mapv_inplace(f64::exp);
 //!
-//!     // You can also use the parallel iterator directly
-//!     a.par_iter_mut().for_each(|x| *x = x.exp());
-//! }
+//! // You can also use the parallel iterator directly
+//! a.par_iter_mut().for_each(|x| *x = x.exp());
 //! ```
 //!
 //! ## Axis iterators
@@ -63,16 +65,14 @@
 //! use ndarray::Axis;
 //! use ndarray::parallel::prelude::*;
 //!
-//! fn main() {
-//!     let a = Array::linspace(0., 63., 64).into_shape((4, 16)).unwrap();
-//!     let mut sums = Vec::new();
-//!     a.axis_iter(Axis(0))
-//!      .into_par_iter()
-//!      .map(|row| row.sum())
-//!      .collect_into_vec(&mut sums);
+//! let a = Array::linspace(0., 63., 64).into_shape((4, 16)).unwrap();
+//! let mut sums = Vec::new();
+//! a.axis_iter(Axis(0))
+//!  .into_par_iter()
+//!  .map(|row| row.sum())
+//!  .collect_into_vec(&mut sums);
 //!
-//!     assert_eq!(sums, [120., 376., 632., 888.]);
-//! }
+//! assert_eq!(sums, [120., 376., 632., 888.]);
 //! ```
 //!
 //! ## Axis chunks iterators
@@ -84,16 +84,14 @@
 //! use ndarray::Axis;
 //! use ndarray::parallel::prelude::*;
 //!
-//! fn main() {
-//!     let a = Array::linspace(0., 63., 64).into_shape((4, 16)).unwrap();
-//!     let mut shapes = Vec::new();
-//!     a.axis_chunks_iter(Axis(0), 3)
-//!      .into_par_iter()
-//!      .map(|chunk| chunk.shape().to_owned())
-//!      .collect_into_vec(&mut shapes);
+//! let a = Array::linspace(0., 63., 64).into_shape((4, 16)).unwrap();
+//! let mut shapes = Vec::new();
+//! a.axis_chunks_iter(Axis(0), 3)
+//!  .into_par_iter()
+//!  .map(|chunk| chunk.shape().to_owned())
+//!  .collect_into_vec(&mut shapes);
 //!
-//!     assert_eq!(shapes, [vec![3, 16], vec![1, 16]]);
-//! }
+//! assert_eq!(shapes, [vec![3, 16], vec![1, 16]]);
 //! ```
 //!
 //! ## Zip
@@ -106,19 +104,17 @@
 //!
 //! type Array3f64 = Array3<f64>;
 //!
-//! fn main() {
-//!     const N: usize = 128;
-//!     let a = Array3f64::from_elem((N, N, N), 1.);
-//!     let b = Array3f64::from_elem(a.dim(), 2.);
-//!     let mut c = Array3f64::zeros(a.dim());
-//!
-//!     Zip::from(&mut c)
-//!         .and(&a)
-//!         .and(&b)
-//!         .par_for_each(|c, &a, &b| {
-//!             *c += a - b;
-//!         });
-//! }
+//! const N: usize = 128;
+//! let a = Array3f64::from_elem((N, N, N), 1.);
+//! let b = Array3f64::from_elem(a.dim(), 2.);
+//! let mut c = Array3f64::zeros(a.dim());
+//!
+//! Zip::from(&mut c)
+//!     .and(&a)
+//!     .and(&b)
+//!     .par_for_each(|c, &a, &b| {
+//!         *c += a - b;
+//!     });
 //! ```
 
 #[allow(unused_imports)] // used by rustdoc links
diff --git a/src/parallel/par.rs b/src/parallel/par.rs
index d9d592af..cc905b5c 100644
--- a/src/parallel/par.rs
+++ b/src/parallel/par.rs
@@ -21,11 +21,14 @@ use crate::split_at::SplitPreference;
 #[derive(Copy, Clone, Debug)]
 pub struct Parallel<I> {
     iter: I,
+    min_len: usize,
 }
 
+const DEFAULT_MIN_LEN: usize = 1;
+
 /// Parallel producer wrapper.
 #[derive(Copy, Clone, Debug)]
-struct ParallelProducer<I>(I);
+struct ParallelProducer<I>(I, usize);
 
 macro_rules! par_iter_wrapper {
     // thread_bounds are either Sync or Send + Sync
@@ -40,6 +43,7 @@ macro_rules! par_iter_wrapper {
         fn into_par_iter(self) -> Self::Iter {
             Parallel {
                 iter: self,
+                min_len: DEFAULT_MIN_LEN,
             }
         }
     }
@@ -67,7 +71,7 @@ macro_rules! par_iter_wrapper {
         fn with_producer<Cb>(self, callback: Cb) -> Cb::Output
             where Cb: ProducerCallback<Self::Item>
         {
-            callback.callback(ParallelProducer(self.iter))
+            callback.callback(ParallelProducer(self.iter, self.min_len))
         }
 
         fn len(&self) -> usize {
@@ -106,7 +110,7 @@ macro_rules! par_iter_wrapper {
 
         fn split_at(self, i: usize) -> (Self, Self) {
             let (a, b) = self.0.split_at(i);
-            (ParallelProducer(a), ParallelProducer(b))
+            (ParallelProducer(a, self.1), ParallelProducer(b, self.1))
         }
     }
 
@@ -131,11 +135,11 @@ macro_rules! par_iter_view_wrapper {
         fn into_par_iter(self) -> Self::Iter {
             Parallel {
                 iter: self,
+                min_len: DEFAULT_MIN_LEN,
             }
         }
     }
 
-
     impl<'a, A, D> ParallelIterator for Parallel<$view_name<'a, A, D>>
         where D: Dimension,
               A: $($thread_bounds)*,
@@ -144,7 +148,7 @@ macro_rules! par_iter_view_wrapper {
         fn drive_unindexed<C>(self, consumer: C) -> C::Result
             where C: UnindexedConsumer<Self::Item>
         {
-            bridge_unindexed(ParallelProducer(self.iter), consumer)
+            bridge_unindexed(ParallelProducer(self.iter, self.min_len), consumer)
         }
 
         fn opt_len(&self) -> Option<usize> {
@@ -152,20 +156,39 @@ macro_rules! par_iter_view_wrapper {
         }
     }
 
+    impl<'a, A, D> Parallel<$view_name<'a, A, D>>
+        where D: Dimension,
+              A: $($thread_bounds)*,
+    {
+        /// Sets the minimum number of elements desired to process in each job. This will not be
+        /// split any smaller than this length, but of course a producer could already be smaller
+        /// to begin with.
+        ///
+        /// ***Panics*** if `min_len` is zero.
+        pub fn with_min_len(self, min_len: usize) -> Self {
+            assert_ne!(min_len, 0, "Minimum number of elements must at least be one to avoid splitting off empty tasks.");
+
+            Self {
+                min_len,
+                ..self
+            }
+        }
+    }
+
     impl<'a, A, D> UnindexedProducer for ParallelProducer<$view_name<'a, A, D>>
         where D: Dimension,
               A: $($thread_bounds)*,
     {
         type Item = <$view_name<'a, A, D> as IntoIterator>::Item;
         fn split(self) -> (Self, Option<Self>) {
-            if self.0.len() <= 1 {
+            if self.0.len() <= self.1 {
                 return (self, None)
             }
             let array = self.0;
             let max_axis = array.max_stride_axis();
             let mid = array.len_of(max_axis) / 2;
             let (a, b) = array.split_at(max_axis, mid);
-            (ParallelProducer(a), Some(ParallelProducer(b)))
+            (ParallelProducer(a, self.1), Some(ParallelProducer(b, self.1)))
         }
 
         fn fold_with<F>(self, folder: F) -> F
@@ -217,6 +240,7 @@ macro_rules! zip_impl {
             fn into_par_iter(self) -> Self::Iter {
                 Parallel {
                     iter: self,
+                    min_len: DEFAULT_MIN_LEN,
                 }
             }
         }
@@ -233,7 +257,7 @@ macro_rules! zip_impl {
             fn drive_unindexed<Cons>(self, consumer: Cons) -> Cons::Result
                 where Cons: UnindexedConsumer<Self::Item>
             {
-                bridge_unindexed(ParallelProducer(self.iter), consumer)
+                bridge_unindexed(ParallelProducer(self.iter, self.min_len), consumer)
             }
 
             fn opt_len(&self) -> Option<usize> {
@@ -251,11 +275,11 @@ macro_rules! zip_impl {
             type Item = ($($p::Item ,)*);
 
             fn split(self) -> (Self, Option<Self>) {
-                if !self.0.can_split() {
+                if self.0.size() <= self.1 {
                     return (self, None)
                 }
                 let (a, b) = self.0.split();
-                (ParallelProducer(a), Some(ParallelProducer(b)))
+                (ParallelProducer(a, self.1), Some(ParallelProducer(b, self.1)))
             }
 
             fn fold_with<Fold>(self, folder: Fold) -> Fold
@@ -284,6 +308,25 @@ zip_impl! {
     [P1 P2 P3 P4 P5 P6],
 }
 
+impl<D, Parts> Parallel<Zip<Parts, D>>
+where
+    D: Dimension,
+{
+    /// Sets the minimum number of elements desired to process in each job. This will not be
+    /// split any smaller than this length, but of course a producer could already be smaller
+    /// to begin with.
+    ///
+    /// ***Panics*** if `min_len` is zero.
+    pub fn with_min_len(self, min_len: usize) -> Self {
+        assert_ne!(min_len, 0, "Minimum number of elements must at least be one to avoid splitting off empty tasks.");
+
+        Self {
+            min_len,
+            ..self
+        }
+    }
+}
+
 /// A parallel iterator (unindexed) that produces the splits of the array
 /// or producer `P`.
 pub(crate) struct ParallelSplits<P> {
diff --git a/src/parallel/zipmacro.rs b/src/parallel/zipmacro.rs
index 95111c93..28188542 100644
--- a/src/parallel/zipmacro.rs
+++ b/src/parallel/zipmacro.rs
@@ -39,18 +39,16 @@
 ///
 /// type M = Array2<f32>;
 ///
-/// fn main() {
-///     let mut a = M::zeros((16, 16));
-///     let b = M::from_elem(a.dim(), 1.);
-///     let c = M::from_elem(a.dim(), 2.);
+/// let mut a = M::zeros((16, 16));
+/// let b = M::from_elem(a.dim(), 1.);
+/// let c = M::from_elem(a.dim(), 2.);
 ///
-///     // Compute a simple ternary operation:
-///     // elementwise addition of b and c, stored in a
+/// // Compute a simple ternary operation:
+/// // elementwise addition of b and c, stored in a
 ///
-///     par_azip!((a in &mut a, &b in &b, &c in &c) *a = b + c);
+/// par_azip!((a in &mut a, &b in &b, &c in &c) *a = b + c);
 ///
-///     assert_eq!(a, &b + &c);
-/// }
+/// assert_eq!(a, &b + &c);
 /// ```
 macro_rules! par_azip {
     ($($t:tt)*) => {
diff --git a/src/partial.rs b/src/partial.rs
index a8146aff..a835081c 100644
--- a/src/partial.rs
+++ b/src/partial.rs
@@ -38,7 +38,7 @@ impl<T> Partial<T> {
 
     #[cfg(feature = "rayon")]
     pub(crate) fn stub() -> Self {
-        Self { len: 0, ptr: 0 as *mut _ }
+        Self { len: 0, ptr: ptr::null_mut() }
     }
 
     #[cfg(feature = "rayon")]
diff --git a/src/slice.rs b/src/slice.rs
index 113b82d0..dd0d162e 100644
--- a/src/slice.rs
+++ b/src/slice.rs
@@ -20,7 +20,7 @@ use std::ops::{Deref, Range, RangeFrom, RangeFull, RangeInclusive, RangeTo, Rang
 /// from the back of the axis. If `end` is `None`, the slice extends to the end
 /// of the axis.
 ///
-/// See also the [`s![]`](macro.s.html) macro.
+/// See also the [`s![]`](s!) macro.
 ///
 /// ## Examples
 ///
@@ -75,13 +75,13 @@ impl Slice {
 
 /// Token to represent a new axis in a slice description.
 ///
-/// See also the [`s![]`](macro.s!.html) macro.
+/// See also the [`s![]`](s!) macro.
 #[derive(Clone, Copy, Debug)]
 pub struct NewAxis;
 
 /// A slice (range with step), an index, or a new axis token.
 ///
-/// See also the [`s![]`](macro.s!.html) macro for a convenient way to create a
+/// See also the [`s![]`](s!) macro for a convenient way to create a
 /// `SliceInfo<[SliceInfoElem; n], Din, Dout>`.
 ///
 /// ## Examples
@@ -295,6 +295,7 @@ impl From<NewAxis> for SliceInfoElem {
 /// that `D`, `Self::OutDim`, `self.in_dim()`, and `self.out_ndim()` are
 /// consistent with the `&[SliceInfoElem]` returned by `self.as_ref()` and that
 /// `self.as_ref()` always returns the same value when called multiple times.
+#[allow(clippy::missing_safety_doc)] // not implementable downstream
 pub unsafe trait SliceArg<D: Dimension>: AsRef<[SliceInfoElem]> {
     /// Dimensionality of the output array.
     type OutDim: Dimension;
@@ -397,7 +398,7 @@ unsafe impl SliceArg<IxDyn> for [SliceInfoElem] {
 /// `SliceInfo` instance can still be used to slice an array with dimension
 /// `IxDyn` as long as the number of axes matches.
 ///
-/// [`.slice()`]: struct.ArrayBase.html#method.slice
+/// [`.slice()`]: crate::ArrayBase::slice
 #[derive(Debug)]
 pub struct SliceInfo<T, Din: Dimension, Dout: Dimension> {
     in_dim: PhantomData<Din>,
@@ -486,7 +487,7 @@ where
     }
 
     /// Returns the number of dimensions of the input array for
-    /// [`.slice()`](struct.ArrayBase.html#method.slice).
+    /// [`.slice()`](crate::ArrayBase::slice).
     ///
     /// If `Din` is a fixed-size dimension type, then this is equivalent to
     /// `Din::NDIM.unwrap()`. Otherwise, the value is calculated by iterating
@@ -500,7 +501,7 @@ where
     }
 
     /// Returns the number of dimensions after calling
-    /// [`.slice()`](struct.ArrayBase.html#method.slice) (including taking
+    /// [`.slice()`](crate::ArrayBase::slice) (including taking
     /// subviews).
     ///
     /// If `Dout` is a fixed-size dimension type, then this is equivalent to
@@ -703,9 +704,6 @@ impl_slicenextdim!((), NewAxis, Ix0, Ix1);
 /// * *new-axis*: a [`NewAxis`] instance that represents the creation of a new axis.
 ///   (Except for [`.slice_collapse()`], which panics on [`NewAxis`] elements.)
 ///
-/// [`Slice`]: struct.Slice.html
-/// [`NewAxis`]: struct.NewAxis.html
-///
 /// The number of *elem*, not including *new-axis*, must match the
 /// number of axes in the array. *index*, *range*, *slice*, *step*, and
 /// *new-axis* can be expressions. *index* must be of type `isize`, `usize`, or
@@ -722,12 +720,12 @@ impl_slicenextdim!((), NewAxis, Ix0, Ix1);
 /// panic. Without the `NewAxis`, i.e. `s![0..4;2, 6, 1..5]`,
 /// [`.slice_collapse()`] would result in an array of shape `[2, 1, 4]`.
 ///
-/// [`.slice()`]: struct.ArrayBase.html#method.slice
-/// [`.slice_mut()`]: struct.ArrayBase.html#method.slice_mut
-/// [`.slice_move()`]: struct.ArrayBase.html#method.slice_move
-/// [`.slice_collapse()`]: struct.ArrayBase.html#method.slice_collapse
+/// [`.slice()`]: crate::ArrayBase::slice
+/// [`.slice_mut()`]: crate::ArrayBase::slice_mut
+/// [`.slice_move()`]: crate::ArrayBase::slice_move
+/// [`.slice_collapse()`]: crate::ArrayBase::slice_collapse
 ///
-/// See also [*Slicing*](struct.ArrayBase.html#slicing).
+/// See also [*Slicing*](crate::ArrayBase#slicing).
 ///
 /// # Example
 ///
diff --git a/src/stacking.rs b/src/stacking.rs
index fb05c196..5a47589d 100644
--- a/src/stacking.rs
+++ b/src/stacking.rs
@@ -178,11 +178,9 @@ where
 
 /// Stack arrays along the new axis.
 ///
-/// Uses the [`stack`][1] function, calling `ArrayView::from(&a)` on each
+/// Uses the [`stack()`] function, calling `ArrayView::from(&a)` on each
 /// argument `a`.
 ///
-/// [1]: fn.stack.html
-///
 /// ***Panics*** if the `stack` function would return an error.
 ///
 /// ```
@@ -192,31 +190,46 @@ where
 ///
 /// # fn main() {
 ///
-/// let a = arr2(&[[2., 2.],
-///                [3., 3.]]);
-/// assert!(
-///     stack![Axis(0), a, a]
-///     == arr3(&[[[2., 2.],
-///                [3., 3.]],
-///               [[2., 2.],
-///                [3., 3.]]])
+/// let a = arr2(&[[1., 2.],
+///                [3., 4.]]);
+/// assert_eq!(
+///     stack![Axis(0), a, a],
+///     arr3(&[[[1., 2.],
+///             [3., 4.]],
+///            [[1., 2.],
+///             [3., 4.]]]),
+/// );
+/// assert_eq!(
+///     stack![Axis(1), a, a,],
+///     arr3(&[[[1., 2.],
+///             [1., 2.]],
+///            [[3., 4.],
+///             [3., 4.]]]),
+/// );
+/// assert_eq!(
+///     stack![Axis(2), a, a],
+///     arr3(&[[[1., 1.],
+///             [2., 2.]],
+///            [[3., 3.],
+///             [4., 4.]]]),
 /// );
 /// # }
 /// ```
 #[macro_export]
 macro_rules! stack {
+    ($axis:expr, $( $array:expr ),+ ,) => {
+        $crate::stack!($axis, $($array),+)
+    };
     ($axis:expr, $( $array:expr ),+ ) => {
         $crate::stack($axis, &[ $($crate::ArrayView::from(&$array) ),* ]).unwrap()
-    }
+    };
 }
 
 /// Concatenate arrays along the given axis.
 ///
-/// Uses the [`concatenate`][1] function, calling `ArrayView::from(&a)` on each
+/// Uses the [`concatenate()`] function, calling `ArrayView::from(&a)` on each
 /// argument `a`.
 ///
-/// [1]: fn.concatenate.html
-///
 /// ***Panics*** if the `concatenate` function would return an error.
 ///
 /// ```
@@ -226,31 +239,37 @@ macro_rules! stack {
 ///
 /// # fn main() {
 ///
-/// let a = arr2(&[[2., 2.],
-///                [3., 3.]]);
-/// assert!(
-///     concatenate![Axis(0), a, a]
-///     == arr2(&[[2., 2.],
-///               [3., 3.],
-///               [2., 2.],
-///               [3., 3.]])
+/// let a = arr2(&[[1., 2.],
+///                [3., 4.]]);
+/// assert_eq!(
+///     concatenate![Axis(0), a, a],
+///     arr2(&[[1., 2.],
+///            [3., 4.],
+///            [1., 2.],
+///            [3., 4.]]),
+/// );
+/// assert_eq!(
+///     concatenate![Axis(1), a, a,],
+///     arr2(&[[1., 2., 1., 2.],
+///            [3., 4., 3., 4.]]),
 /// );
 /// # }
 /// ```
 #[macro_export]
 macro_rules! concatenate {
+    ($axis:expr, $( $array:expr ),+ ,) => {
+        $crate::concatenate!($axis, $($array),+)
+    };
     ($axis:expr, $( $array:expr ),+ ) => {
         $crate::concatenate($axis, &[ $($crate::ArrayView::from(&$array) ),* ]).unwrap()
-    }
+    };
 }
 
 /// Stack arrays along the new axis.
 ///
-/// Uses the [`stack_new_axis`][1] function, calling `ArrayView::from(&a)` on each
+/// Uses the [`stack_new_axis()`] function, calling `ArrayView::from(&a)` on each
 /// argument `a`.
 ///
-/// [1]: fn.stack_new_axis.html
-///
 /// ***Panics*** if the `stack` function would return an error.
 ///
 /// ```
diff --git a/src/zip/mod.rs b/src/zip/mod.rs
index d611fabb..a24adc2c 100644
--- a/src/zip/mod.rs
+++ b/src/zip/mod.rs
@@ -37,9 +37,7 @@ macro_rules! fold_while {
 
 /// Broadcast an array so that it acts like a larger size and/or shape array.
 ///
-/// See [broadcasting][1] for more information.
-///
-/// [1]: struct.ArrayBase.html#broadcasting
+/// See [broadcasting](ArrayBase#broadcasting) for more information.
 trait Broadcast<E>
 where
     E: IntoDimension,
@@ -120,7 +118,7 @@ trait ZippableTuple: Sized {
 /// a time).
 ///
 /// In general, the zip uses a tuple of producers
-/// ([`NdProducer`](trait.NdProducer.html) trait) that all have to be of the
+/// ([`NdProducer`] trait) that all have to be of the
 /// same shape. The NdProducer implementation defines what its item type is
 /// (for example if it's a shared reference, mutable reference or an array
 /// view etc).
@@ -135,11 +133,9 @@ trait ZippableTuple: Sized {
 /// `fold_while`. The zip object can be split, which allows parallelization.
 /// A read-only zip object (no mutable producers) can be cloned.
 ///
-/// See also the [`azip!()` macro][az] which offers a convenient shorthand
+/// See also the [`azip!()`] which offers a convenient shorthand
 /// to common ways to use `Zip`.
 ///
-/// [az]: macro.azip.html
-///
 /// ```
 /// use ndarray::Zip;
 /// use ndarray::Array2;
@@ -190,6 +186,7 @@ trait ZippableTuple: Sized {
 /// assert_eq!(totals, totals2);
 /// ```
 #[derive(Debug, Clone)]
+#[must_use = "zipping producers is lazy and does nothing unless consumed"]
 pub struct Zip<Parts, D> {
     parts: Parts,
     dimension: D,
@@ -252,7 +249,7 @@ where
     P: NdProducer<Dim = D>,
 {
     ndassert!(
-        part.equal_dim(&dimension),
+        part.equal_dim(dimension),
         "Zip: Producer dimension mismatch, expected: {:?}, got: {:?}",
         dimension,
         part.raw_dim()
diff --git a/src/zip/ndproducer.rs b/src/zip/ndproducer.rs
index 7c201741..619fadcc 100644
--- a/src/zip/ndproducer.rs
+++ b/src/zip/ndproducer.rs
@@ -36,8 +36,8 @@ where
 /// for example an array view, mutable array view or an iterator
 /// that yields chunks.
 ///
-/// Producers are used as a arguments to [`Zip`](struct.Zip.html) and
-/// [`azip!()`](macro.azip.html).
+/// Producers are used as a arguments to [`Zip`](crate::Zip) and
+/// [`azip!()`].
 ///
 /// # Comparison to `IntoIterator`
 ///
@@ -51,7 +51,7 @@ where
 /// (`AxisIter` traverses a one dimensional sequence, along an axis, while
 /// *producing* multidimensional items).
 ///
-/// See also [`IntoNdProducer`](trait.IntoNdProducer.html)
+/// See also [`IntoNdProducer`]
 pub trait NdProducer {
     /// The element produced per iteration.
     type Item;
diff --git a/src/zip/zipmacro.rs b/src/zip/zipmacro.rs
index f25e33d3..0bbe956b 100644
--- a/src/zip/zipmacro.rs
+++ b/src/zip/zipmacro.rs
@@ -1,7 +1,7 @@
 /// Array zip macro: lock step function application across several arrays and
 /// producers.
 ///
-/// This is a shorthand for [`Zip`](struct.Zip.html).
+/// This is a shorthand for [`Zip`](crate::Zip).
 ///
 /// This example:
 ///
@@ -42,60 +42,57 @@
 ///
 /// type M = Array2<f32>;
 ///
-/// fn main() {
-///     // Setup example arrays
-///     let mut a = M::zeros((16, 16));
-///     let mut b = M::zeros(a.dim());
-///     let mut c = M::zeros(a.dim());
+/// // Setup example arrays
+/// let mut a = M::zeros((16, 16));
+/// let mut b = M::zeros(a.dim());
+/// let mut c = M::zeros(a.dim());
 ///
-///     // assign values
-///     b.fill(1.);
-///     for ((i, j), elt) in c.indexed_iter_mut() {
-///         *elt = (i + 10 * j) as f32;
-///     }
-///
-///     // Example 1: Compute a simple ternary operation:
-///     // elementwise addition of b and c, stored in a
-///     azip!((a in &mut a, &b in &b, &c in &c) *a = b + c);
+/// // assign values
+/// b.fill(1.);
+/// for ((i, j), elt) in c.indexed_iter_mut() {
+///     *elt = (i + 10 * j) as f32;
+/// }
 ///
-///     assert_eq!(a, &b + &c);
+/// // Example 1: Compute a simple ternary operation:
+/// // elementwise addition of b and c, stored in a
+/// azip!((a in &mut a, &b in &b, &c in &c) *a = b + c);
 ///
-///     // Example 2: azip!() with index
-///     azip!((index (i, j), &b in &b, &c in &c) {
-///         a[[i, j]] = b - c;
-///     });
+/// assert_eq!(a, &b + &c);
 ///
-///     assert_eq!(a, &b - &c);
+/// // Example 2: azip!() with index
+/// azip!((index (i, j), &b in &b, &c in &c) {
+///     a[[i, j]] = b - c;
+/// });
 ///
+/// assert_eq!(a, &b - &c);
 ///
-///     // Example 3: azip!() on references
-///     // See the definition of the function below
-///     borrow_multiply(&mut a, &b, &c);
 ///
-///     assert_eq!(a, &b * &c);
+/// // Example 3: azip!() on references
+/// // See the definition of the function below
+/// borrow_multiply(&mut a, &b, &c);
 ///
+/// assert_eq!(a, &b * &c);
 ///
-///     // Since this function borrows its inputs, the `IntoNdProducer`
-///     // expressions don't need to explicitly include `&mut` or `&`.
-///     fn borrow_multiply(a: &mut M, b: &M, c: &M) {
-///         azip!((a in a, &b in b, &c in c) *a = b * c);
-///     }
 ///
+/// // Since this function borrows its inputs, the `IntoNdProducer`
+/// // expressions don't need to explicitly include `&mut` or `&`.
+/// fn borrow_multiply(a: &mut M, b: &M, c: &M) {
+///     azip!((a in a, &b in b, &c in c) *a = b * c);
+/// }
 ///
-///     // Example 4: using azip!() without dereference in pattern.
-///     //
-///     // Create a new array `totals` with one entry per row of `a`.
-///     // Use azip to traverse the rows of `a` and assign to the corresponding
-///     // entry in `totals` with the sum across each row.
-///     //
-///     // The row is an array view; it doesn't need to be dereferenced.
-///     let mut totals = Array1::zeros(a.nrows());
-///     azip!((totals in &mut totals, row in a.rows()) *totals = row.sum());
 ///
-///     // Check the result against the built in `.sum_axis()` along axis 1.
-///     assert_eq!(totals, a.sum_axis(Axis(1)));
-/// }
+/// // Example 4: using azip!() without dereference in pattern.
+/// //
+/// // Create a new array `totals` with one entry per row of `a`.
+/// // Use azip to traverse the rows of `a` and assign to the corresponding
+/// // entry in `totals` with the sum across each row.
+/// //
+/// // The row is an array view; it doesn't need to be dereferenced.
+/// let mut totals = Array1::zeros(a.nrows());
+/// azip!((totals in &mut totals, row in a.rows()) *totals = row.sum());
 ///
+/// // Check the result against the built in `.sum_axis()` along axis 1.
+/// assert_eq!(totals, a.sum_axis(Axis(1)));
 /// ```
 #[macro_export]
 macro_rules! azip {


Test patch:
diff --git a/tests/array.rs b/tests/array.rs
index edd58adb..d0fc67de 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -990,8 +990,22 @@ fn map1() {
 }
 
 #[test]
-fn as_slice_memory_order() {
-    // test that mutation breaks sharing
+fn mapv_into_any_same_type() {
+    let a: Array<f64, _> = array![[1., 2., 3.], [4., 5., 6.]];
+    let a_plus_one: Array<f64, _> = array![[2., 3., 4.], [5., 6., 7.]];
+    assert_eq!(a.mapv_into_any(|a| a + 1.), a_plus_one);
+}
+
+#[test]
+fn mapv_into_any_diff_types() {
+    let a: Array<f64, _> = array![[1., 2., 3.], [4., 5., 6.]];
+    let a_even: Array<bool, _> = array![[false, true, false], [true, false, true]];
+    assert_eq!(a.mapv_into_any(|a| a.round() as i32 % 2 == 0), a_even);
+}
+
+#[test]
+fn as_slice_memory_order_mut_arcarray() {
+    // Test that mutation breaks sharing for `ArcArray`.
     let a = rcarr2(&[[1., 2.], [3., 4.0f32]]);
     let mut b = a.clone();
     for elt in b.as_slice_memory_order_mut().unwrap() {
@@ -1000,6 +1014,38 @@ fn as_slice_memory_order() {
     assert!(a != b, "{:?} != {:?}", a, b);
 }
 
+#[test]
+fn as_slice_memory_order_mut_cowarray() {
+    // Test that mutation breaks sharing for `CowArray`.
+    let a = arr2(&[[1., 2.], [3., 4.0f32]]);
+    let mut b = CowArray::from(a.view());
+    for elt in b.as_slice_memory_order_mut().unwrap() {
+        *elt = 0.;
+    }
+    assert!(a != b, "{:?} != {:?}", a, b);
+}
+
+#[test]
+fn as_slice_memory_order_mut_contiguous_arcarray() {
+    // Test that unsharing preserves the strides in the contiguous case for `ArcArray`.
+    let a = rcarr2(&[[0, 5], [1, 6], [2, 7], [3, 8], [4, 9]]).reversed_axes();
+    let mut b = a.clone().slice_move(s![.., ..2]);
+    assert_eq!(b.strides(), &[1, 2]);
+    b.as_slice_memory_order_mut().unwrap();
+    assert_eq!(b.strides(), &[1, 2]);
+}
+
+#[test]
+fn as_slice_memory_order_mut_contiguous_cowarray() {
+    // Test that unsharing preserves the strides in the contiguous case for `CowArray`.
+    let a = arr2(&[[0, 5], [1, 6], [2, 7], [3, 8], [4, 9]]).reversed_axes();
+    let mut b = CowArray::from(a.slice(s![.., ..2]));
+    assert!(b.is_view());
+    assert_eq!(b.strides(), &[1, 2]);
+    b.as_slice_memory_order_mut().unwrap();
+    assert_eq!(b.strides(), &[1, 2]);
+}
+
 #[test]
 fn array0_into_scalar() {
     // With this kind of setup, the `Array`'s pointer is not the same as the
@@ -1788,6 +1834,33 @@ fn map_memory_order() {
     assert_eq!(amap.strides(), v.strides());
 }
 
+#[test]
+fn map_mut_with_unsharing() {
+    // Fortran-layout `ArcArray`.
+    let a = rcarr2(&[[0, 5], [1, 6], [2, 7], [3, 8], [4, 9]]).reversed_axes();
+    assert_eq!(a.shape(), &[2, 5]);
+    assert_eq!(a.strides(), &[1, 2]);
+    assert_eq!(
+        a.as_slice_memory_order(),
+        Some(&[0, 5, 1, 6, 2, 7, 3, 8, 4, 9][..])
+    );
+
+    // Shared reference of a portion of `a`.
+    let mut b = a.clone().slice_move(s![.., ..2]);
+    assert_eq!(b.shape(), &[2, 2]);
+    assert_eq!(b.strides(), &[1, 2]);
+    assert_eq!(b.as_slice_memory_order(), Some(&[0, 5, 1, 6][..]));
+    assert_eq!(b, array![[0, 1], [5, 6]]);
+
+    // `.map_mut()` unshares the data. Earlier versions of `ndarray` failed
+    // this assertion. See #1018.
+    assert_eq!(b.map_mut(|&mut x| x + 10), array![[10, 11], [15, 16]]);
+
+    // The strides should be preserved.
+    assert_eq!(b.shape(), &[2, 2]);
+    assert_eq!(b.strides(), &[1, 2]);
+}
+
 #[test]
 fn test_view_from_shape() {
     let s = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11];
diff --git a/tests/windows.rs b/tests/windows.rs
index b0482e4b..664616f6 100644
--- a/tests/windows.rs
+++ b/tests/windows.rs
@@ -6,7 +6,7 @@
 )]
 
 use ndarray::prelude::*;
-use ndarray::Zip;
+use ndarray::{arr3, Zip};
 
 // Edge Cases for Windows iterator:
 //
@@ -79,7 +79,6 @@ fn windows_iterator_2d() {
 /// Simple test for iterating 3d-arrays via `Windows`.
 #[test]
 fn windows_iterator_3d() {
-    use ndarray::arr3;
     let a = Array::from_iter(10..37).into_shape((3, 3, 3)).unwrap();
     itertools::assert_equal(
         a.windows(Dim((2, 2, 2))),
@@ -117,6 +116,87 @@ fn test_window_zip() {
     }
 }
 
+/// Test verifies that non existant Axis results in panic
+#[test]
+#[should_panic]
+fn axis_windows_outofbound() {
+    let a = Array::from_iter(10..37).into_shape((3, 3, 3)).unwrap();
+    a.axis_windows(Axis(4), 2);
+}
+
+/// Test verifies that zero sizes results in panic
+#[test]
+#[should_panic]
+fn axis_windows_zero_size() {
+    let a = Array::from_iter(10..37).into_shape((3, 3, 3)).unwrap();
+    a.axis_windows(Axis(0), 0);
+}
+
+/// Test verifies that over sized windows yield nothing
+#[test]
+fn axis_windows_oversized() {
+    let a = Array::from_iter(10..37).into_shape((3, 3, 3)).unwrap();
+    let mut iter = a.axis_windows(Axis(2), 4).into_iter();
+    assert_eq!(iter.next(), None);
+}
+
+/// Simple test for iterating 1d-arrays via `Axis Windows`.
+#[test]
+fn test_axis_windows_1d() {
+    let a = Array::from_iter(10..20).into_shape(10).unwrap();
+
+    itertools::assert_equal(
+        a.axis_windows(Axis(0), 5),
+        vec![
+            arr1(&[10, 11, 12, 13, 14]),
+            arr1(&[11, 12, 13, 14, 15]),
+            arr1(&[12, 13, 14, 15, 16]),
+            arr1(&[13, 14, 15, 16, 17]),
+            arr1(&[14, 15, 16, 17, 18]),
+            arr1(&[15, 16, 17, 18, 19]),
+        ],
+    );
+}
+
+/// Simple test for iterating 2d-arrays via `Axis Windows`.
+#[test]
+fn test_axis_windows_2d() {
+    let a = Array::from_iter(10..30).into_shape((5, 4)).unwrap();
+
+    itertools::assert_equal(
+        a.axis_windows(Axis(0), 2),
+        vec![
+            arr2(&[[10, 11, 12, 13], [14, 15, 16, 17]]),
+            arr2(&[[14, 15, 16, 17], [18, 19, 20, 21]]),
+            arr2(&[[18, 19, 20, 21], [22, 23, 24, 25]]),
+            arr2(&[[22, 23, 24, 25], [26, 27, 28, 29]]),
+        ],
+    );
+}
+
+/// Simple test for iterating 3d-arrays via `Axis Windows`.
+#[test]
+fn test_axis_windows_3d() {
+    let a = Array::from_iter(0..27).into_shape((3, 3, 3)).unwrap();
+
+    itertools::assert_equal(
+        a.axis_windows(Axis(1), 2),
+        vec![
+            arr3(&[
+                [[0, 1, 2], [3, 4, 5]],
+                [[9, 10, 11], [12, 13, 14]],
+                [[18, 19, 20], [21, 22, 23]],
+            ]),
+            arr3(&[
+                [[3, 4, 5], [6, 7, 8]],
+                [[12, 13, 14], [15, 16, 17]],
+                [[21, 22, 23], [24, 25, 26]],
+            ]),
+        ],
+    );
+}
+
+
 #[test]
 fn test_window_neg_stride() {
     let array = Array::from_iter(1..10).into_shape((3, 3)).unwrap();
diff --git a/xtest-blas/tests/oper.rs b/xtest-blas/tests/oper.rs
index 0aeb4768..95ab4ebb 100644
--- a/xtest-blas/tests/oper.rs
+++ b/xtest-blas/tests/oper.rs
@@ -21,6 +21,25 @@ fn mat_vec_product_1d() {
     assert_eq!(a.t().dot(&b), ans);
 }
 
+#[test]
+fn mat_vec_product_1d_broadcast() {
+    let a = arr2(&[[1.], [2.], [3.]]);
+    let b = arr1(&[1.]);
+    let b = b.broadcast(3).unwrap();
+    let ans = arr1(&[6.]);
+    assert_eq!(a.t().dot(&b), ans);
+}
+
+#[test]
+fn mat_vec_product_1d_inverted_axis() {
+    let a = arr2(&[[1.], [2.], [3.]]);
+    let mut b = arr1(&[1., 2., 3.]);
+    b.invert_axis(Axis(0));
+
+    let ans = arr1(&[3. + 4. + 3.]);
+    assert_eq!(a.t().dot(&b), ans);
+}
+
 fn range_mat(m: Ix, n: Ix) -> Array2<f32> {
     Array::linspace(0., (m * n) as f32 - 1., m * n)
         .into_shape((m, n))