Repository: rust-ndarray/ndarray
Pull request #1025 (2021-06-05 01:08:17)

Add approx-0_5 feature for approx 0.5 support
`approx` 0.5 updates its `num-complex` dependency to 0.4, which is the same version of `num-complex` which `ndarray` depends on. So, it's very useful for `ndarray` to support `approx` 0.5.

See also #1024, which attempts to support `approx` 0.5 by changing the version range to `">=0.4, <0.6"`. Unfortunately, that doesn't work.


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..56de806b 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\"",
@@ -34,6 +34,7 @@ num-complex = { version = "0.4", default-features = false }
 rayon_ = { version = "1.0.3", optional = true, package = "rayon" }
 
 approx = { version = "0.4", optional = true , default-features = false }
+approx-0_5 = { package = "approx", version = "0.5", optional = true , default-features = false }
 
 # Use via the `blas` crate feature!
 cblas-sys = { version = "0.1.4", optional = true, default-features = false }
@@ -46,7 +47,7 @@ rawpointer = { version = "0.2" }
 
 [dev-dependencies]
 defmac = "0.2"
-quickcheck = { version = "0.9", default-features = false }
+quickcheck = { version = "1.0", default-features = false }
 approx = "0.4"
 itertools = { version = "0.10.0", default-features = false, features = ["use_std"] }
 
@@ -64,7 +65,7 @@ serde-1 = ["serde"]
 test = []
 
 # This feature is used for docs
-docs = ["approx", "serde", "rayon"]
+docs = ["approx", "approx-0_5", "serde", "rayon"]
 
 std = ["num-traits/std", "matrixmultiply/std"]
 rayon = ["rayon_", "std"]
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/Cargo.toml b/ndarray-rand/Cargo.toml
index 4127d563..4e9bf265 100644
--- a/ndarray-rand/Cargo.toml
+++ b/ndarray-rand/Cargo.toml
@@ -16,7 +16,7 @@ keywords = ["multidimensional", "matrix", "rand", "ndarray"]
 [dependencies]
 ndarray = { version = "0.15", path = ".." }
 rand_distr = "0.4.0"
-quickcheck = { version = "0.9", default-features = false, optional = true }
+quickcheck = { version = "1.0", default-features = false, optional = true }
 
 [dependencies.rand]
 version = "0.8.0"
@@ -24,7 +24,7 @@ features = ["small_rng"]
 
 [dev-dependencies]
 rand_isaac = "0.3.0"
-quickcheck = { version = "0.9", default-features = false }
+quickcheck = { version = "1.0", default-features = false }
 
 [package.metadata.release]
 no-dev-version = true
diff --git a/ndarray-rand/src/lib.rs b/ndarray-rand/src/lib.rs
index e293eb56..448eb10d 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,
@@ -304,7 +304,7 @@ pub enum SamplingStrategy {
 // `Arbitrary` enables `quickcheck` to generate random `SamplingStrategy` values for testing.
 #[cfg(feature = "quickcheck")]
 impl Arbitrary for SamplingStrategy {
-    fn arbitrary<G: Gen>(g: &mut G) -> Self {
+    fn arbitrary(g: &mut Gen) -> Self {
         if bool::arbitrary(g) {
             SamplingStrategy::WithReplacement
         } else {
diff --git a/src/array_approx.rs b/src/array_approx.rs
index bb6804a2..a40982a5 100644
--- a/src/array_approx.rs
+++ b/src/array_approx.rs
@@ -1,181 +1,202 @@
-use crate::imp_prelude::*;
-use crate::Zip;
-use approx::{AbsDiffEq, RelativeEq, UlpsEq};
-
-/// **Requires crate feature `"approx"`**
-impl<A, B, S, S2, D> AbsDiffEq<ArrayBase<S2, D>> for ArrayBase<S, D>
-where
-    A: AbsDiffEq<B>,
-    A::Epsilon: Clone,
-    S: Data<Elem = A>,
-    S2: Data<Elem = B>,
-    D: Dimension,
-{
-    type Epsilon = A::Epsilon;
-
-    fn default_epsilon() -> A::Epsilon {
-        A::default_epsilon()
-    }
+#[cfg(feature = "approx")]
+mod approx_methods {
+    use crate::imp_prelude::*;
 
-    fn abs_diff_eq(&self, other: &ArrayBase<S2, D>, epsilon: A::Epsilon) -> bool {
-        if self.shape() != other.shape() {
-            return false;
+    impl<A, S, D> ArrayBase<S, D>
+    where
+        S: Data<Elem = A>,
+        D: Dimension,
+    {
+        /// A test for equality that uses the elementwise absolute difference to compute the
+        /// approximate equality of two arrays.
+        ///
+        /// **Requires crate feature `"approx"`**
+        pub fn abs_diff_eq<S2>(&self, other: &ArrayBase<S2, D>, epsilon: A::Epsilon) -> bool
+        where
+            A: ::approx::AbsDiffEq<S2::Elem>,
+            A::Epsilon: Clone,
+            S2: Data,
+        {
+            <Self as ::approx::AbsDiffEq<_>>::abs_diff_eq(self, other, epsilon)
         }
 
-        Zip::from(self)
-            .and(other)
-            .all(move |a, b| A::abs_diff_eq(a, b, epsilon.clone()))
-    }
-}
-
-/// **Requires crate feature `"approx"`**
-impl<A, B, S, S2, D> RelativeEq<ArrayBase<S2, D>> for ArrayBase<S, D>
-where
-    A: RelativeEq<B>,
-    A::Epsilon: Clone,
-    S: Data<Elem = A>,
-    S2: Data<Elem = B>,
-    D: Dimension,
-{
-    fn default_max_relative() -> A::Epsilon {
-        A::default_max_relative()
-    }
-
-    fn relative_eq(
-        &self,
-        other: &ArrayBase<S2, D>,
-        epsilon: A::Epsilon,
-        max_relative: A::Epsilon,
-    ) -> bool {
-        if self.shape() != other.shape() {
-            return false;
+        /// A test for equality that uses an elementwise relative comparison if the values are far
+        /// apart; and the absolute difference otherwise.
+        ///
+        /// **Requires crate feature `"approx"`**
+        pub fn relative_eq<S2>(
+            &self,
+            other: &ArrayBase<S2, D>,
+            epsilon: A::Epsilon,
+            max_relative: A::Epsilon,
+        ) -> bool
+        where
+            A: ::approx::RelativeEq<S2::Elem>,
+            A::Epsilon: Clone,
+            S2: Data,
+        {
+            <Self as ::approx::RelativeEq<_>>::relative_eq(self, other, epsilon, max_relative)
         }
-
-        Zip::from(self)
-            .and(other)
-            .all(move |a, b| A::relative_eq(a, b, epsilon.clone(), max_relative.clone()))
     }
 }
 
-/// **Requires crate feature `"approx"`**
-impl<A, B, S, S2, D> UlpsEq<ArrayBase<S2, D>> for ArrayBase<S, D>
-where
-    A: UlpsEq<B>,
-    A::Epsilon: Clone,
-    S: Data<Elem = A>,
-    S2: Data<Elem = B>,
-    D: Dimension,
-{
-    fn default_max_ulps() -> u32 {
-        A::default_max_ulps()
-    }
-
-    fn ulps_eq(&self, other: &ArrayBase<S2, D>, epsilon: A::Epsilon, max_ulps: u32) -> bool {
-        if self.shape() != other.shape() {
-            return false;
+macro_rules! impl_approx_traits {
+    ($approx:ident, $doc:expr) => {
+        mod $approx {
+            use crate::imp_prelude::*;
+            use crate::Zip;
+            use $approx::{AbsDiffEq, RelativeEq, UlpsEq};
+
+            #[doc = $doc]
+            impl<A, B, S, S2, D> AbsDiffEq<ArrayBase<S2, D>> for ArrayBase<S, D>
+            where
+                A: AbsDiffEq<B>,
+                A::Epsilon: Clone,
+                S: Data<Elem = A>,
+                S2: Data<Elem = B>,
+                D: Dimension,
+            {
+                type Epsilon = A::Epsilon;
+
+                fn default_epsilon() -> A::Epsilon {
+                    A::default_epsilon()
+                }
+
+                fn abs_diff_eq(&self, other: &ArrayBase<S2, D>, epsilon: A::Epsilon) -> bool {
+                    if self.shape() != other.shape() {
+                        return false;
+                    }
+
+                    Zip::from(self)
+                        .and(other)
+                        .all(move |a, b| A::abs_diff_eq(a, b, epsilon.clone()))
+                }
+            }
+
+            #[doc = $doc]
+            impl<A, B, S, S2, D> RelativeEq<ArrayBase<S2, D>> for ArrayBase<S, D>
+            where
+                A: RelativeEq<B>,
+                A::Epsilon: Clone,
+                S: Data<Elem = A>,
+                S2: Data<Elem = B>,
+                D: Dimension,
+            {
+                fn default_max_relative() -> A::Epsilon {
+                    A::default_max_relative()
+                }
+
+                fn relative_eq(
+                    &self,
+                    other: &ArrayBase<S2, D>,
+                    epsilon: A::Epsilon,
+                    max_relative: A::Epsilon,
+                ) -> bool {
+                    if self.shape() != other.shape() {
+                        return false;
+                    }
+
+                    Zip::from(self).and(other).all(move |a, b| {
+                        A::relative_eq(a, b, epsilon.clone(), max_relative.clone())
+                    })
+                }
+            }
+
+            #[doc = $doc]
+            impl<A, B, S, S2, D> UlpsEq<ArrayBase<S2, D>> for ArrayBase<S, D>
+            where
+                A: UlpsEq<B>,
+                A::Epsilon: Clone,
+                S: Data<Elem = A>,
+                S2: Data<Elem = B>,
+                D: Dimension,
+            {
+                fn default_max_ulps() -> u32 {
+                    A::default_max_ulps()
+                }
+
+                fn ulps_eq(
+                    &self,
+                    other: &ArrayBase<S2, D>,
+                    epsilon: A::Epsilon,
+                    max_ulps: u32,
+                ) -> bool {
+                    if self.shape() != other.shape() {
+                        return false;
+                    }
+
+                    Zip::from(self)
+                        .and(other)
+                        .all(move |a, b| A::ulps_eq(a, b, epsilon.clone(), max_ulps))
+                }
+            }
+
+            #[cfg(test)]
+            mod tests {
+                use crate::prelude::*;
+                use alloc::vec;
+                use $approx::{
+                    assert_abs_diff_eq, assert_abs_diff_ne, assert_relative_eq, assert_relative_ne,
+                    assert_ulps_eq, assert_ulps_ne,
+                };
+
+                #[test]
+                fn abs_diff_eq() {
+                    let a: Array2<f32> = array![[0., 2.], [-0.000010001, 100000000.]];
+                    let mut b: Array2<f32> = array![[0., 1.], [-0.000010002, 100000001.]];
+                    assert_abs_diff_ne!(a, b);
+                    b[(0, 1)] = 2.;
+                    assert_abs_diff_eq!(a, b);
+
+                    // Check epsilon.
+                    assert_abs_diff_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
+                    assert_abs_diff_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
+
+                    // Make sure we can compare different shapes without failure.
+                    let c = array![[1., 2.]];
+                    assert_abs_diff_ne!(a, c);
+                }
+
+                #[test]
+                fn relative_eq() {
+                    let a: Array2<f32> = array![[1., 2.], [-0.000010001, 100000000.]];
+                    let mut b: Array2<f32> = array![[1., 1.], [-0.000010002, 100000001.]];
+                    assert_relative_ne!(a, b);
+                    b[(0, 1)] = 2.;
+                    assert_relative_eq!(a, b);
+
+                    // Check epsilon.
+                    assert_relative_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
+                    assert_relative_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
+
+                    // Make sure we can compare different shapes without failure.
+                    let c = array![[1., 2.]];
+                    assert_relative_ne!(a, c);
+                }
+
+                #[test]
+                fn ulps_eq() {
+                    let a: Array2<f32> = array![[1., 2.], [-0.000010001, 100000000.]];
+                    let mut b: Array2<f32> = array![[1., 1.], [-0.000010002, 100000001.]];
+                    assert_ulps_ne!(a, b);
+                    b[(0, 1)] = 2.;
+                    assert_ulps_eq!(a, b);
+
+                    // Check epsilon.
+                    assert_ulps_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
+                    assert_ulps_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
+
+                    // Make sure we can compare different shapes without failure.
+                    let c = array![[1., 2.]];
+                    assert_ulps_ne!(a, c);
+                }
+            }
         }
-
-        Zip::from(self)
-            .and(other)
-            .all(move |a, b| A::ulps_eq(a, b, epsilon.clone(), max_ulps))
-    }
-}
-
-impl<A, S, D> ArrayBase<S, D>
-where
-    S: Data<Elem = A>,
-    D: Dimension,
-{
-    /// A test for equality that uses the elementwise absolute difference to compute the
-    /// approximate equality of two arrays.
-    ///
-    /// **Requires crate feature `"approx"`**
-    pub fn abs_diff_eq<S2>(&self, other: &ArrayBase<S2, D>, epsilon: A::Epsilon) -> bool
-    where
-        A: AbsDiffEq<S2::Elem>,
-        A::Epsilon: Clone,
-        S2: Data,
-    {
-        <Self as AbsDiffEq<_>>::abs_diff_eq(self, other, epsilon)
-    }
-
-    /// A test for equality that uses an elementwise relative comparison if the values are far
-    /// apart; and the absolute difference otherwise.
-    ///
-    /// **Requires crate feature `"approx"`**
-    pub fn relative_eq<S2>(
-        &self,
-        other: &ArrayBase<S2, D>,
-        epsilon: A::Epsilon,
-        max_relative: A::Epsilon,
-    ) -> bool
-    where
-        A: RelativeEq<S2::Elem>,
-        A::Epsilon: Clone,
-        S2: Data
-    {
-        <Self as RelativeEq<_>>::relative_eq(self, other, epsilon, max_relative)
-    }
-}
-
-
-#[cfg(test)]
-mod tests {
-    use crate::prelude::*;
-    use alloc::vec;
-    use approx::{
-        assert_abs_diff_eq, assert_abs_diff_ne, assert_relative_eq, assert_relative_ne,
-        assert_ulps_eq, assert_ulps_ne,
     };
+}
 
-    #[test]
-    fn abs_diff_eq() {
-        let a: Array2<f32> = array![[0., 2.], [-0.000010001, 100000000.]];
-        let mut b: Array2<f32> = array![[0., 1.], [-0.000010002, 100000001.]];
-        assert_abs_diff_ne!(a, b);
-        b[(0, 1)] = 2.;
-        assert_abs_diff_eq!(a, b);
-
-        // Check epsilon.
-        assert_abs_diff_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
-        assert_abs_diff_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
-
-        // Make sure we can compare different shapes without failure.
-        let c = array![[1., 2.]];
-        assert_abs_diff_ne!(a, c);
-    }
+#[cfg(feature = "approx")]
+impl_approx_traits!(approx, "**Requires crate feature `\"approx\"`.**");
 
-    #[test]
-    fn relative_eq() {
-        let a: Array2<f32> = array![[1., 2.], [-0.000010001, 100000000.]];
-        let mut b: Array2<f32> = array![[1., 1.], [-0.000010002, 100000001.]];
-        assert_relative_ne!(a, b);
-        b[(0, 1)] = 2.;
-        assert_relative_eq!(a, b);
-
-        // Check epsilon.
-        assert_relative_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
-        assert_relative_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
-
-        // Make sure we can compare different shapes without failure.
-        let c = array![[1., 2.]];
-        assert_relative_ne!(a, c);
-    }
-
-    #[test]
-    fn ulps_eq() {
-        let a: Array2<f32> = array![[1., 2.], [-0.000010001, 100000000.]];
-        let mut b: Array2<f32> = array![[1., 1.], [-0.000010002, 100000001.]];
-        assert_ulps_ne!(a, b);
-        b[(0, 1)] = 2.;
-        assert_ulps_eq!(a, b);
-
-        // Check epsilon.
-        assert_ulps_eq!(array![0.0f32], array![1e-40f32], epsilon = 1e-40f32);
-        assert_ulps_ne!(array![0.0f32], array![1e-40f32], epsilon = 1e-41f32);
-
-        // Make sure we can compare different shapes without failure.
-        let c = array![[1., 2.]];
-        assert_ulps_ne!(a, c);
-    }
-}
+#[cfg(feature = "approx-0_5")]
+impl_approx_traits!(approx_0_5, "**Requires crate feature `\"approx-0_5\"`.**");
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/mod.rs b/src/dimension/mod.rs
index 7bb92892..f0752142 100644
--- a/src/dimension/mod.rs
+++ b/src/dimension/mod.rs
@@ -920,12 +920,16 @@ mod test {
     }
 
     quickcheck! {
-        fn extended_gcd_solves_eq(a: isize, b: isize) -> bool {
+        // FIXME: This test can't handle larger values at the moment
+        fn extended_gcd_solves_eq(a: i16, b: i16) -> bool {
+            let (a, b) = (a as isize, b as isize);
             let (g, (x, y)) = extended_gcd(a, b);
             a * x + b * y == g
         }
 
-        fn extended_gcd_correct_gcd(a: isize, b: isize) -> bool {
+        // FIXME: This test can't handle larger values at the moment
+        fn extended_gcd_correct_gcd(a: i16, b: i16) -> bool {
+            let (a, b) = (a as isize, b as isize);
             let (g, _) = extended_gcd(a, b);
             g == gcd(a, b)
         }
@@ -941,9 +945,12 @@ mod test {
     }
 
     quickcheck! {
+        // FIXME: This test can't handle larger values at the moment
         fn solve_linear_diophantine_eq_solution_existence(
-            a: isize, b: isize, c: isize
+            a: i16, b: i16, c: i16
         ) -> TestResult {
+            let (a, b, c) = (a as isize, b as isize, c as isize);
+
             if a == 0 || b == 0 {
                 TestResult::discard()
             } else {
@@ -953,9 +960,12 @@ mod test {
             }
         }
 
+        // FIXME: This test can't handle larger values at the moment
         fn solve_linear_diophantine_eq_correct_solution(
-            a: isize, b: isize, c: isize, t: isize
+            a: i8, b: i8, c: i8, t: i8
         ) -> TestResult {
+            let (a, b, c, t) = (a as isize, b as isize, c as isize, t as isize);
+
             if a == 0 || b == 0 {
                 TestResult::discard()
             } else {
@@ -972,17 +982,23 @@ mod test {
     }
 
     quickcheck! {
+        // FIXME: This test is extremely slow, even with i16 values, investigate
         fn arith_seq_intersect_correct(
-            first1: isize, len1: isize, step1: isize,
-            first2: isize, len2: isize, step2: isize
+            first1: i8, len1: i8, step1: i8,
+            first2: i8, len2: i8, step2: i8
         ) -> TestResult {
             use std::cmp;
 
+            let (len1, len2) = (len1 as isize, len2 as isize);
+            let (first1, step1) = (first1 as isize, step1 as isize);
+            let (first2, step2) = (first2 as isize, step2 as isize);
+
             if len1 == 0 || len2 == 0 {
                 // This case is impossible to reach in `arith_seq_intersect()`
                 // because the `min*` and `max*` arguments are inclusive.
                 return TestResult::discard();
             }
+
             let len1 = len1.abs();
             let len2 = len2.abs();
 
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..22d5c78b 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,
@@ -53,3 +51,34 @@ where
         }
     }
 }
+
+impl<'a, A, Slice: ?Sized> From<&'a Slice> for CowArray<'a, A, Ix1>
+where
+    Slice: AsRef<[A]>,
+{
+    /// Create a one-dimensional clone-on-write view of the data in `slice`.
+    ///
+    /// **Panics** if the slice length is greater than [`isize::MAX`].
+    ///
+    /// ```
+    /// use ndarray::{array, CowArray};
+    ///
+    /// let array = CowArray::from(&[1., 2., 3., 4.]);
+    /// assert!(array.is_view());
+    /// assert_eq!(array, array![1., 2., 3., 4.]);
+    /// ```
+    fn from(slice: &'a Slice) -> Self {
+        Self::from(ArrayView1::from(slice))
+    }
+}
+
+impl<'a, A, S, D> From<&'a ArrayBase<S, D>> for CowArray<'a, A, D>
+where
+    S: Data<Elem = A>,
+    D: Dimension,
+{
+    /// Create a read-only clone-on-write view of the array.
+    fn from(array: &'a ArrayBase<S, D>) -> Self {
+        Self::from(array.view())
+    }
+}
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..425af96c 100644
--- a/src/impl_raw_views.rs
+++ b/src/impl_raw_views.rs
@@ -1,3 +1,4 @@
+use num_complex::Complex;
 use std::mem;
 use std::ptr::NonNull;
 
@@ -75,7 +76,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();
             }
@@ -149,6 +150,73 @@ where
     }
 }
 
+impl<T, D> RawArrayView<Complex<T>, D>
+where
+    D: Dimension,
+{
+    /// Splits the view into views of the real and imaginary components of the
+    /// elements.
+    pub fn split_re_im(self) -> Complex<RawArrayView<T, D>> {
+        // Check that the size and alignment of `Complex<T>` are as expected.
+        // These assertions should always pass, for arbitrary `T`.
+        assert_eq!(
+            mem::size_of::<Complex<T>>(),
+            mem::size_of::<T>().checked_mul(2).unwrap()
+        );
+        assert_eq!(mem::align_of::<Complex<T>>(), mem::align_of::<T>());
+
+        let dim = self.dim.clone();
+
+        // Double the strides. In the zero-sized element case and for axes of
+        // length <= 1, we leave the strides as-is to avoid possible overflow.
+        let mut strides = self.strides.clone();
+        if mem::size_of::<T>() != 0 {
+            for ax in 0..strides.ndim() {
+                if dim[ax] > 1 {
+                    strides[ax] = (strides[ax] as isize * 2) as usize;
+                }
+            }
+        }
+
+        let ptr_re: *mut T = self.ptr.as_ptr().cast();
+        let ptr_im: *mut T = if self.is_empty() {
+            // In the empty case, we can just reuse the existing pointer since
+            // it won't be dereferenced anyway. It is not safe to offset by
+            // one, since the allocation may be empty.
+            ptr_re
+        } else {
+            // In the nonempty case, we can safely offset into the first
+            // (complex) element.
+            unsafe { ptr_re.add(1) }
+        };
+
+        // `Complex` is `repr(C)` with only fields `re: T` and `im: T`. So, the
+        // real components of the elements start at the same pointer, and the
+        // imaginary components start at the pointer offset by one, with
+        // exactly double the strides. The new, doubled strides still meet the
+        // overflow constraints:
+        //
+        // - For the zero-sized element case, the strides are unchanged in
+        //   units of bytes and in units of the element type.
+        //
+        // - For the nonzero-sized element case:
+        //
+        //   - In units of bytes, the strides are unchanged. The only exception
+        //     is axes of length <= 1, but those strides are irrelevant anyway.
+        //
+        //   - Since `Complex<T>` for nonzero `T` is always at least 2 bytes,
+        //     and the original strides did not overflow in units of bytes, we
+        //     know that the new, doubled strides will not overflow in units of
+        //     `T`.
+        unsafe {
+            Complex {
+                re: RawArrayView::new_(ptr_re, dim.clone(), strides.clone()),
+                im: RawArrayView::new_(ptr_im, dim, strides),
+            }
+        }
+    }
+}
+
 impl<A, D> RawArrayViewMut<A, D>
 where
     D: Dimension,
@@ -217,7 +285,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();
             }
@@ -300,3 +368,20 @@ where
         unsafe { RawArrayViewMut::new(ptr, self.dim, self.strides) }
     }
 }
+
+impl<T, D> RawArrayViewMut<Complex<T>, D>
+where
+    D: Dimension,
+{
+    /// Splits the view into views of the real and imaginary components of the
+    /// elements.
+    pub fn split_re_im(self) -> Complex<RawArrayViewMut<T, D>> {
+        let Complex { re, im } = self.into_raw_view().split_re_im();
+        unsafe {
+            Complex {
+                re: RawArrayViewMut::new(re.ptr, re.dim, re.strides),
+                im: RawArrayViewMut::new(im.ptr, im.dim, im.strides),
+            }
+        }
+    }
+}
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..84d24038 100644
--- a/src/impl_views/splitting.rs
+++ b/src/impl_views/splitting.rs
@@ -8,6 +8,7 @@
 
 use crate::imp_prelude::*;
 use crate::slice::MultiSliceArg;
+use num_complex::Complex;
 
 /// Methods for read-only array views.
 impl<'a, A, D> ArrayView<'a, A, D>
@@ -95,6 +96,37 @@ where
     }
 }
 
+impl<'a, T, D> ArrayView<'a, Complex<T>, D>
+where
+    D: Dimension,
+{
+    /// Splits the view into views of the real and imaginary components of the
+    /// elements.
+    ///
+    /// ```
+    /// use ndarray::prelude::*;
+    /// use num_complex::{Complex, Complex64};
+    ///
+    /// let arr = array![
+    ///     [Complex64::new(1., 2.), Complex64::new(3., 4.)],
+    ///     [Complex64::new(5., 6.), Complex64::new(7., 8.)],
+    ///     [Complex64::new(9., 10.), Complex64::new(11., 12.)],
+    /// ];
+    /// let Complex { re, im } = arr.view().split_re_im();
+    /// assert_eq!(re, array![[1., 3.], [5., 7.], [9., 11.]]);
+    /// assert_eq!(im, array![[2., 4.], [6., 8.], [10., 12.]]);
+    /// ```
+    pub fn split_re_im(self) -> Complex<ArrayView<'a, T, D>> {
+        unsafe {
+            let Complex { re, im } = self.into_raw_view().split_re_im();
+            Complex {
+                re: re.deref_into_view(),
+                im: im.deref_into_view(),
+            }
+        }
+    }
+}
+
 /// Methods for read-write array views.
 impl<'a, A, D> ArrayViewMut<'a, A, D>
 where
@@ -121,7 +153,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:
     ///
@@ -135,3 +167,41 @@ where
         info.multi_slice_move(self)
     }
 }
+
+impl<'a, T, D> ArrayViewMut<'a, Complex<T>, D>
+where
+    D: Dimension,
+{
+    /// Splits the view into views of the real and imaginary components of the
+    /// elements.
+    ///
+    /// ```
+    /// use ndarray::prelude::*;
+    /// use num_complex::{Complex, Complex64};
+    ///
+    /// let mut arr = array![
+    ///     [Complex64::new(1., 2.), Complex64::new(3., 4.)],
+    ///     [Complex64::new(5., 6.), Complex64::new(7., 8.)],
+    ///     [Complex64::new(9., 10.), Complex64::new(11., 12.)],
+    /// ];
+    ///
+    /// let Complex { mut re, mut im } = arr.view_mut().split_re_im();
+    /// assert_eq!(re, array![[1., 3.], [5., 7.], [9., 11.]]);
+    /// assert_eq!(im, array![[2., 4.], [6., 8.], [10., 12.]]);
+    ///
+    /// re[[0, 1]] = 13.;
+    /// im[[2, 0]] = 14.;
+    ///
+    /// assert_eq!(arr[[0, 1]], Complex64::new(13., 4.));
+    /// assert_eq!(arr[[2, 0]], Complex64::new(9., 14.));
+    /// ```
+    pub fn split_re_im(self) -> Complex<ArrayViewMut<'a, T, D>> {
+        unsafe {
+            let Complex { re, im } = self.into_raw_view_mut().split_re_im();
+            Complex {
+                re: re.deref_into_view_mut(),
+                im: im.deref_into_view_mut(),
+            }
+        }
+    }
+}
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..8cf23e91 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
 //!
@@ -89,7 +91,9 @@
 //!   - Enables parallel iterators, parallelized methods, the [`parallel`] module and [`par_azip!`].
 //!   - Implies std
 //! - `approx`
-//!   - Enables implementations of traits from the [`approx`] crate.
+//!   - Enables implementations of traits from version 0.4 of the [`approx`] crate.
+//! - `approx-0_5`
+//!   - Enables implementations of traits from version 0.5 of the [`approx`] crate.
 //! - `blas`
 //!   - Enable transparent BLAS support for matrix multiplication.
 //!     Uses ``blas-src`` for pluggable backend, which needs to be configured
@@ -99,13 +103,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 +262,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 +284,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 +303,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 +377,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 +400,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 +475,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 +494,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 +506,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 +518,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 +627,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 +719,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 +742,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 +888,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 +903,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 +1037,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 +1062,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 +1113,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 +1136,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 +1301,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 +1325,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 +1345,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 +1366,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 +1379,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 +1393,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 +1416,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 +1435,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 +1448,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 +1465,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 +1482,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>),
@@ -1616,7 +1598,7 @@ pub mod linalg;
 mod impl_ops;
 pub use crate::impl_ops::ScalarOperand;
 
-#[cfg(feature = "approx")]
+#[cfg(any(feature = "approx", feature = "approx-0_5"))]
 mod array_approx;
 
 // Array view methods
diff --git a/src/linalg/impl_linalg.rs b/src/linalg/impl_linalg.rs
index 61b91eae..0e6700d6 100644
--- a/src/linalg/impl_linalg.rs
+++ b/src/linalg/impl_linalg.rs
@@ -7,25 +7,32 @@
 // except according to those terms.
 
 use crate::imp_prelude::*;
+
+#[cfg(feature = "blas")]
+use crate::dimension::offset_from_low_addr_ptr_to_logical_ptr;
 use crate::numeric_util;
 
 use crate::{LinalgScalar, Zip};
 
 use std::any::TypeId;
+use std::mem::MaybeUninit;
 use alloc::vec::Vec;
 
+#[cfg(feature = "blas")]
+use libc::c_int;
 #[cfg(feature = "blas")]
 use std::cmp;
 #[cfg(feature = "blas")]
 use std::mem::swap;
-#[cfg(feature = "blas")]
-use libc::c_int;
 
 #[cfg(feature = "blas")]
 use cblas_sys as blas_sys;
 #[cfg(feature = "blas")]
 use cblas_sys::{CblasNoTrans, CblasRowMajor, CblasTrans, CBLAS_LAYOUT};
 
+#[cfg(feature = "blas")]
+use num_complex::{Complex32 as c32, Complex64 as c64};
+
 /// len of vector before we use blas
 #[cfg(feature = "blas")]
 const DOT_BLAS_CUTOFF: usize = 32;
@@ -374,7 +381,12 @@ fn mat_mul_impl<A>(
     // size cutoff for using BLAS
     let cut = GEMM_BLAS_CUTOFF;
     let ((mut m, a), (_, mut n)) = (lhs.dim(), rhs.dim());
-    if !(m > cut || n > cut || a > cut) || !(same_type::<A, f32>() || same_type::<A, f64>()) {
+    if !(m > cut || n > cut || a > cut)
+        || !(same_type::<A, f32>()
+        || same_type::<A, f64>()
+        || same_type::<A, c32>()
+        || same_type::<A, c64>())
+    {
         return mat_mul_general(alpha, lhs, rhs, beta, c);
     }
     {
@@ -404,8 +416,23 @@ fn mat_mul_impl<A>(
             rhs_trans = CblasTrans;
         }
 
+        macro_rules! gemm_scalar_cast {
+            (f32, $var:ident) => {
+                cast_as(&$var)
+            };
+            (f64, $var:ident) => {
+                cast_as(&$var)
+            };
+            (c32, $var:ident) => {
+                &$var as *const A as *const _
+            };
+            (c64, $var:ident) => {
+                &$var as *const A as *const _
+            };
+        }
+
         macro_rules! gemm {
-            ($ty:ty, $gemm:ident) => {
+            ($ty:tt, $gemm:ident) => {
                 if blas_row_major_2d::<$ty, _>(&lhs_)
                     && blas_row_major_2d::<$ty, _>(&rhs_)
                     && blas_row_major_2d::<$ty, _>(&c_)
@@ -425,7 +452,7 @@ fn mat_mul_impl<A>(
                     let lhs_stride = cmp::max(lhs_.strides()[0] as blas_index, k as blas_index);
                     let rhs_stride = cmp::max(rhs_.strides()[0] as blas_index, n as blas_index);
                     let c_stride = cmp::max(c_.strides()[0] as blas_index, n as blas_index);
-
+                    
                     // gemm is C ← αA^Op B^Op + βC
                     // Where Op is notrans/trans/conjtrans
                     unsafe {
@@ -433,17 +460,17 @@ fn mat_mul_impl<A>(
                             CblasRowMajor,
                             lhs_trans,
                             rhs_trans,
-                            m as blas_index,               // m, rows of Op(a)
-                            n as blas_index,               // n, cols of Op(b)
-                            k as blas_index,               // k, cols of Op(a)
-                            cast_as(&alpha),               // alpha
-                            lhs_.ptr.as_ptr() as *const _, // a
-                            lhs_stride,                    // lda
-                            rhs_.ptr.as_ptr() as *const _, // b
-                            rhs_stride,                    // ldb
-                            cast_as(&beta),                // beta
-                            c_.ptr.as_ptr() as *mut _,     // c
-                            c_stride,                      // ldc
+                            m as blas_index,                 // m, rows of Op(a)
+                            n as blas_index,                 // n, cols of Op(b)
+                            k as blas_index,                 // k, cols of Op(a)
+                            gemm_scalar_cast!($ty, alpha),   // alpha
+                            lhs_.ptr.as_ptr() as *const _,   // a
+                            lhs_stride,                      // lda
+                            rhs_.ptr.as_ptr() as *const _,   // b
+                            rhs_stride,                      // ldb
+                            gemm_scalar_cast!($ty, beta),    // beta
+                            c_.ptr.as_ptr() as *mut _,       // c
+                            c_stride,                        // ldc
                         );
                     }
                     return;
@@ -452,6 +479,9 @@ fn mat_mul_impl<A>(
         }
         gemm!(f32, cblas_sgemm);
         gemm!(f64, cblas_dgemm);
+
+        gemm!(c32, cblas_cgemm);
+        gemm!(c64, cblas_zgemm);
     }
     mat_mul_general(alpha, lhs, rhs, beta, c)
 }
@@ -600,9 +630,7 @@ pub fn general_mat_vec_mul<A, S1, S2, S3>(
     S3: DataMut<Elem = A>,
     A: LinalgScalar,
 {
-    unsafe {
-        general_mat_vec_mul_impl(alpha, a, x, beta, y.raw_view_mut())
-    }
+    unsafe { general_mat_vec_mul_impl(alpha, a, x, beta, y.raw_view_mut()) }
 }
 
 /// General matrix-vector multiplication
@@ -649,6 +677,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 +694,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;
@@ -691,6 +725,39 @@ unsafe fn general_mat_vec_mul_impl<A, S1, S2>(
     }
 }
 
+
+/// Kronecker product of 2D matrices.
+///
+/// The kronecker product of a LxN matrix A and a MxR matrix B is a (L*M)x(N*R)
+/// matrix K formed by the block multiplication A_ij * B.
+pub fn kron<A, S1, S2>(a: &ArrayBase<S1, Ix2>, b: &ArrayBase<S2, Ix2>) -> Array<A, Ix2>
+where
+    S1: Data<Elem = A>,
+    S2: Data<Elem = A>,
+    A: LinalgScalar,
+{
+    let dimar = a.shape()[0];
+    let dimac = a.shape()[1];
+    let dimbr = b.shape()[0];
+    let dimbc = b.shape()[1];
+    let mut out: Array2<MaybeUninit<A>> = Array2::uninit((
+        dimar
+            .checked_mul(dimbr)
+            .expect("Dimensions of kronecker product output array overflows usize."),
+        dimac
+            .checked_mul(dimbc)
+            .expect("Dimensions of kronecker product output array overflows usize."),
+    ));
+    Zip::from(out.exact_chunks_mut((dimbr, dimbc)))
+        .and(a)
+        .for_each(|out, &a| {
+            Zip::from(out).and(b).for_each(|out, &b| {
+                *out = MaybeUninit::new(a * b);
+            })
+        });
+    unsafe { out.assume_init() }
+}
+
 #[inline(always)]
 /// Return `true` if `A` and `B` are the same type
 fn same_type<A: 'static, B: 'static>() -> bool {
@@ -719,7 +786,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/linalg/mod.rs b/src/linalg/mod.rs
index 8575905c..abd7b2b9 100644
--- a/src/linalg/mod.rs
+++ b/src/linalg/mod.rs
@@ -11,5 +11,6 @@
 pub use self::impl_linalg::general_mat_mul;
 pub use self::impl_linalg::general_mat_vec_mul;
 pub use self::impl_linalg::Dot;
+pub use self::impl_linalg::kron;
 
 mod impl_linalg;
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/ndarray-rand/tests/tests.rs b/ndarray-rand/tests/tests.rs
index d8b18fe5..5a0cd6c1 100644
--- a/ndarray-rand/tests/tests.rs
+++ b/ndarray-rand/tests/tests.rs
@@ -5,7 +5,7 @@ use ndarray_rand::rand::{distributions::Distribution, thread_rng};
 use ndarray::ShapeBuilder;
 use ndarray_rand::rand_distr::Uniform;
 use ndarray_rand::{RandomExt, SamplingStrategy};
-use quickcheck::quickcheck;
+use quickcheck::{quickcheck, TestResult};
 
 #[test]
 fn test_dim() {
@@ -51,7 +51,8 @@ fn oversampling_without_replacement_should_panic() {
 }
 
 quickcheck! {
-    fn oversampling_with_replacement_is_fine(m: usize, n: usize) -> bool {
+    fn oversampling_with_replacement_is_fine(m: u8, n: u8) -> TestResult {
+        let (m, n) = (m as usize, n as usize);
         let a = Array::random((m, n), Uniform::new(0., 2.));
         // Higher than the length of both axes
         let n_samples = m + n + 1;
@@ -59,24 +60,28 @@ quickcheck! {
         // We don't want to deal with sampling from 0-length axes in this test
         if m != 0 {
             if !sampling_works(&a, SamplingStrategy::WithReplacement, Axis(0), n_samples) {
-                return false;
+                return TestResult::failed();
             }
+        } else {
+            return TestResult::discard();
         }
 
         // We don't want to deal with sampling from 0-length axes in this test
         if n != 0 {
             if !sampling_works(&a, SamplingStrategy::WithReplacement, Axis(1), n_samples) {
-                return false;
+                return TestResult::failed();
             }
+        } else {
+            return TestResult::discard();
         }
-
-        true
+        TestResult::passed()
     }
 }
 
 #[cfg(feature = "quickcheck")]
 quickcheck! {
-    fn sampling_behaves_as_expected(m: usize, n: usize, strategy: SamplingStrategy) -> bool {
+    fn sampling_behaves_as_expected(m: u8, n: u8, strategy: SamplingStrategy) -> TestResult {
+        let (m, n) = (m as usize, n as usize);
         let a = Array::random((m, n), Uniform::new(0., 2.));
         let mut rng = &mut thread_rng();
 
@@ -84,19 +89,23 @@ quickcheck! {
         if m != 0 {
             let n_row_samples = Uniform::from(1..m+1).sample(&mut rng);
             if !sampling_works(&a, strategy.clone(), Axis(0), n_row_samples) {
-                return false;
+                return TestResult::failed();
             }
+        } else {
+            return TestResult::discard();
         }
 
         // We don't want to deal with sampling from 0-length axes in this test
         if n != 0 {
             let n_col_samples = Uniform::from(1..n+1).sample(&mut rng);
             if !sampling_works(&a, strategy, Axis(1), n_col_samples) {
-                return false;
+                return TestResult::failed();
             }
+        } else {
+            return TestResult::discard();
         }
 
-        true
+        TestResult::passed()
     }
 }
 
diff --git a/tests/array.rs b/tests/array.rs
index edd58adb..16d90156 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -7,12 +7,14 @@
     clippy::float_cmp
 )]
 
+use approx::assert_relative_eq;
 use defmac::defmac;
 use itertools::{zip, Itertools};
 use ndarray::prelude::*;
 use ndarray::{arr3, rcarr2};
 use ndarray::indices;
 use ndarray::{Slice, SliceInfo, SliceInfoElem};
+use num_complex::Complex;
 use std::convert::TryFrom;
 
 macro_rules! assert_panics {
@@ -990,8 +992,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 +1016,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 +1836,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];
@@ -2428,3 +2503,71 @@ fn test_remove_index_oob3() {
     let mut a = array![[10], [4], [1]];
     a.remove_index(Axis(2), 0);
 }
+
+#[test]
+fn test_split_re_im_view() {
+    let a = Array3::from_shape_fn((3, 4, 5), |(i, j, k)| {
+        Complex::<f32>::new(i as f32 * j as f32, k as f32)
+    });
+    let Complex { re, im } = a.view().split_re_im();
+    assert_relative_eq!(re.sum(), 90.);
+    assert_relative_eq!(im.sum(), 120.);
+}
+
+#[test]
+fn test_split_re_im_view_roundtrip() {
+    let a_re = Array3::from_shape_fn((3,1,5), |(i, j, _k)| {
+        i * j
+    });
+    let a_im = Array3::from_shape_fn((3,1,5), |(_i, _j, k)| {
+        k
+    });
+    let a = Array3::from_shape_fn((3,1,5), |(i,j,k)| {
+        Complex::new(a_re[[i,j,k]], a_im[[i,j,k]])
+    });
+    let Complex { re, im } = a.view().split_re_im();
+    assert_eq!(a_re, re);
+    assert_eq!(a_im, im);
+}
+
+#[test]
+fn test_split_re_im_view_mut() {
+    let eye_scalar = Array2::<u32>::eye(4);
+    let eye_complex = Array2::<Complex<u32>>::eye(4);
+    let mut a = Array2::<Complex<u32>>::zeros((4, 4));
+    let Complex { mut re, im } = a.view_mut().split_re_im();
+    re.assign(&eye_scalar);
+    assert_eq!(im.sum(), 0);
+    assert_eq!(a, eye_complex);
+}
+
+#[test]
+fn test_split_re_im_zerod() {
+    let mut a = Array0::from_elem((), Complex::new(42, 32));
+    let Complex { re, im } = a.view().split_re_im();
+    assert_eq!(re.get(()), Some(&42));
+    assert_eq!(im.get(()), Some(&32));
+    let cmplx = a.view_mut().split_re_im();
+    cmplx.re.assign_to(cmplx.im);
+    assert_eq!(a.get(()).unwrap().im, 42);
+}
+
+#[test]
+fn test_split_re_im_permuted() {
+    let a = Array3::from_shape_fn((3, 4, 5), |(i, j, k)| {
+        Complex::new(i * k + j, k)
+    });
+    let permuted = a.view().permuted_axes([1,0,2]);
+    let Complex { re, im } = permuted.split_re_im();
+    assert_eq!(re.get((3,2,4)).unwrap(), &11);
+    assert_eq!(im.get((3,2,4)).unwrap(), &4);
+}
+
+#[test]
+fn test_split_re_im_invert_axis() {
+    let mut a = Array::from_shape_fn((2, 3, 2), |(i, j, k)| Complex::new(i as f64 + j as f64, i as f64 + k as f64));
+    a.invert_axis(Axis(1));
+    let cmplx = a.view().split_re_im();
+    assert_eq!(cmplx.re, a.mapv(|z| z.re));
+    assert_eq!(cmplx.im, a.mapv(|z| z.im));
+}
diff --git a/tests/oper.rs b/tests/oper.rs
index ed612bad..05172868 100644
--- a/tests/oper.rs
+++ b/tests/oper.rs
@@ -6,6 +6,7 @@
 )]
 #![cfg(feature = "std")]
 use ndarray::linalg::general_mat_mul;
+use ndarray::linalg::kron;
 use ndarray::prelude::*;
 use ndarray::{rcarr1, rcarr2};
 use ndarray::{Data, LinalgScalar};
@@ -820,3 +821,65 @@ fn vec_mat_mul() {
         }
     }
 }
+
+#[test]
+fn kron_square_f64() {
+    let a = arr2(&[[1.0, 0.0], [0.0, 1.0]]);
+    let b = arr2(&[[0.0, 1.0], [1.0, 0.0]]);
+
+    assert_eq!(
+        kron(&a, &b),
+        arr2(&[
+            [0.0, 1.0, 0.0, 0.0],
+            [1.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 1.0],
+            [0.0, 0.0, 1.0, 0.0]
+        ]),
+    );
+
+    assert_eq!(
+        kron(&b, &a),
+        arr2(&[
+            [0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 0.0, 1.0],
+            [1.0, 0.0, 0.0, 0.0],
+            [0.0, 1.0, 0.0, 0.0]
+        ]),
+    )
+}
+
+#[test]
+fn kron_square_i64() {
+    let a = arr2(&[[1, 0], [0, 1]]);
+    let b = arr2(&[[0, 1], [1, 0]]);
+
+    assert_eq!(
+        kron(&a, &b),
+        arr2(&[[0, 1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]]),
+    );
+
+    assert_eq!(
+        kron(&b, &a),
+        arr2(&[[0, 0, 1, 0], [0, 0, 0, 1], [1, 0, 0, 0], [0, 1, 0, 0]]),
+    )
+}
+
+#[test]
+fn kron_i64() {
+    let a = arr2(&[[1, 0]]);
+    let b = arr2(&[[0, 1], [1, 0]]);
+    let r = arr2(&[[0, 1, 0, 0], [1, 0, 0, 0]]);
+    assert_eq!(kron(&a, &b), r);
+
+    let a = arr2(&[[1, 0], [0, 0], [0, 1]]);
+    let b = arr2(&[[0, 1], [1, 0]]);
+    let r = arr2(&[
+        [0, 1, 0, 0],
+        [1, 0, 0, 0],
+        [0, 0, 0, 0],
+        [0, 0, 0, 0],
+        [0, 0, 0, 1],
+        [0, 0, 1, 0],
+    ]);
+    assert_eq!(kron(&a, &b), r);
+}
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/Cargo.toml b/xtest-blas/Cargo.toml
index 463be686..aee8c3ce 100644
--- a/xtest-blas/Cargo.toml
+++ b/xtest-blas/Cargo.toml
@@ -11,6 +11,7 @@ test = false
 approx = "0.4"
 defmac = "0.2"
 num-traits = "0.2"
+num-complex = { version = "0.4", default-features = false }
 
 [dependencies]
 ndarray = { path = "../", features = ["approx", "blas"] }
diff --git a/xtest-blas/tests/oper.rs b/xtest-blas/tests/oper.rs
index 0aeb4768..120626d0 100644
--- a/xtest-blas/tests/oper.rs
+++ b/xtest-blas/tests/oper.rs
@@ -1,8 +1,9 @@
 extern crate approx;
+extern crate blas_src;
 extern crate defmac;
 extern crate ndarray;
+extern crate num_complex;
 extern crate num_traits;
-extern crate blas_src;
 
 use ndarray::prelude::*;
 
@@ -12,6 +13,8 @@ use ndarray::{Data, Ix, LinalgScalar};
 
 use approx::assert_relative_eq;
 use defmac::defmac;
+use num_complex::Complex32;
+use num_complex::Complex64;
 
 #[test]
 fn mat_vec_product_1d() {
@@ -21,6 +24,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))
@@ -33,6 +55,20 @@ fn range_mat64(m: Ix, n: Ix) -> Array2<f64> {
         .unwrap()
 }
 
+fn range_mat_complex(m: Ix, n: Ix) -> Array2<Complex32> {
+    Array::linspace(0., (m * n) as f32 - 1., m * n)
+        .into_shape((m, n))
+        .unwrap()
+        .map(|&f| Complex32::new(f, 0.))
+}
+
+fn range_mat_complex64(m: Ix, n: Ix) -> Array2<Complex64> {
+    Array::linspace(0., (m * n) as f64 - 1., m * n)
+        .into_shape((m, n))
+        .unwrap()
+        .map(|&f| Complex64::new(f, 0.))
+}
+
 fn range1_mat64(m: Ix) -> Array1<f64> {
     Array::linspace(0., m as f64 - 1., m)
 }
@@ -231,6 +267,77 @@ fn gemm_64_1_f() {
     assert_relative_eq!(y, answer, epsilon = 1e-12, max_relative = 1e-7);
 }
 
+#[test]
+fn gemm_c64_1_f() {
+    let a = range_mat_complex64(64, 64).reversed_axes();
+    let (m, n) = a.dim();
+    // m x n  times n x 1  == m x 1
+    let x = range_mat_complex64(n, 1);
+    let mut y = range_mat_complex64(m, 1);
+    let answer = reference_mat_mul(&a, &x) + &y;
+    general_mat_mul(
+        Complex64::new(1.0, 0.),
+        &a,
+        &x,
+        Complex64::new(1.0, 0.),
+        &mut y,
+    );
+    assert_relative_eq!(
+        y.mapv(|i| i.norm_sqr()),
+        answer.mapv(|i| i.norm_sqr()),
+        epsilon = 1e-12,
+        max_relative = 1e-7
+    );
+}
+
+#[test]
+fn gemm_c32_1_f() {
+    let a = range_mat_complex(64, 64).reversed_axes();
+    let (m, n) = a.dim();
+    // m x n  times n x 1  == m x 1
+    let x = range_mat_complex(n, 1);
+    let mut y = range_mat_complex(m, 1);
+    let answer = reference_mat_mul(&a, &x) + &y;
+    general_mat_mul(
+        Complex32::new(1.0, 0.),
+        &a,
+        &x,
+        Complex32::new(1.0, 0.),
+        &mut y,
+    );
+    assert_relative_eq!(
+        y.mapv(|i| i.norm_sqr()),
+        answer.mapv(|i| i.norm_sqr()),
+        epsilon = 1e-12,
+        max_relative = 1e-7
+    );
+}
+
+#[test]
+fn gemm_c64_actually_complex() {
+    let mut a = range_mat_complex64(4,4);
+    a = a.map(|&i| if i.re > 8. { i.conj() } else { i });
+    let mut b = range_mat_complex64(4,6);
+    b = b.map(|&i| if i.re > 4. { i.conj() } else {i});
+    let mut y = range_mat_complex64(4,6);
+    let alpha = Complex64::new(0., 1.0);
+    let beta = Complex64::new(1.0, 1.0);
+    let answer = alpha * reference_mat_mul(&a, &b) + beta * &y;
+    general_mat_mul(
+        alpha.clone(),
+        &a,
+        &b,
+        beta.clone(),
+        &mut y,
+    );
+    assert_relative_eq!(
+        y.mapv(|i| i.norm_sqr()),
+        answer.mapv(|i| i.norm_sqr()),
+        epsilon = 1e-12,
+        max_relative = 1e-7
+    );
+}
+
 #[test]
 fn gen_mat_vec_mul() {
     let alpha = -2.3;