Repository: rust-ndarray/ndarray
Pull request #982 (2021-04-17 21:19:06)

Better reshaping, part 1: .to_shape()
New features:

- `.to_shape(impl IntoDimension) -> CowArray`
- `.to_shape((impl IntoDimension, Order)) -> CowArray`
- `to_shape` can return a view in all cases where it's technically possible to do so, saving copying. This is a big step forwards compared to previous `into_shape` and similar (and this improvement will eventually land in all the reshaping methods.)



```rust
use ndarray::{array, Order};

assert!(
    array![1., 2., 3., 4., 5., 6.].to_shape(((2, 3), Order::ColumnMajor)).unwrap()
    == array![[1., 3., 5.],
              [2., 4., 6.]]
);
```

- `to_shape` returns a view if possible otherwise copies into an owned array. Always succeeds if the number of elements is correct.


Method parameters:

- `.to_shape<Sh>(shape: Sh) where Sh: ShapeArg`
  - Takes a dimension pattern like `(1, 2, 3)` like before. This uses the default Order.
  - Also takes a tuple of dimension and order like `((1, 2, 3), Order::RowMajor)` to specify order
 
- Order enum:
  - `RowMajor` (with alias `C`; this is the default if not specified)
  - `ColumnMajor` (with alias `F`)

### Notes

It's important to have reshaping methods that can do so without any clone bound or similar. For example a reshape operation that can be used on everything from raw views, views and other arrays. `into_shape` or similar should fill this need.

###  Coming in a later change:

A version of `into_shape` that uses the `Order` parameter and that can always reshape owned arrays.

- `.into_shape() -> "Self"`  (Same array storage type)

- `into_shape` transforms the input but can only do so (for views) if the memory layout is exactly correct. For `Array` and uniquely owned `ArcArray`, it always succeeds if the number of elements is correct.

### Order::Automatic

`Order::Automatic` has been explored (and implemented) but is not part of this version of the change. No use case has been presented.

Partially addresses #390 

bluss: @jturner314 What do you think about these aspects of the shape + order argument?

1. Use the same shape + order argument in `.to_shape()` as we do in for example `Array::zeros`
2. The vehicle for ordering is going to be the `Order { RowMajor, ColumnMajor, Automatic }` enum - with aliases `Order::C` and `Order::F`
  - Automatic could be removed if we don't like it. For construction it would have to mean "default".
  - I think it is ok to "overload" the meaning of this to mean both memory layout and logical index order - we can understand those as two sides of the same coin
3. Important quesiton if the function is going to be called like **a** or **b**:
  a. `.to_shape((2, 3).order(Order::RowMajor))`
  b. `.to_shape(((2, 3), Order::RowMajor)`
  - **a** is easier to read but requires us to continue the trait that adds methods to tuples, and `.order(Order::F)` is a bit redundant.
  - **b** is simple and easy to understand - drawback is the forest of `()`'s.
  - At the moment I'm leaning towards **b**
  - For construction, custom strides needs to be factored in. Either `(2, 3).strides(..)` or `((2, 3), CustomStrides(..))` would work.
jturner314: > Use the same shape + order argument in `.to_shape()` as we do in for example `Array::zeros`

I think this is fine. While it is true that the order parameter is being used in two different ways -- `.to_shape()` uses it as "the mapping from an n-D index to a 1-D  *logical index* of the *input array and output array*" and `Array::zeros` uses it as "the mapping from an n-D index to a 1-D *memory offset* of the *output array*" -- I think that's okay.

The primary case where I could see us wanting different options for these two cases would be if we wanted to support directly specifying strides for `Array::zeros` for convenience (with a check that there are no holes). Allowing strides as a parameter for `.to_shape()` would be too confusing, IMO. However, we could handle this case in a nice way by adding a conversion from contiguous strides to the `Order` type.

Another case is that we may want to define a `ContiguousLayout` struct (with shape and strides/order fields) for various purposes, including constructors. The name `ContiguousLayout` would be confusing as a parameter for `.to_shape()`. I generally associate "layout" with "memory layout", but `.to_shape()` isn't really using the order to specify a memory layout. I don't think this is a significant concern, though.

> The vehicle for ordering can be the `enum Order { RowMajor, ColumnMajor, Automatic }` - with aliases `Order::C`, `Order::F`, `Order::A` (associated consts).

I agree that an enum works well for this. By the way, the idea of using associated consts for enum variant aliases is pretty clever. I haven't seen that before.

> * Automatic could be removed if we don't like it. For construction it would have to mean "default".

Personally, I don't think I'd ever use an `Automatic` option with reshaping, since it makes the output values depend on the memory layout of the input. NumPy has this option, though, so I'd guess that it has some users. My preference would be to leave it out for now, and add it only if someone requests it.

> * I think it is ok to "overload" the meaning of this to mean both memory layout and logical index order - we can understand those as two sides of the same coin

Yes, I think it's okay to "overload" the meaning of the order. In both cases, it represents a mapping from an n-D index to a 1-D index. It's basically equivalent to strides, except for the lack of holes in the 1-D space.

I think that it would be worth adding a more flexible `Axes` order variant which has a list of axis indices, so, for example, `[0, 1, 2]` would be C layout for 3 dimensions and `[2, 1, 0]` would be F layout. I'd actually suggest adding an `AxesOrder` struct as a wrapper around `D: Dimension` to enforce uniqueness and in-boundedness. The `Order::Axes` enum variant would wrap the struct. However, this variant isn't necessary for an initial implementation of `.to_shape()`.

We could also consider allowing axis inversions in the order, but that's less important than the "axes order" variant, since inversions don't change the shape. (It's simple to call `invert_axis` *after* constructing an array, for example. In contrast, a workaround for the lack of an `Axes` order variant requires reordering the shape passed to the constructor.) That said, it would be nice to support axis inversions in the order so that we could have an `.order()` method which would return the order of an array's elements.

> An important question is if the function is going to be called like **a** or **b**:
>
> a. `.to_shape((2, 3).order(Order::RowMajor))`
> b. `.to_shape(((2, 3), Order::RowMajor))`

First, it's worth noting that these will be written as

a. `.to_shape([2, 3].order(Order::RowMajor))`
b. `.to_shape(([2, 3], Order::RowMajor))`

once we drop support for using tuples to represent shapes, so the parentheses in the `b` case aren't too bad. I don't have a strong preference. I like `b` a bit better because depending on an extension trait implemented for standard library types is a little unusual.

I wouldn't mind requiring the order parameter for `.to_shape()`, which would simplify things for it.

We could remove support for using `[usize; N]` directly as a shape, and instead require a `Shape` struct, which would remove the need for an extension trait to add an `.order()` method. I dislike this option, though, because it's much less concise than just using `[usize; N]` directly.
bluss: - `to_shape`/`into_shape` calls are so useful and common that I want to be able to use them without explicit ordering arguments.
- We'll use lots of examples where ordering is explicitly specified, to make it clear that this is also a normal use case.
- To avoid having builder structs for so many methods, let's try to use a ShapeArg or similar trait i.e. where the diversity of the method argument sits in the argument position. The upside is that the simple case is simple - just `.to_shape((2, 3))` etc.
- I agree that we will get to a place where shapes are *canonically* specified with arrays - we will map `[usize; N]` <-> *N*-dimensional array shape one to one. However, why not keep the `IntoDimension` or similar implementations for tuples. I don't want to break all ndarray code in existance, and right now it looks like we can easily keep the tuple conversions too.
jturner314: > * `to_shape`/`into_shape` calls are so useful and common that I want to be able to use them without explicit ordering arguments.

Okay, that's fine with me. I think the default should definitely be C order rather than e.g. an "automatic" (layout-dependent) order.

> right now it looks like we can easily keep the tuple conversions too.

Yes, that will probably always be true for small ndims. I just assumed that we'd drop support for tuples at some point, because I don't see much reason to continue supporting them except for ease of upgrading to the next breaking release. Personally, I don't mind breaking changes in a breaking release as long as (1) there's at least a compiler warning (preferably an error) about the breakage and (2) it's easy to update the code. In the case of dropping tuple support, there will be a compiler error for each tuple, and the fix is to change `()` to `[]` wherever there's an error.

There are a few reasons to drop support for tuples:

1. It's less code to support.
2. If we support arbitrary const ndims once Rust has full const generics, we won't be able to support conversions from tuples for all supported ndims (unless Rust also adds some kind of variadic tuple support).
3. IMO, tuples are really intended for heterogeneous elements, so it's a bit awkward to use them as shapes.
4. It would force downstream crates to switch from tuples + arrays to just arrays, which may help slightly with compatibility between downstream crates.
bluss: I'm setting up this PR to be just the non-breaking addition of `.to_shape()`. I've implemented the `into_shape()` changes and prototyped future `ShapeArg` usages for constructors, so I'm ready for the future breaking changes too, to those.

What you said now brings up the future breaking change to `into_shape()`. I would suggest that `.into_shape((2, 3))` changes meaning in the coming breaking change - from default "automatic" layout inference to the new "default row-major"; this has the potential for being a silent breaking change for people, but I don't see a great way of fixing that.
jturner314: > What you said now brings up the future breaking change to `into_shape()`. I would suggest that `.into_shape((2, 3))` changes meaning - from the current "automatic" layout inference to the new "default row-major"; this has the potential for being a silent breaking change for people, but I don't see a great way of fixing that.

Oh, that's a good point. I agree that the order for `into_shape` should change from "automatic" to "default row-major".

We should do what we can to avoid a silent breaking change, because that could cause a real headache for someone using Fortran-layout arrays if they don't notice the relevant item in the release notes. I see a couple of options:

1. For the next breaking release or two, panic or return an error if the user doesn't explicitly specify an order and the result would be different from what it is now. (In other words, check for what would otherwise be a silent change.) After a couple of releases, remove this check. The disadvantage of this option is that the panic/error is at runtime rather than at compile time.

2. For the next breaking release or two, deprecate `.into_shape()` and add a method with a new name. After a couple of releases, rename the new method to `.into_shape()` and deprecate the old name. The disadvantage of this option is that it requires the users to make *two* changes to their code (switch to the new method name for a while, and then switch back to `.into_shape()`).

Hopefully, if we keep the mitigation in place for a couple of breaking releases, we'll catch the vast majority of users.
bluss: It seems there is no perfect solution for `into_shape`, that both has the breaking change but also helps the users.

I'm contemplating now your solution (1) with the addition of an opt-in feature flag - `into_shape_row_major_default` that the user can enable to explicitly enable the new default. Without the feature flag, only explicit order works. They can acknowledge the breaking change by enabling the feature flag, and then we deprecate and remove the feature it in a later release. Convoluted, but maybe that's the best possible solution?
jturner314: The feature flag is a good idea. The primary issue is that if *any* crate in the build enables the feature, *all* crates depending on that version of `ndarray` will be affected. So, if `ndarray-linalg` enables the feature, for example, users of `ndarray-linalg` and `ndarray` won't have to explicitly enable the feature. So, if we go this route, we should strongly recommend that only the final binary crate (i.e. not library crates) enable the feature.
SparrowLii: Just have a question: as for non-standard C or Fortran arrays, can Order become smarter? Look at the following examples:

    use ndarray::{Array, ShapeBuilder};
    fn main() {
        let data = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16];
        let a = Array::from_shape_vec((2, 4, 2).strides((1, 4, 2)), data.to_vec()).unwrap();

        // Unchanged
        let a1 = a.to_shape((2, 4, 2)).unwrap();

        // Add a dimension
        let a2 = a.to_shape((1, 2, 4, 2)).unwrap();

        // Expand a certain dimension
        let a3 = a.to_shape((2, 2, 2, 2)).unwrap();
    }

In these cases, it is completely feasible to output an ArrayView instead of a new Array


SparrowLii: I thought about this problem again, and I think it can be implemented according to the following rules:
Take Order::RawMajor for example, we traversing the incoming and outgoing dimensions from back to front. This matter will encounter the following situations:
1. incoming or outcoming is 1

2. incoming and outgoing are equal

3. incoming is divisible by outgoing
For example, [4] -> [2, 2]. It means that we need to expand the original dimension

4. incoming cannot be divisible by outgoing
For example, [2, 2] -> [4] or [2, 3] -> [3, 2]. 
At this time, we try to merge the last two dimensions of incoming,(If it cannot be merged, it means that reshape cannot be performed). Then we enter the next cycle.
Because the total number of elements of incoming and outgoing is the same, finally incoming can always be divisible by outgoing.
jturner314: I'd suggest a preprocessing step on the input array to merge as many axes as possible toward the "inner" axes (according to the desired order). I think @bluss is adding functionality like this in #979. Then, you don't have to worry about merging axes of the input at the same time as dealing with the output shape/strides.
bluss: Actually as soon as I had mentioned some details I could resist trying to implement it, so I have a work in progress going. i think it's alreadyw working for both C/F order and quite general (can handle things like 3, 4, 5 to 15, 4 and so on. We can compare code when I have time to post. It's quite cool.

Patch:
diff --git a/benches/to_shape.rs b/benches/to_shape.rs
new file mode 100644
index 00000000..a048eb77
--- /dev/null
+++ b/benches/to_shape.rs
@@ -0,0 +1,106 @@
+#![feature(test)]
+
+extern crate test;
+use test::Bencher;
+
+use ndarray::prelude::*;
+use ndarray::Order;
+
+#[bench]
+fn to_shape2_1(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape(4 * 5).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape2_2_same(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((4, 5)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape2_2_flip(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((5, 4)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape2_3(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((2, 5, 2)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_1(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape(3 * 4 * 5).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_2_order(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((12, 5)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_2_outoforder(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((4, 15)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_3c(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape((3, 4, 5)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_3f(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5).f());
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape(((3, 4, 5), Order::F)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_4c(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5));
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape(((2, 3, 2, 5), Order::C)).unwrap()
+    });
+}
+
+#[bench]
+fn to_shape3_4f(bench: &mut Bencher) {
+    let a = Array::<f32, _>::zeros((3, 4, 5).f());
+    let view = a.view();
+    bench.iter(|| {
+        view.to_shape(((2, 3, 2, 5), Order::F)).unwrap()
+    });
+}
diff --git a/src/dimension/mod.rs b/src/dimension/mod.rs
index f4f46e76..4aa7c664 100644
--- a/src/dimension/mod.rs
+++ b/src/dimension/mod.rs
@@ -9,10 +9,10 @@
 use crate::error::{from_kind, ErrorKind, ShapeError};
 use crate::slice::SliceArg;
 use crate::{Ix, Ixs, Slice, SliceInfoElem};
+use crate::shape_builder::Strides;
 use num_integer::div_floor;
 
 pub use self::axes::{Axes, AxisDescription};
-pub(crate) use self::axes::axes_of;
 pub use self::axis::Axis;
 pub use self::broadcast::DimMax;
 pub use self::conversion::IntoDimension;
@@ -23,7 +23,8 @@ pub use self::ndindex::NdIndex;
 pub use self::ops::DimAdd;
 pub use self::remove_axis::RemoveAxis;
 
-use crate::shape_builder::Strides;
+pub(crate) use self::axes::axes_of;
+pub(crate) use self::reshape::reshape_dim;
 
 use std::isize;
 use std::mem;
@@ -40,6 +41,8 @@ mod dynindeximpl;
 mod ndindex;
 mod ops;
 mod remove_axis;
+pub(crate) mod reshape;
+mod sequence;
 
 /// Calculate offset from `Ix` stride converting sign properly
 #[inline(always)]
diff --git a/src/dimension/reshape.rs b/src/dimension/reshape.rs
new file mode 100644
index 00000000..c6e08848
--- /dev/null
+++ b/src/dimension/reshape.rs
@@ -0,0 +1,241 @@
+
+use crate::{Dimension, Order, ShapeError, ErrorKind};
+use crate::dimension::sequence::{Sequence, SequenceMut, Forward, Reverse};
+
+#[inline]
+pub(crate) fn reshape_dim<D, E>(from: &D, strides: &D, to: &E, order: Order)
+    -> Result<E, ShapeError>
+where
+    D: Dimension,
+    E: Dimension,
+{
+    debug_assert_eq!(from.ndim(), strides.ndim());
+    let mut to_strides = E::zeros(to.ndim());
+    match order {
+        Order::RowMajor => {
+            reshape_dim_c(&Forward(from), &Forward(strides),
+                          &Forward(to), Forward(&mut to_strides))?;
+        }
+        Order::ColumnMajor => {
+            reshape_dim_c(&Reverse(from), &Reverse(strides),
+                          &Reverse(to), Reverse(&mut to_strides))?;
+        }
+    }
+    Ok(to_strides)
+}
+
+/// Try to reshape an array with dimensions `from_dim` and strides `from_strides` to the new
+/// dimension `to_dim`, while keeping the same layout of elements in memory. The strides needed
+/// if this is possible are stored into `to_strides`.
+///
+/// This function uses RowMajor index ordering if the inputs are read in the forward direction
+/// (index 0 is axis 0 etc) and ColumnMajor index ordering if the inputs are read in reversed
+/// direction (as made possible with the Sequence trait).
+/// 
+/// Preconditions:
+///
+/// 1. from_dim and to_dim are valid dimensions (product of all non-zero axes
+/// fits in isize::MAX).
+/// 2. from_dim and to_dim are don't have any axes that are zero (that should be handled before
+///    this function).
+/// 3. `to_strides` should be an all-zeros or all-ones dimension of the right dimensionality
+/// (but it will be overwritten after successful exit of this function).
+///
+/// This function returns:
+///
+/// - IncompatibleShape if the two shapes are not of matching number of elements
+/// - IncompatibleLayout if the input shape and stride can not be remapped to the output shape
+///   without moving the array data into a new memory layout.
+/// - Ok if the from dim could be mapped to the new to dim.
+fn reshape_dim_c<D, E, E2>(from_dim: &D, from_strides: &D, to_dim: &E, mut to_strides: E2)
+    -> Result<(), ShapeError>
+where
+    D: Sequence<Output=usize>,
+    E: Sequence<Output=usize>,
+    E2: SequenceMut<Output=usize>,
+{
+    // cursor indexes into the from and to dimensions
+    let mut fi = 0;  // index into `from_dim`
+    let mut ti = 0;  // index into `to_dim`.
+
+    while fi < from_dim.len() && ti < to_dim.len() {
+        let mut fd = from_dim[fi];
+        let mut fs = from_strides[fi] as isize;
+        let mut td = to_dim[ti];
+
+        if fd == td {
+            to_strides[ti] = from_strides[fi];
+            fi += 1;
+            ti += 1;
+            continue
+        }
+
+        if fd == 1 {
+            fi += 1;
+            continue;
+        }
+
+        if td == 1 {
+            to_strides[ti] = 1;
+            ti += 1;
+            continue;
+        }
+
+        if fd == 0 || td == 0 {
+            debug_assert!(false, "zero dim not handled by this function");
+            return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
+        }
+
+        // stride times element count is to be distributed out over a combination of axes.
+        let mut fstride_whole = fs * (fd as isize);
+        let mut fd_product = fd;  // cumulative product of axis lengths in the combination (from)
+        let mut td_product = td;  // cumulative product of axis lengths in the combination (to)
+
+        // The two axis lengths are not a match, so try to combine multiple axes
+        // to get it to match up.
+        while fd_product != td_product {
+            if fd_product < td_product {
+                // Take another axis on the from side
+                fi += 1;
+                if fi >= from_dim.len() {
+                    return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
+                }
+                fd = from_dim[fi];
+                fd_product *= fd;
+                if fd > 1 {
+                    let fs_old = fs;
+                    fs = from_strides[fi] as isize;
+                    // check if this axis and the next are contiguous together
+                    if fs_old != fd as isize * fs {
+                        return Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout));
+                    }
+                }
+            } else {
+                // Take another axis on the `to` side
+                // First assign the stride to the axis we leave behind
+                fstride_whole /= td as isize;
+                to_strides[ti] = fstride_whole as usize;
+                ti += 1;
+                if ti >= to_dim.len() {
+                    return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
+                }
+
+                td = to_dim[ti];
+                td_product *= td;
+            }
+        }
+
+        fstride_whole /= td as isize;
+        to_strides[ti] = fstride_whole as usize;
+
+        fi += 1;
+        ti += 1;
+    }
+
+    // skip past 1-dims at the end
+    while fi < from_dim.len() && from_dim[fi] == 1 {
+        fi += 1;
+    }
+
+    while ti < to_dim.len() && to_dim[ti] == 1 {
+        to_strides[ti] = 1;
+        ti += 1;
+    }
+
+    if fi < from_dim.len() || ti < to_dim.len() {
+        return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
+    }
+
+    Ok(())
+}
+
+#[cfg(feature = "std")]
+#[test]
+fn test_reshape() {
+    use crate::Dim;
+
+    macro_rules! test_reshape {
+        (fail $order:ident from $from:expr, $stride:expr, to $to:expr) => {
+            let res = reshape_dim(&Dim($from), &Dim($stride), &Dim($to), Order::$order);
+            println!("Reshape {:?} {:?} to {:?}, order {:?}\n  => {:?}",
+                     $from, $stride, $to, Order::$order, res);
+            let _res = res.expect_err("Expected failed reshape");
+        };
+        (ok $order:ident from $from:expr, $stride:expr, to $to:expr, $to_stride:expr) => {{
+            let res = reshape_dim(&Dim($from), &Dim($stride), &Dim($to), Order::$order);
+            println!("Reshape {:?} {:?} to {:?}, order {:?}\n  => {:?}",
+                     $from, $stride, $to, Order::$order, res);
+            println!("default stride for from dim: {:?}", Dim($from).default_strides());
+            println!("default stride for to dim: {:?}", Dim($to).default_strides());
+            let res = res.expect("Expected successful reshape");
+            assert_eq!(res, Dim($to_stride), "mismatch in strides");
+        }};
+    }
+
+    test_reshape!(ok C from [1, 2, 3], [6, 3, 1], to [1, 2, 3], [6, 3, 1]);
+    test_reshape!(ok C from [1, 2, 3], [6, 3, 1], to [2, 3], [3, 1]);
+    test_reshape!(ok C from [1, 2, 3], [6, 3, 1], to [6], [1]);
+    test_reshape!(fail C from [1, 2, 3], [6, 3, 1], to [1]);
+    test_reshape!(fail F from [1, 2, 3], [6, 3, 1], to [1]);
+
+    test_reshape!(ok C from [6], [1], to [3, 2], [2, 1]);
+    test_reshape!(ok C from [3, 4, 5], [20, 5, 1], to [4, 15], [15, 1]);
+
+    test_reshape!(ok C from [4, 4, 4], [16, 4, 1], to [16, 4], [4, 1]);
+
+    test_reshape!(ok C from [4, 4], [4, 1], to [2, 2, 4, 1], [8, 4, 1, 1]);
+    test_reshape!(ok C from [4, 4], [4, 1], to [2, 2, 4], [8, 4, 1]);
+    test_reshape!(ok C from [4, 4], [4, 1], to [2, 2, 2, 2], [8, 4, 2, 1]);
+
+    test_reshape!(ok C from [4, 4], [4, 1], to [2, 2, 1, 4], [8, 4, 1, 1]);
+
+    test_reshape!(ok C from [4, 4, 4], [16, 4, 1], to [16, 4], [4, 1]);
+    test_reshape!(ok C from [3, 4, 4], [16, 4, 1], to [3, 16], [16, 1]);
+
+    test_reshape!(ok C from [4, 4], [8, 1], to [2, 2, 2, 2], [16, 8, 2, 1]);
+
+    test_reshape!(fail C from [4, 4], [8, 1], to [2, 1, 4, 2]);
+
+    test_reshape!(ok C from [16], [4], to [2, 2, 4], [32, 16, 4]);
+    test_reshape!(ok C from [16], [-4isize as usize], to [2, 2, 4],
+                  [-32isize as usize, -16isize as usize, -4isize as usize]);
+    test_reshape!(ok F from [16], [4], to [2, 2, 4], [4, 8, 16]);
+    test_reshape!(ok F from [16], [-4isize as usize], to [2, 2, 4],
+                  [-4isize as usize, -8isize as usize, -16isize as usize]);
+
+    test_reshape!(ok C from [3, 4, 5], [20, 5, 1], to [12, 5], [5, 1]);
+    test_reshape!(ok C from [3, 4, 5], [20, 5, 1], to [4, 15], [15, 1]);
+    test_reshape!(fail F from [3, 4, 5], [20, 5, 1], to [4, 15]);
+    test_reshape!(ok C from [3, 4, 5, 7], [140, 35, 7, 1], to [28, 15], [15, 1]);
+
+    // preserve stride if shape matches
+    test_reshape!(ok C from [10], [2], to [10], [2]);
+    test_reshape!(ok F from [10], [2], to [10], [2]);
+    test_reshape!(ok C from [2, 10], [1, 2], to [2, 10], [1, 2]);
+    test_reshape!(ok F from [2, 10], [1, 2], to [2, 10], [1, 2]);
+    test_reshape!(ok C from [3, 4, 5], [20, 5, 1], to [3, 4, 5], [20, 5, 1]);
+    test_reshape!(ok F from [3, 4, 5], [20, 5, 1], to [3, 4, 5], [20, 5, 1]);
+
+    test_reshape!(ok C from [3, 4, 5], [4, 1, 1], to [12, 5], [1, 1]);
+    test_reshape!(ok F from [3, 4, 5], [1, 3, 12], to [12, 5], [1, 12]);
+    test_reshape!(ok F from [3, 4, 5], [1, 3, 1], to [12, 5], [1, 1]);
+
+    // broadcast shapes
+    test_reshape!(ok C from [3, 4, 5, 7], [0, 0, 7, 1], to [12, 35], [0, 1]);
+    test_reshape!(fail C from [3, 4, 5, 7], [0, 0, 7, 1], to [28, 15]);
+
+    // one-filled shapes
+    test_reshape!(ok C from [10], [1], to [1, 10, 1, 1, 1], [1, 1, 1, 1, 1]);
+    test_reshape!(ok F from [10], [1], to [1, 10, 1, 1, 1], [1, 1, 1, 1, 1]);
+    test_reshape!(ok C from [1, 10], [10, 1], to [1, 10, 1, 1, 1], [10, 1, 1, 1, 1]);
+    test_reshape!(ok F from [1, 10], [10, 1], to [1, 10, 1, 1, 1], [10, 1, 1, 1, 1]);
+    test_reshape!(ok C from [1, 10], [1, 1], to [1, 5, 1, 1, 2], [1, 2, 2, 2, 1]);
+    test_reshape!(ok F from [1, 10], [1, 1], to [1, 5, 1, 1, 2], [1, 1, 5, 5, 5]);
+    test_reshape!(ok C from [10, 1, 1, 1, 1], [1, 1, 1, 1, 1], to [10], [1]);
+    test_reshape!(ok F from [10, 1, 1, 1, 1], [1, 1, 1, 1, 1], to [10], [1]);
+    test_reshape!(ok C from [1, 5, 1, 2, 1], [1, 2, 1, 1, 1], to [10], [1]);
+    test_reshape!(fail F from [1, 5, 1, 2, 1], [1, 2, 1, 1, 1], to [10]);
+    test_reshape!(ok F from [1, 5, 1, 2, 1], [1, 1, 1, 5, 1], to [10], [1]);
+    test_reshape!(fail C from [1, 5, 1, 2, 1], [1, 1, 1, 5, 1], to [10]);
+}
+
diff --git a/src/dimension/sequence.rs b/src/dimension/sequence.rs
new file mode 100644
index 00000000..835e00d1
--- /dev/null
+++ b/src/dimension/sequence.rs
@@ -0,0 +1,109 @@
+use std::ops::Index;
+use std::ops::IndexMut;
+
+use crate::dimension::Dimension;
+
+pub(in crate::dimension) struct Forward<D>(pub(crate) D);
+pub(in crate::dimension) struct Reverse<D>(pub(crate) D);
+
+impl<D> Index<usize> for Forward<&D>
+where
+    D: Dimension,
+{
+    type Output = usize;
+
+    #[inline]
+    fn index(&self, index: usize) -> &usize {
+        &self.0[index]
+    }
+}
+
+impl<D> Index<usize> for Forward<&mut D>
+where
+    D: Dimension,
+{
+    type Output = usize;
+
+    #[inline]
+    fn index(&self, index: usize) -> &usize {
+        &self.0[index]
+    }
+}
+
+impl<D> IndexMut<usize> for Forward<&mut D>
+where
+    D: Dimension,
+{
+    #[inline]
+    fn index_mut(&mut self, index: usize) -> &mut usize {
+        &mut self.0[index]
+    }
+}
+
+impl<D> Index<usize> for Reverse<&D>
+where
+    D: Dimension,
+{
+    type Output = usize;
+
+    #[inline]
+    fn index(&self, index: usize) -> &usize {
+        &self.0[self.len() - index - 1]
+    }
+}
+
+impl<D> Index<usize> for Reverse<&mut D>
+where
+    D: Dimension,
+{
+    type Output = usize;
+
+    #[inline]
+    fn index(&self, index: usize) -> &usize {
+        &self.0[self.len() - index - 1]
+    }
+}
+
+impl<D> IndexMut<usize> for Reverse<&mut D>
+where
+    D: Dimension,
+{
+    #[inline]
+    fn index_mut(&mut self, index: usize) -> &mut usize {
+        let len = self.len();
+        &mut self.0[len - index - 1]
+    }
+}
+
+/// Indexable sequence with length
+pub(in crate::dimension) trait Sequence: Index<usize> {
+    fn len(&self) -> usize;
+}
+
+/// Indexable sequence with length (mut)
+pub(in crate::dimension) trait SequenceMut: Sequence + IndexMut<usize> { }
+
+impl<D> Sequence for Forward<&D> where D: Dimension {
+    #[inline]
+    fn len(&self) -> usize { self.0.ndim() }
+}
+
+impl<D> Sequence for Forward<&mut D> where D: Dimension {
+    #[inline]
+    fn len(&self) -> usize { self.0.ndim() }
+}
+
+impl<D> SequenceMut for Forward<&mut D> where D: Dimension { }
+
+impl<D> Sequence for Reverse<&D> where D: Dimension {
+    #[inline]
+    fn len(&self) -> usize { self.0.ndim() }
+}
+
+impl<D> Sequence for Reverse<&mut D> where D: Dimension {
+    #[inline]
+    fn len(&self) -> usize { self.0.ndim() }
+}
+
+impl<D> SequenceMut for Reverse<&mut D> where D: Dimension { }
+
diff --git a/src/impl_methods.rs b/src/impl_methods.rs
index 9ef4277d..6c51b451 100644
--- a/src/impl_methods.rs
+++ b/src/impl_methods.rs
@@ -23,11 +23,14 @@ use crate::dimension::{
     offset_from_ptr_to_memory, size_of_shape_checked, stride_offset, Axes,
 };
 use crate::dimension::broadcast::co_broadcast;
+use crate::dimension::reshape_dim;
 use crate::error::{self, ErrorKind, ShapeError, from_kind};
 use crate::math_cell::MathCell;
 use crate::itertools::zip;
-use crate::zip::{IntoNdProducer, Zip};
 use crate::AxisDescription;
+use crate::order::Order;
+use crate::shape_builder::ShapeArg;
+use crate::zip::{IntoNdProducer, Zip};
 
 use crate::iter::{
     AxisChunksIter, AxisChunksIterMut, AxisIter, AxisIterMut, ExactChunks, ExactChunksMut,
@@ -1577,6 +1580,101 @@ where
         }
     }
 
+    /// Transform the array into `new_shape`; any shape with the same number of elements is
+    /// accepted.
+    ///
+    /// `order` specifies the *logical* order in which the array is to be read and reshaped.
+    /// The array is returned as a `CowArray`; a view if possible, otherwise an owned array.
+    ///
+    /// For example, when starting from the one-dimensional sequence 1 2 3 4 5 6, it would be
+    /// understood as a 2 x 3 array in row major ("C") order this way:
+    ///
+    /// ```text
+    /// 1 2 3
+    /// 4 5 6
+    /// ```
+    ///
+    /// and as 2 x 3 in column major ("F") order this way:
+    ///
+    /// ```text
+    /// 1 3 5
+    /// 2 4 6
+    /// ```
+    ///
+    /// This example should show that any time we "reflow" the elements in the array to a different
+    /// number of rows and columns (or more axes if applicable), it is important to pick an index
+    /// ordering, and that's the reason for the function parameter for `order`.
+    ///
+    /// **Errors** if the new shape doesn't have the same number of elements as the array's current
+    /// shape.
+    ///
+    /// ```
+    /// use ndarray::array;
+    /// use ndarray::Order;
+    ///
+    /// assert!(
+    ///     array![1., 2., 3., 4., 5., 6.].to_shape(((2, 3), Order::RowMajor)).unwrap()
+    ///     == array![[1., 2., 3.],
+    ///               [4., 5., 6.]]
+    /// );
+    ///
+    /// assert!(
+    ///     array![1., 2., 3., 4., 5., 6.].to_shape(((2, 3), Order::ColumnMajor)).unwrap()
+    ///     == array![[1., 3., 5.],
+    ///               [2., 4., 6.]]
+    /// );
+    /// ```
+    pub fn to_shape<E>(&self, new_shape: E) -> Result<CowArray<'_, A, E::Dim>, ShapeError>
+    where
+        E: ShapeArg,
+        A: Clone,
+        S: Data,
+    {
+        let (shape, order) = new_shape.into_shape_and_order();
+        self.to_shape_order(shape, order.unwrap_or(Order::RowMajor))
+    }
+
+    fn to_shape_order<E>(&self, shape: E, order: Order)
+        -> Result<CowArray<'_, A, E>, ShapeError>
+    where
+        E: Dimension,
+        A: Clone,
+        S: Data,
+    {
+        let len = self.dim.size();
+        if size_of_shape_checked(&shape) != Ok(len) {
+            return Err(error::incompatible_shapes(&self.dim, &shape));
+        }
+
+        // Create a view if the length is 0, safe because the array and new shape is empty.
+        if len == 0 {
+            unsafe {
+                return Ok(CowArray::from(ArrayView::from_shape_ptr(shape, self.as_ptr())));
+            }
+        }
+
+        // Try to reshape the array as a view into the existing data
+        match reshape_dim(&self.dim, &self.strides, &shape, order) {
+            Ok(to_strides) => unsafe {
+                return Ok(CowArray::from(ArrayView::new(self.ptr, shape, to_strides)));
+            }
+            Err(err) if err.kind() == ErrorKind::IncompatibleShape => {
+                return Err(error::incompatible_shapes(&self.dim, &shape));
+            }
+            _otherwise => { }
+        }
+
+        // otherwise create a new array and copy the elements
+        unsafe {
+            let (shape, view) = match order {
+                Order::RowMajor => (shape.set_f(false), self.view()),
+                Order::ColumnMajor => (shape.set_f(true), self.t()),
+            };
+            Ok(CowArray::from(Array::from_shape_trusted_iter_unchecked(
+                        shape, view.into_iter(), A::clone)))
+        }
+    }
+
     /// Transform the array into `shape`; any shape with the same number of
     /// elements is accepted, but the source array or view must be in standard
     /// or column-major (Fortran) layout.
diff --git a/src/lib.rs b/src/lib.rs
index a4c59eb6..dfce924e 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -145,6 +145,7 @@ pub use crate::dimension::IxDynImpl;
 pub use crate::dimension::NdIndex;
 pub use crate::error::{ErrorKind, ShapeError};
 pub use crate::indexes::{indices, indices_of};
+pub use crate::order::Order;
 pub use crate::slice::{
     MultiSliceArg, NewAxis, Slice, SliceArg, SliceInfo, SliceInfoElem, SliceNextDim,
 };
@@ -162,7 +163,7 @@ pub use crate::stacking::{concatenate, stack, stack_new_axis};
 
 pub use crate::math_cell::MathCell;
 pub use crate::impl_views::IndexLonger;
-pub use crate::shape_builder::{Shape, ShapeBuilder, StrideShape};
+pub use crate::shape_builder::{Shape, ShapeBuilder, ShapeArg, StrideShape};
 
 #[macro_use]
 mod macro_utils;
@@ -202,6 +203,7 @@ mod linspace;
 mod logspace;
 mod math_cell;
 mod numeric_util;
+mod order;
 mod partial;
 mod shape_builder;
 #[macro_use]
diff --git a/src/order.rs b/src/order.rs
new file mode 100644
index 00000000..e8d9c8db
--- /dev/null
+++ b/src/order.rs
@@ -0,0 +1,83 @@
+
+/// Array order
+///
+/// Order refers to indexing order, or how a linear sequence is translated
+/// into a two-dimensional or multi-dimensional array.
+///
+/// - `RowMajor` means that the index along the row is the most rapidly changing
+/// - `ColumnMajor` means that the index along the column is the most rapidly changing
+///
+/// Given a sequence like: 1, 2, 3, 4, 5, 6
+///
+/// If it is laid it out in a 2 x 3 matrix using row major ordering, it results in:
+///
+/// ```text
+/// 1  2  3
+/// 4  5  6
+/// ```
+///
+/// If it is laid using column major ordering, it results in:
+///
+/// ```text
+/// 1  3  5
+/// 2  4  6
+/// ```
+///
+/// It can be seen as filling in "rows first" or "columns first".
+///
+/// `Order` can be used both to refer to logical ordering as well as memory ordering or memory
+/// layout. The orderings have common short names, also seen in other environments, where
+/// row major is called "C" order (after the C programming language) and column major is called "F"
+/// or "Fortran" order.
+#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+#[non_exhaustive]
+pub enum Order {
+    /// Row major or "C" order
+    RowMajor,
+    /// Column major or "F" order
+    ColumnMajor,
+}
+
+impl Order {
+    /// "C" is an alias for row major ordering
+    pub const C: Order = Order::RowMajor;
+
+    /// "F" (for Fortran) is an alias for column major ordering
+    pub const F: Order = Order::ColumnMajor;
+
+    /// Return true if input is Order::RowMajor, false otherwise
+    #[inline]
+    pub fn is_row_major(self) -> bool {
+        match self {
+            Order::RowMajor => true,
+            Order::ColumnMajor => false,
+        }
+    }
+
+    /// Return true if input is Order::ColumnMajor, false otherwise
+    #[inline]
+    pub fn is_column_major(self) -> bool {
+        !self.is_row_major()
+    }
+
+    /// Return Order::RowMajor if the input is true, Order::ColumnMajor otherwise
+    #[inline]
+    pub fn row_major(row_major: bool) -> Order {
+        if row_major { Order::RowMajor } else { Order::ColumnMajor }
+    }
+
+    /// Return Order::ColumnMajor if the input is true, Order::RowMajor otherwise
+    #[inline]
+    pub fn column_major(column_major: bool) -> Order {
+        Self::row_major(!column_major)
+    }
+
+    /// Return the transpose: row major becomes column major and vice versa.
+    #[inline]
+    pub fn transpose(self) -> Order {
+        match self {
+            Order::RowMajor => Order::ColumnMajor,
+            Order::ColumnMajor => Order::RowMajor,
+        }
+    }
+}
diff --git a/src/shape_builder.rs b/src/shape_builder.rs
index dcfddc1b..47037407 100644
--- a/src/shape_builder.rs
+++ b/src/shape_builder.rs
@@ -1,5 +1,6 @@
 use crate::dimension::IntoDimension;
 use crate::Dimension;
+use crate::order::Order;
 
 /// A contiguous array shape of n dimensions.
 ///
@@ -184,3 +185,33 @@ where
         self.dim.size()
     }
 }
+
+
+/// Array shape argument with optional order parameter
+///
+/// Shape or array dimension argument, with optional [`Order`] parameter.
+///
+/// This is an argument conversion trait that is used to accept an array shape and
+/// (optionally) an ordering argument.
+///
+/// See for example [`.to_shape()`](crate::ArrayBase::to_shape).
+pub trait ShapeArg {
+    type Dim: Dimension;
+    fn into_shape_and_order(self) -> (Self::Dim, Option<Order>);
+}
+
+impl<T> ShapeArg for T where T: IntoDimension {
+    type Dim = T::Dim;
+
+    fn into_shape_and_order(self) -> (Self::Dim, Option<Order>) {
+        (self.into_dimension(), None)
+    }
+}
+
+impl<T> ShapeArg for (T, Order) where T: IntoDimension {
+    type Dim = T::Dim;
+
+    fn into_shape_and_order(self) -> (Self::Dim, Option<Order>) {
+        (self.0.into_dimension(), Some(self.1))
+    }
+}


Test patch:
diff --git a/tests/array.rs b/tests/array.rs
index 976824df..edd58adb 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -8,7 +8,7 @@
 )]
 
 use defmac::defmac;
-use itertools::{enumerate, zip, Itertools};
+use itertools::{zip, Itertools};
 use ndarray::prelude::*;
 use ndarray::{arr3, rcarr2};
 use ndarray::indices;
@@ -1370,64 +1370,6 @@ fn transpose_view_mut() {
     assert_eq!(at, arr2(&[[1, 4], [2, 5], [3, 7]]));
 }
 
-#[test]
-fn reshape() {
-    let data = [1, 2, 3, 4, 5, 6, 7, 8];
-    let v = aview1(&data);
-    let u = v.into_shape((3, 3));
-    assert!(u.is_err());
-    let u = v.into_shape((2, 2, 2));
-    assert!(u.is_ok());
-    let u = u.unwrap();
-    assert_eq!(u.shape(), &[2, 2, 2]);
-    let s = u.into_shape((4, 2)).unwrap();
-    assert_eq!(s.shape(), &[4, 2]);
-    assert_eq!(s, aview2(&[[1, 2], [3, 4], [5, 6], [7, 8]]));
-}
-
-#[test]
-#[should_panic(expected = "IncompatibleShape")]
-fn reshape_error1() {
-    let data = [1, 2, 3, 4, 5, 6, 7, 8];
-    let v = aview1(&data);
-    let _u = v.into_shape((2, 5)).unwrap();
-}
-
-#[test]
-#[should_panic(expected = "IncompatibleLayout")]
-fn reshape_error2() {
-    let data = [1, 2, 3, 4, 5, 6, 7, 8];
-    let v = aview1(&data);
-    let mut u = v.into_shape((2, 2, 2)).unwrap();
-    u.swap_axes(0, 1);
-    let _s = u.into_shape((2, 4)).unwrap();
-}
-
-#[test]
-fn reshape_f() {
-    let mut u = Array::zeros((3, 4).f());
-    for (i, elt) in enumerate(u.as_slice_memory_order_mut().unwrap()) {
-        *elt = i as i32;
-    }
-    let v = u.view();
-    println!("{:?}", v);
-
-    // noop ok
-    let v2 = v.into_shape((3, 4));
-    assert!(v2.is_ok());
-    assert_eq!(v, v2.unwrap());
-
-    let u = v.into_shape((3, 2, 2));
-    assert!(u.is_ok());
-    let u = u.unwrap();
-    println!("{:?}", u);
-    assert_eq!(u.shape(), &[3, 2, 2]);
-    let s = u.into_shape((4, 3)).unwrap();
-    println!("{:?}", s);
-    assert_eq!(s.shape(), &[4, 3]);
-    assert_eq!(s, aview2(&[[0, 4, 8], [1, 5, 9], [2, 6, 10], [3, 7, 11]]));
-}
-
 #[test]
 #[allow(clippy::cognitive_complexity)]
 fn insert_axis() {
diff --git a/tests/reshape.rs b/tests/reshape.rs
new file mode 100644
index 00000000..21fe407e
--- /dev/null
+++ b/tests/reshape.rs
@@ -0,0 +1,232 @@
+use ndarray::prelude::*;
+
+use itertools::enumerate;
+
+use ndarray::Order;
+
+#[test]
+fn reshape() {
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let u = v.into_shape((3, 3));
+    assert!(u.is_err());
+    let u = v.into_shape((2, 2, 2));
+    assert!(u.is_ok());
+    let u = u.unwrap();
+    assert_eq!(u.shape(), &[2, 2, 2]);
+    let s = u.into_shape((4, 2)).unwrap();
+    assert_eq!(s.shape(), &[4, 2]);
+    assert_eq!(s, aview2(&[[1, 2], [3, 4], [5, 6], [7, 8]]));
+}
+
+#[test]
+#[should_panic(expected = "IncompatibleShape")]
+fn reshape_error1() {
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let _u = v.into_shape((2, 5)).unwrap();
+}
+
+#[test]
+#[should_panic(expected = "IncompatibleLayout")]
+fn reshape_error2() {
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let mut u = v.into_shape((2, 2, 2)).unwrap();
+    u.swap_axes(0, 1);
+    let _s = u.into_shape((2, 4)).unwrap();
+}
+
+#[test]
+fn reshape_f() {
+    let mut u = Array::zeros((3, 4).f());
+    for (i, elt) in enumerate(u.as_slice_memory_order_mut().unwrap()) {
+        *elt = i as i32;
+    }
+    let v = u.view();
+    println!("{:?}", v);
+
+    // noop ok
+    let v2 = v.into_shape((3, 4));
+    assert!(v2.is_ok());
+    assert_eq!(v, v2.unwrap());
+
+    let u = v.into_shape((3, 2, 2));
+    assert!(u.is_ok());
+    let u = u.unwrap();
+    println!("{:?}", u);
+    assert_eq!(u.shape(), &[3, 2, 2]);
+    let s = u.into_shape((4, 3)).unwrap();
+    println!("{:?}", s);
+    assert_eq!(s.shape(), &[4, 3]);
+    assert_eq!(s, aview2(&[[0, 4, 8], [1, 5, 9], [2, 6, 10], [3, 7, 11]]));
+}
+
+
+#[test]
+fn to_shape_easy() {
+    // 1D -> C -> C
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let u = v.to_shape(((3, 3), Order::RowMajor));
+    assert!(u.is_err());
+
+    let u = v.to_shape(((2, 2, 2), Order::C));
+    assert!(u.is_ok());
+
+    let u = u.unwrap();
+    assert!(u.is_view());
+    assert_eq!(u.shape(), &[2, 2, 2]);
+    assert_eq!(u, array![[[1, 2], [3, 4]], [[5, 6], [7, 8]]]);
+
+    let s = u.to_shape((4, 2)).unwrap();
+    assert_eq!(s.shape(), &[4, 2]);
+    assert_eq!(s, aview2(&[[1, 2], [3, 4], [5, 6], [7, 8]]));
+
+    // 1D -> F -> F
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let u = v.to_shape(((3, 3), Order::ColumnMajor));
+    assert!(u.is_err());
+
+    let u = v.to_shape(((2, 2, 2), Order::ColumnMajor));
+    assert!(u.is_ok());
+
+    let u = u.unwrap();
+    assert!(u.is_view());
+    assert_eq!(u.shape(), &[2, 2, 2]);
+    assert_eq!(u, array![[[1, 5], [3, 7]], [[2, 6], [4, 8]]]);
+
+    let s = u.to_shape(((4, 2), Order::ColumnMajor)).unwrap();
+    assert_eq!(s.shape(), &[4, 2]);
+    assert_eq!(s, array![[1, 5], [2, 6], [3, 7], [4, 8]]);
+}
+
+#[test]
+fn to_shape_copy() {
+    // 1D -> C -> F
+    let v = ArrayView::from(&[1, 2, 3, 4, 5, 6, 7, 8]);
+    let u = v.to_shape(((4, 2), Order::RowMajor)).unwrap();
+    assert_eq!(u.shape(), &[4, 2]);
+    assert_eq!(u, array![[1, 2], [3, 4], [5, 6], [7, 8]]);
+
+    let u = u.to_shape(((2, 4), Order::ColumnMajor)).unwrap();
+    assert_eq!(u.shape(), &[2, 4]);
+    assert_eq!(u, array![[1, 5, 2, 6], [3, 7, 4, 8]]);
+
+    // 1D -> F -> C
+    let v = ArrayView::from(&[1, 2, 3, 4, 5, 6, 7, 8]);
+    let u = v.to_shape(((4, 2), Order::ColumnMajor)).unwrap();
+    assert_eq!(u.shape(), &[4, 2]);
+    assert_eq!(u, array![[1, 5], [2, 6], [3, 7], [4, 8]]);
+
+    let u = u.to_shape((2, 4)).unwrap();
+    assert_eq!(u.shape(), &[2, 4]);
+    assert_eq!(u, array![[1, 5, 2, 6], [3, 7, 4, 8]]);
+}
+
+#[test]
+fn to_shape_add_axis() {
+    // 1D -> C -> C
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let u = v.to_shape(((4, 2), Order::RowMajor)).unwrap();
+
+    assert!(u.to_shape(((1, 4, 2), Order::RowMajor)).unwrap().is_view());
+    assert!(u.to_shape(((1, 4, 2), Order::ColumnMajor)).unwrap().is_view());
+}
+
+
+#[test]
+fn to_shape_copy_stride() {
+    let v = array![[1, 2, 3, 4], [5, 6, 7, 8]];
+    let vs = v.slice(s![.., ..3]);
+    let lin1 = vs.to_shape(6).unwrap();
+    assert_eq!(lin1, array![1, 2, 3, 5, 6, 7]);
+    assert!(lin1.is_owned());
+
+    let lin2 = vs.to_shape((6, Order::ColumnMajor)).unwrap();
+    assert_eq!(lin2, array![1, 5, 2, 6, 3, 7]);
+    assert!(lin2.is_owned());
+}
+
+
+#[test]
+fn to_shape_zero_len() {
+    let v = array![[1, 2, 3, 4], [5, 6, 7, 8]];
+    let vs = v.slice(s![.., ..0]);
+    let lin1 = vs.to_shape(0).unwrap();
+    assert_eq!(lin1, array![]);
+    assert!(lin1.is_view());
+}
+
+#[test]
+#[should_panic(expected = "IncompatibleShape")]
+fn to_shape_error1() {
+    let data = [1, 2, 3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let _u = v.to_shape((2, 5)).unwrap();
+}
+
+#[test]
+#[should_panic(expected = "IncompatibleShape")]
+fn to_shape_error2() {
+    // overflow
+    let data = [3, 4, 5, 6, 7, 8];
+    let v = aview1(&data);
+    let _u = v.to_shape((2, usize::MAX)).unwrap();
+}
+
+#[test]
+fn to_shape_discontig() {
+    for &create_order in &[Order::C, Order::F] {
+        let a = Array::from_iter(0..64);
+        let mut a1 = a.to_shape(((4, 4, 4), create_order)).unwrap();
+        a1.slice_collapse(s![.., ..;2, ..]); // now shape (4, 2, 4)
+        assert!(a1.as_slice_memory_order().is_none());
+
+        for &order in &[Order::C, Order::F] {
+            let v1 = a1.to_shape(((2, 2, 2, 2, 2), order)).unwrap();
+            assert!(v1.is_view());
+            let v1 = a1.to_shape(((4, 1, 2, 1, 2, 2), order)).unwrap();
+            assert!(v1.is_view());
+            let v1 = a1.to_shape(((4, 2, 4), order)).unwrap();
+            assert!(v1.is_view());
+            let v1 = a1.to_shape(((8, 4), order)).unwrap();
+            assert_eq!(v1.is_view(), order == create_order && create_order == Order::C,
+                       "failed for {:?}, {:?}", create_order, order);
+            let v1 = a1.to_shape(((4, 8), order)).unwrap();
+            assert_eq!(v1.is_view(), order == create_order && create_order == Order::F,
+                       "failed for {:?}, {:?}", create_order, order);
+            let v1 = a1.to_shape((32, order)).unwrap();
+            assert!(!v1.is_view());
+        }
+    }
+}
+
+#[test]
+fn to_shape_broadcast() {
+    for &create_order in &[Order::C, Order::F] {
+        let a = Array::from_iter(0..64);
+        let mut a1 = a.to_shape(((4, 4, 4), create_order)).unwrap();
+        a1.slice_collapse(s![.., ..1, ..]); // now shape (4, 1, 4)
+        let v1 = a1.broadcast((4, 4, 4)).unwrap(); // Now shape (4, 4, 4)
+        assert!(v1.as_slice_memory_order().is_none());
+
+        for &order in &[Order::C, Order::F] {
+            let v2 = v1.to_shape(((2, 2, 2, 2, 2, 2), order)).unwrap();
+            assert_eq!(v2.strides(), match (create_order, order) {
+                (Order::C, Order::C) => { &[32, 16, 0, 0, 2, 1] }
+                (Order::C, Order::F) => { &[16, 32, 0, 0, 1, 2] }
+                (Order::F, Order::C) => { &[2, 1, 0, 0, 32, 16] }
+                (Order::F, Order::F) => { &[1, 2, 0, 0, 16, 32] }
+                _other => unreachable!()
+            });
+
+            let v2 = v1.to_shape(((4, 4, 4), order)).unwrap();
+            assert!(v2.is_view());
+            let v2 = v1.to_shape(((8, 8), order)).unwrap();
+            assert!(v2.is_owned());
+        }
+    }
+}