Reinterpret between SVector{N,ComplexF64} and SVector{2N,Float64} without accessing internal fields

General Usage
,
1

Hi, everyone. I want to ask you if you know a better way of doing the following that does not imply accessing internal fields.

Imagine I have a vector of complex numbers:

julia> v1 = [1.0 + 0.5im, 4.8 - 9.0im]
2-element Vector{ComplexF64}:
 1.0 + 0.5im
 4.8 - 9.0im

I can convert it to a vector of floats with:

julia> v2 = reinterpret(Float64, v1)
4-element reinterpret(Float64, ::Vector{ComplexF64}):
  1.0
  0.5
  4.8
 -9.0

And go back likewise:

julia> v3 = reinterpret(ComplexF64, v2)
2-element Vector{ComplexF64}:
 1.0 + 0.5im
 4.8 - 9.0im

Everything works fine as I like it. Now I would like to do the same with SVectors:

julia> s1 = SVector(1.0 + 0.5im, 4.8 - 9.0im)
2-element SVector{2, ComplexF64} with indices SOneTo(2):
 1.0 + 0.5im
 4.8 - 9.0im

The naĂŻve implementation fails:

julia> s2 = reinterpret(SVector{4,Float64}, s1)
Error showing value of type Base.ReinterpretArray{SVector{4, Float64}, 1, ComplexF64, SVector{2, ComplexF64}, false}:

SYSTEM (REPL): showing an error caused an error
ERROR: 1-element ExceptionStack:
DimensionMismatch: 1:1 is inconsistent with SOneTo{2}
Stacktrace:
...

However, if access the internal field data of s1, then it works.

julia> s2 = reinterpret(SVector{4,Float64}, s1.data)
4-element SVector{4, Float64} with indices SOneTo(4):
  1.0
  0.5
  4.8
 -9.0

If I wanted to go back, I could do so accessing again .data:

julia> s3 = reinterpret(SVector{2,ComplexF64}, s2.data)
2-element SVector{2, ComplexF64} with indices SOneTo(2):
 1.0 + 0.5im
 4.8 - 9.0im

Is there a way of doing this without accessing .data? What I show here is MWE of what I need. I cannot afford to allocate intermediate arrays in the process.

Likewise, in another part of the program, I need to convert SMatrix and SHermitianCompact arrays back to the vector that built them, and the only way to do so is accessing those internal fields. I know it is a bad practice in case the api changes in the future, but I haven’t found a better way for now.

2

Just destructure and construct a new ‘SArray’.

If the array is an ‘SArray’, it is supposed not to be heap allocated.

The problem here is that you consider interface violation as “a way”.

1 Like
3

IIRC Tuple(::SArray) is public interface, so s2 = reinterpret(SVector{4,Float64}, Tuple(s1)) is legitimate.

For SMatrix, use vec(matr).

For a Hermitian matrix, you can use StaticArrays.sacollect.

2 Likes
4

I suspect the reason direct reinterpret does not work is that in your s2 = reinterpret(SVector{4,Float64}, s1), s1 <: AbstractArray. So it goes through the ReinterpretArray machinery instead of the path for isbits values. One could try to fix this in StaticArrays.jl, though getting the dispatch right might be difficult.

I think the above suggestion to use Tuple is the right idea here.

5

You’re right, I forgot to call Tuple directly o the SVector. It works like a charm, thanks!

Oh, nice. Even though this method is not exported, it’s very useful. I’ve come up with something like this:

import StaticArrays: sacollect

function f(::Val{N}, A) where N
    L = div(N*(N+1), 2)
    return sacollect(SVector{L, eltype(A)}, A[i,j] for j in 1:N for i in j:N)
end

It is fast, it doesn’t allocate and does not depend on internal fields.

1 Like
6

Yes, that’s what I was thinking at the beginning, but I couldn’t figure out the way of doing it withoug calling to those internal fields. The solution proposed by @Vasily_Pisarev avoids this problem, too.

What I meant is that I could always transform the SArray to Array, do the reinterpret and go back to SArray, all without internal fields, but at the cost of allocating.

1 Like
7

I realized there’s an even simpler method by reinterpreting Scalar{<:SArray} which works with the same syntax for all of the cases:

julia> using StaticArrays

julia> s_2cf64 = SVector(1.0 + 0.5im, 4.8 - 9.0im);

julia> reinterpret(SVector{4,Float64}, Scalar(s_2cf64))[]
4-element SVector{4, Float64} with indices SOneTo(4):
  1.0
  0.5
  4.8
 -9.0

julia> s_2x2fp64 = @SMatrix [1.0 2.0; 3.0 4.0];

julia> reinterpret(SVector{4,Float64}, Scalar(s_2x2fp64))[]
4-element SVector{4, Float64} with indices SOneTo(4):
 1.0
 3.0
 2.0
 4.0

julia> s_2x2hf64 = SHermitianCompact{2,Float64}([1, 2, 3]);

julia> reinterpret(SVector{3,Float64}, Scalar(s_2x2hf64))[]
3-element SVector{3, Float64} with indices SOneTo(3):
 1.0
 2.0
 3.0
1 Like
8

Wow, that’s cool! I guess in the case of the SHermitianCompact what is happening with Scalar is that it is treating everything as a block and then reinterpret can read the .lowertriangle data? I don’t completely understand it, though. It even allows me to directly reinterpret a complex hermitian matrix into a vector of floats:

julia> sh = SHermitianCompact{2, ComplexF64, 3}([1.3 - 6.8im, 0.5 + 2im, -7.3 + 1.4im])
2×2 SHermitianCompact{2, ComplexF64, 3} with indices SOneTo(2)×SOneTo(2):
 1.3-6.8im   0.5-2.0im
 0.5+2.0im  -7.3+1.4im

julia> reinterpret(SVector{6,Float64}, Scalar(sh))[]
6-element SVector{6, Float64} with indices SOneTo(6):
  1.3
 -6.8
  0.5
  2.0
 -7.3
  1.4

I’ve tested both approachs with @btime and for 3x3 matrices they seem equal.

9

Yes. Scalar is just an alias for 0-rank SArray. Basically, reinterpret(S, ::SArray{N,T}) has issues whenever sizes of S and T don’t match, but we can add another layer of container on top so that the outer container does not change its size. The trick here is to use a static array as the outer container to avoid heap allocations.

1 Like