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.

1 Like
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.

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.