Reinterpret an Array of Float64 as an Array of SVector{Float64} ?

[LaurentPlagne]

Hi,
I wonder if it possible to reinterpret an array of T as an array of SVector of T ?
something like
a=rand(Float64,4*n)
b=Array{SVector{4,Float64}(Ptr(a),n)

[mschauer]

It possible in Julia 0.6, see eg. https://github.com/mschauer/Bridge.jl/blob/master/src/misc.jl#L77). In Julia 0.7 there will be changes,

• "The return
  type from reinterpret is now ReinterpretArray, which implements
  the AbstractArray interface and does the reinterpreting lazily on
  demand." Implement ReinterpretArray JuliaLang/julia#23750 .

but it should still be possible with an layer of indirection.

[LaurentPlagne]

So you say that:
b=rand(Float64,(4*20))
c=reinterpret(SVector{4,Float64},b)
will not work anymore on 0.7 ?
P.S. I just realized that it worked in 0.6 via your link.

[mschauer]

You can do it on 0.7, it will give a different return type, not an Array{SArray{Tuple{4},Float64,1,4},1} but an

dump(c)
Base.ReinterpretArray{SArray{Tuple{4},Float64,1,4},1,Float64,Array{Float64,1}}
  parent: Array{Float64}((80,)) [0.232836,  … , 0.538685]

which introduces an layer of indirection. You can account for that when writing functions with typed arguments: Functions should expect arguments that are AbstractArrays, not only Arrays.

[andyferris]

I am very interested in the performance impact of this change. I'm building the latest Julia master (new optimizer) now to compare v0.6 to v0.7 performance.

[andyferris]

OK this does indeed seem to be a performance problem!

Current master

   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.7.0-DEV.5152 (2018-05-21 21:19 UTC)
 _/ |\__'_|_|_|\__'_|  |  Commit dc30e38 (0 days old master)
|__/                   |  x86_64-linux-gnu

julia> using StaticArrays

julia> m = rand(3, 1_000_000);

julia> v = reinterpret(SVector{3,Float64}, m, (1_000_000,));
┌ Warning: `reinterpret(::Type{T}, a::Array{S}, dims::NTuple{N, Int}) where {T, S, N}` is deprecated, use `reshape(reinterpret(T, vec(a)), dims)` instead.
│   caller = top-level scope
└ @ Core :0

julia> sum(v);

julia> @time sum(v);
  0.050470 seconds (9 allocations: 320 bytes)

julia> typeof(v)
Base.ReinterpretArray{SArray{Tuple{3},Float64,1,3},1,Float64,Array{Float64,1}}

versus v0.6.2

  _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.2 (2017-12-13 18:08 UTC)
 _/ |\__'_|_|_|\__'_|  |  
|__/                   |  x86_64-linux-gnu

julia> using StaticArrays

julia> m = rand(3, 1_000_000);

julia> v = reinterpret(SVector{3,Float64}, m, (1_000_000,));

julia> sum(v);

julia> @time sum(v);
  0.004232 seconds (84 allocations: 6.123 KiB)

julia> typeof(v)
Array{StaticArrays.SArray{Tuple{3},Float64,1,3},1}

@Keno is this known about?

[andyferris]

See also

JuliaLang/julia#25014

https://discourse.julialang.org/t/big-overhead-with-the-new-lazy-reshape-reinterpret/7635/5

[mschauer]

This is known. Typically it is not so bad, often you can create the array in the form you need it and fill its reinterpreted form,

reinterpret(Float64, vec(v)) .= rand(3*1_000_000)

instead of constructing something and reinterpreting it into the form in which you need it. There was also a discussion here JuliaLang/julia#22849

[andyferris]

My use case is to load (large amounts of) memory-mapped data from e.g. HDF5, as a 3 x N Matrix{Float64}, say, and I have no control of the type unless I fork the HDF5.jl package or something. I'm using memory mapping because I want to avoid copies. The only way I can get the performance I want on v0.7 at the moment appears to be via pointer hacks. 😞

[andyferris]

EDIT: Sorry wrong issue

[Keno]

I'll have a PR for reinterpret shortly.

[mschauer]

Link: JuliaLang/julia#27213

[andyferris]

OK, I think this is resolved satisfactorily for now.

[chethega]

I think this is not resolved yet, please reopen.

We might consider exporting an unsafe_fast_reinterpret based on unsafe_wrap. I don't know whether this would tbaa-unsound, because we effectively only cheat between Ptr{T} and Ptr{SVector{N,T}} (Keno would know?).

It should be possible to make gc work (keeping the buffer or parent alive) with some finalizer-based trickery. But regardless, I would prefer if Base offered unsafe (wrapping-based) reshape and reinterpret variants that deal with all gc issues.

Also see JuliaLang/julia#28980.

On Julia Version 1.0.0 and Julia Version 1.1.0-DEV.265:

julia> using StaticArrays, BenchmarkTools, Random
julia> M = rand(20, 10_000);
julia> M_sv = reshape(reinterpret(SVector{size(M,1), eltype(M)},M), size(M,2));
julia> @btime copy($M);
  184.841 μs (2 allocations: 1.53 MiB)
julia> @btime copy($M_sv);
  1.242 ms (2 allocations: 1.53 MiB)
julia> sv = rand(SVector{20,Float64}, 10_000);
julia> sv_M = reshape(reinterpret(Float64, sv), length(eltype(sv)), length(sv));
julia> @btime copy($sv);
  190.262 μs (2 allocations: 1.53 MiB)
julia> @btime copy($sv_M);
  3.321 ms (2 allocations: 1.53 MiB) #Julia Version 1.0.0
  1.578 ms (2 allocations: 1.53 MiB) #Julia Version 1.1.0-DEV.265.