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.

5

Is it possible to abuse Array constructor from pointer to make a view as separate array?
IIRC, the docs say to not do that, but are there real implications if the view has a limited scope and the parent (or maybe “donor” would be more fitting for such a case) is not resized?

6

Yes:

julia> x3 = unsafe_wrap(Array, pointer(v2), (4, 5));

julia> x3 == x2
true

julia> x3 isa StridedMatrix
true
7

You can also use Base.wrap to wrap the memory of some vector:

julia> v = float.(1:8);

julia> Base.wrap(Array, v.ref, (2, 4))
2×4 Matrix{Float64}:
 1.0  3.0  5.0  7.0
 2.0  4.0  6.0  8.0
8

Yes, using unsafe_wrap to make two Memory with different eltypes alias is very much undefined behavior, so IIUC julia may give LLVM nonsense input. If LLVM gets nonsense input it is hard to predict exactly what could go wrong.

9

But also it doesn’t have to be that way. Maybe in a future version of Julia using unsafe_wrap like that is fine.

10

This is not two arbitrary eltypes, however, but rather Float64 and ComplexF64, and the relationship between their memory layouts is known and should not change (many, many things rely on this; that’s why e.g. C and C++ included it in their standards).

11

I think @nhz2 is referring to how it can break TBAA (type based alias analysis), not the potential for layout differences.

unsafe_load and unsafe_store! disable TBAA in the function body they’re used in to defend against the sort of problems that can be caused here, but I’m not sure if unsafe_wrap would do the same or not.

Edit: from the docstring of unsafe_wrap:

Unlike unsafe_load and unsafe_store!, the programmer is responsible also for ensuring that the underlying data is not accessed through two arrays of different element type, similar to the strict aliasing rule in C.

So if you want to use the pointer to convert between arrays of Floats / ComplexFloats, you’ll need to use a special array type that unsafe_loads / unsafe_stores.

12

Perhaps an alternative solution is if SciML/OrdinaryDiffEq can co-evolve multiple quantities? i.e. if f!(du, u, p, t) can be implemented for u::Tuple of concrete array types? I’m not sure how then to feed this structure to the integrator in OrdinaryDiffEq.solve.