Compiler mis-recognizes contiguous reinterpreted matrices?

General Usage
, , ,
1

I am implementing an ODE to use with OrdinaryDiffEq. I am simultaneously evolving in time multiple quantities, which are combinations of vectors and matrices, some real and some complex. The approach I took has been to group them all in a long Vector{Float64} and and my function f!(du, u, p ,t) which computes derivates, I slice into that long vector. However, I noticed that this leads to poor matrix-matrix multiplication performance on the obtained slices for some methods of slicing. I tracked down the issue to the fact gemm is not called for certain matrix types, though it seems that it should be possible. A MWE is the following:

v1 = rand(Float64, 100);
v2 = reinterpret(ComplexF64, v1);
x1 = view(v2, 1:20);
x2 = reshape(view(v2, 1:20), 4, 5);
isa.((x1, x2), StridedVecOrMat)
# (true, false)

This suggests that though x2 and x1 refer to the same contiguous memory, the compiler cannot reason that x2 is such, and thus calls generic matmul methods instead of gemm. Can this be considered a bug? Is there a better way to annotate the types?

(I do have a temporary workaround: first slice into the float vector, then reinterpret. Is there a rule of thumb to know when these operations do not "commute?)

2

This could be fixed by Julia having strided array traits:

julia> Base.is_contiguous(typeof(x2))
true

julia> Base.is_contiguous(typeof(x1))
true
3

It’s not a limitation of the compiler, it’s a limitation of how the StridedArray type is defined here.

x1 is a StridedVector because it is a StridedSubArray, which includes a SubArray of a StridedReinterpretArray such as v2.

However, the ReshapedArray, x2, is not a StridedReshapedArray because that definition does not include reshaped arrays made from a SubArray of a reinterpreted array, only subarrays of DenseArray.

The basic problem here is that these types are defined as a hierarchy and couldn’t be mutually recursive (until recently, at least?), which is what you would really want here for a StridedArray type (e.g. a StridedArray should include any strided subarray of any StridedArray). Perhaps this could be fixed by the new typegroup facility in Julia 1.14?

A trait would give more flexibility, I guess.

4

Yes, the main flexibility a trait gives is for array types not part of Base. Since the StridedArray type is defined in Base it can only reference other types defined in Base.