Bug in `strides` of reshaped adjoint/transposed arrays

[ranocha]

Reported by @chriselrod in https://github.com/JuliaArrays/ArrayInterface.jl/pull/161#discussion_r643907939:

julia> using ArrayInterface

julia> A = rand(5,12);

julia> B = reshape(A', (3,4,5));

julia> ArrayInterface.strides(B) # @chriselrod: should be `(5,15,StaticInt{1}())`
(static(1), 3, 12)

[ranocha]

Another example (continuing the one above):

julia> ArrayInterface.strides(vec(A'))
(static(1),)

This doesn't make sense at all since vec(A') isn't a plain strided vector.

[Tokazama]

The first example here is fixed. Did we want an error on the second one? This is what we get now:

julia> ArrayInterface.strides(vec(A'))
ERROR: ArgumentError: Input is not strided.

[ranocha]

Did we want an error on the second one?

Yes, this sounds reasonable to me since the input is indeed not a plain strided vector.

The first example here is fixed.

I get

(@v1.7) pkg> activate --temp
  Activating new project at `/tmp/jl_4vH1Gz`

(jl_4vH1Gz) pkg> add ArrayInterface
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
    Updating `/tmp/jl_4vH1Gz/Project.toml`
  [4fba245c] + ArrayInterface v6.0.14

julia> using ArrayInterface

julia> A = rand(5,12);

julia> B = reshape(A', (3,4,5));

julia> ArrayInterface.strides(B) # @chriselrod: should be `(5,15,StaticInt{1}())`
(5, 15, 1)

When reporting this, @chriselrod mentioned it should be (5, 15, StaticInt{1}()). What do you think about this last part - static or not?

[N5N3]

It should not be static.

julia> A = rand(5,12);

julia> B = reshape(A', (3,4,5));

julia> C = randn(1,60);

julia> D = reshape(C', (3,4,5));

julia> typeof(B) == typeof(D)
true

julia> using ArrayInterface

julia> ArrayInterface.strides(B)
(5, 15, 1)

julia> ArrayInterface.strides(D)
(1, 3, 12)

[chriselrod]

Both 1s above should be static, otherwise performance with LoopVectorization.jl will be garbage.

static(1) is needed for avoiding the extremely slow gather and scatter instructions.

[N5N3]

I know non-static result might cause performance issue. But, the above example shows that we can't derive the dense dimension with the current type information.

One possible solution is adding some SizedArray wrappers, and make strides(::SizedArray{S, T, N, <:Base.ReshapedArray{T,N}}) static if the parent also have static size.

[chriselrod]

Why is D's first stride 1? I'll have to look at this later.

Does base have special casing based on array size changing how cartesian indices are mapped to linear? That's rather bizarre behavior.

[chriselrod]

Interesting, those strides are indeed correct.

So, when you reshape an axis, the reshape is in column-major order. Hence

julia> C = randn(1,60);

julia> D = reshape(C', (3,4,5));

D is basically column-major. While A isn't

julia> A = rand(5,12);

julia> B = reshape(A', (3,4,5));

Maybe LV should just do a dynamic dispatch (or, hopefully a union-split) in a case like this. A dynamic dispatch would be cheaper than gather/scatters, unless the arrays are small.