# Convert Vector{SVector} to Matrix

I have some code where I need to store coordinates `x` of a mesh. I’m using a vector of points (static vectors), i.e.

``````using StaticArrays
Point{dim} = SVector{dim, Float64}
x = rand(Point{3}, 10)
``````

but some parts of the code (e.g. using WriteVTK) require the coordinates to be in a matrix instead of a vector. Is there anything wrong with adding the following constructor method to convert the above vector to a matrix

``````Base.Matrix(v::Vector{SVector{S, T}}) where {S, T} = reduce(hcat, v)
``````

and using

``````m = Matrix(x)
mT = transpose(Matrix(x))
``````

to get the vector as a matrix and matrix transpose? (EDIT: Yes, because it creates a copy it is not efficient)

Is there a more efficient way to do this? (EDIT: Yes, use `reinterpret` and `reshape` as explained below)

Is there any way to do this without making a copy so that modifying elements of matrix `m` will change the corresponding elements of `x`? (EDIT: Yes, as above)

Checkout the amazing `reinterpret`:

``````julia> reinterpret(Float64, x)
30-element reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}):
0.8552062516276724
0.6691039952034261
...
``````

which you can then `reshape`. That could maybe be somewhere in the staticarrays doc.

7 Likes

Possibly even a `convert` method? So an array of svectors would be automatically converted to a matrix. But possibly that’d be a bit confusing…

Thanks @antoine-levitt, I wasn’t aware of the `reinterpret` function. This does appear to be much faster for large arrays:

``````using StaticArrays

Base.Matrix(v::Vector{SVector{S, T}}) where {S, T} = reduce(hcat, v)

eltocols(v::Vector{SVector{dim, T}}) where {dim, T} =
reshape(reinterpret(Float64, v), dim, :)

eltorows(v::Vector{SVector{dim, T}}) where {dim, T} =
transpose(reshape(reinterpret(Float64, v), dim, :))

Point{dim} = SVector{dim, Float64}

x = rand(Point{3}, 1000000)

using BenchmarkTools

@benchmark m = Matrix(\$x)
@benchmark mT = transpose(Matrix(\$x))
@benchmark m = eltocols(\$x)
@benchmark mT = eltorows(\$x)
``````

Results:

``````julia> @benchmark m = Matrix(\$x)
BenchmarkTools.Trial:
memory estimate:  22.89 MiB
allocs estimate:  2
--------------
minimum time:     10.874 ms (0.00% GC)
median time:      17.070 ms (0.00% GC)
mean time:        15.128 ms (4.16% GC)
maximum time:     53.498 ms (78.78% GC)
--------------
samples:          331
evals/sample:     1

julia> @benchmark mT = transpose(Matrix(\$x))
BenchmarkTools.Trial:
memory estimate:  22.89 MiB
allocs estimate:  3
--------------
minimum time:     10.695 ms (0.00% GC)
median time:      11.154 ms (0.00% GC)
mean time:        11.346 ms (2.48% GC)
maximum time:     52.555 ms (78.71% GC)
--------------
samples:          441
evals/sample:     1

julia> @benchmark m = eltocols(\$x)
BenchmarkTools.Trial:
memory estimate:  64 bytes
allocs estimate:  2
--------------
minimum time:     13.549 ns (0.00% GC)
median time:      15.117 ns (0.00% GC)
mean time:        17.357 ns (6.48% GC)
maximum time:     1.122 μs (97.98% GC)
--------------
samples:          10000
evals/sample:     998

julia> @benchmark mT = eltorows(\$x)
BenchmarkTools.Trial:
memory estimate:  80 bytes
allocs estimate:  3
--------------
minimum time:     16.939 ns (0.00% GC)
median time:      17.799 ns (0.00% GC)
mean time:        20.796 ns (8.00% GC)
maximum time:     1.255 μs (98.28% GC)
--------------
samples:          10000
evals/sample:     998
``````

Moreover, it doesn’t create a copy so the original `x` vector can be modified by changing the entries of the matrix `m`.

``````julia> x = rand(Point{3}, 1000000)
1000000-element Array{SArray{Tuple{3},Float64,1,3},1}:
[0.4185398323656475, 0.501898669373563, 0.41183968543604865]

julia> m = eltocols(x)
3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64:
0.41854   0.86096   0.762052  0.801235  …  0.293778  0.1493     0.462913
0.501899  0.964937  0.796707  0.379564     0.932526  0.0662906  0.149578
0.41184   0.790648  0.52458   0.837944     0.921923  0.477261   0.267126

julia> m[1,1] = 0

julia> x
1000000-element Array{SArray{Tuple{3},Float64,1,3},1}:
[0.0, 0.501898669373563, 0.41183968543604865]
``````
1 Like
4 Likes

The doc explains how to convert from a dynamically sized `AbstractArray` to one of the statically sized array types, and how to `reinterpret` a `dim` x `N` `Matrix{T}` as a `Vector{SVector{dim,T}}` but not how to `reinterpret` a `Vector{SVector{dim,T}}` as a `Matrix{T}`.

However, it is still very helpful.

Thank you so much!

`reinterpret` does not give you a `Matrix` though.

Thanks @kristoffer.carlsson, I didn’t think of that. In practice they appear to behave the same (see below) which is what I need. What should I call this “matrix-like” type: `3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64`? Could it be a an `AbstractMatrix`? And is there any difference in practice apart from the different type (for dispatch etc.)?

Comparing `Matrix(x)` with `eltocols(x)`:

``````julia> m = Matrix(x)
3×1000000 Array{Float64,2}:
0.468326  0.568942  0.599504  0.569535  …  0.917145   0.440911  0.739385
0.470099  0.442289  0.811437  0.809661     0.559429   0.42691   0.382616
0.200765  0.235891  0.985741  0.960975     0.0103481  0.41207   0.507787

julia> r = eltocols(x)
3×1000000 reshape(reinterpret(Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}), 3, 1000000) with eltype Float64:
0.468326  0.568942  0.599504  0.569535  …  0.917145   0.440911  0.739385
0.470099  0.442289  0.811437  0.809661     0.559429   0.42691   0.382616
0.200765  0.235891  0.985741  0.960975     0.0103481  0.41207   0.507787

julia> using LinearAlgebra

julia> norm(m - r)
0.0

julia> norm(m[1,1:10] - r[1,1:10])
0.0``````
1 Like

I guess I just answered my own question:

``````julia> typeof(m)
Array{Float64,2}

julia> supertype(typeof(m))
DenseArray{Float64,2}

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

julia> supertype(typeof(r))
AbstractArray{Float64,2}

julia> typeof(r) <: AbstractMatrix
true
``````

Thanks! I was looking for a method to convert a `Vector` of `SVector`s into a matrix without allocation, and your technique is the solution.

But as you mentioned, I notice that modifying the resulting `ReshapedArray` entries modifies the entries of the original `SVector`s. How is this possible? I thought the entries of `SVector`s are unmutable. Does this mean that we can always bypass the unmutability of an unmutable type (like `SVector` in the present example) by `reinterpret`ing it as a mutable type (like `ReshapedArray` in the present example)?

2 Likes

Note you can now also `reinterpret(reshape, ...`:

``````julia> x = [@SVector(rand(3)) for _ in 1:5];

julia> reinterpret(reshape, eltype(eltype(x)), x)
3×5 reinterpret(reshape, Float64, ::Vector{SVector{3, Float64}}) with eltype Float64:
0.12019   0.0829133  0.580475  0.873409  0.710157
0.97712   0.556687   0.696489  0.419789  0.0425647
0.506971  0.106399   0.887721  0.116239  0.486331
``````
1 Like

They’re elements of a mutable array.
You can always edit the mutable array `x` in my above example, by replacing one immutable `SVector` with another immutable `SVector`.
The `reinterpret(reshape` array does alias `x`, so mutating it will also change the contents of the mutable `x`.

1 Like

For example, consider

``````julia> a = [SVector(1,2,3), SVector(4,5,6)]
2-element Vector{SVector{3, Int64}}:
[1, 2, 3]
[4, 5, 6]

julia> A = reinterpret(reshape, Int64, a)
3×2 reinterpret(reshape, Int64, ::Vector{SVector{3, Int64}}) with eltype Int64:
1  4
2  5
3  6
``````

Here, if I set `A[1,1] = 0`, then `a` also changes:

``````julia> A[1,1] = 0
0

julia> a
2-element Vector{SVector{3, Int64}}:
[0, 2, 3]
[4, 5, 6]
``````

To me, replacing `A[1,1]` seems equivalent to replacing the first entry to `a[1]`. Because `a[1]` is `SVector`, being able to replace its first entry seems strange to me.

But you say that replacing `A[1,1]` is actually equivalent to replacing the entire vector `a[1]`. I am not sure how replacing one element `A[1,1]` could be interpreted as replacing the entire `a[1]` with three entries. Still, I think you are correct, because if we make `a` an `SVector` of `SVector`s, we cannot replace `A[1,1]`:

``````julia> a = SVector(SVector(1,2,3), SVector(4,5,6))
2-element SVector{2, SVector{3, Int64}} with indices SOneTo(2):
[1, 2, 3]
[4, 5, 6]

julia> A = reinterpret(reshape, Int64, a)
3×2 reinterpret(reshape, Int64, ::SVector{2, SVector{3, Int64}}) with eltype Int64 with indices Base.OneTo(3)×SOneTo(2):
1  4
2  5
3  6

julia> A[1,1] = 0
ERROR: setindex!(::SVector{2, SVector{3, Int64}}, value, ::Int) is not defined.
Hint: Use `MArray` or `SizedArray` to create a mutable static array
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] setindex!(a::SVector{2, SVector{3, Int64}}, value::SVector{3, Int64}, i::Int64)
@ StaticArrays ~/Programming/StaticArrays/src/indexing.jl:3
[3] _setindex_ra!
@ ./reinterpretarray.jl:648 [inlined]
[4] setindex!(::Base.ReinterpretArray{Int64, 2, SVector{3, Int64}, SVector{2, SVector{3, Int64}}, true}, ::Int64, ::Int64, ::Int64)
@ Base ./reinterpretarray.jl:507
[5] top-level scope
@ REPL[33]:1
``````
1 Like

I think the `SVector` here is distracting you. We can do this even just with `Int`s:

``````julia> b = [1, 3];

julia> B = reinterpret(Bool, b)
16-element reinterpret(Bool, ::Vector{Int64}):
1
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
``````
``````julia> B[1] = false; B[2] = true;

julia> b
2-element Vector{Int64}:
256
3
``````

Just like an `Int`, an `SVector{Int}` is just a collection of bits. Once we put those bits in an array, the memory of the array is mutable so we can change their values.

1 Like

You are not mutating an immutable data structure, if that’s what you are worrying about:

``````as = SVector(1,2,3)
a = [as, SVector(4,5,6)]
A = reinterpret(reshape, Int64, a)
A[1,1] = 0
``````

Now,

``````julia> a[1]
3-element SVector{3, Int64} with indices SOneTo(3):
0
2
3
``````

as you know, but `as` is unchanged.

``````julia> as
3-element SVector{3, Int64} with indices SOneTo(3):
1
2
3
``````

You have only mutated the mutable variable `a`.

@Mason, this doesn’t explain why we cannot mutate `A[1,1]` reinterpreted from `a` that is `SVector` of `SVector`s as shown in the last example in my earlier posting above.

In your final example, you have an `SVector` of `SVector`s as the base memory. That’s all immutable so of course you can’t mutate it, unlike the case where you have an mutable array of `SVector`s

Sure, I’m just just pointing out that my final example contradicts your following claim:

Actually I had the same thought as your above claim, but I was puzzled because my final example didn’t allow me to change the bit array.

No it doesnt. Note that I wrote Once we put those bits in an array, the memory of the array is mutable so we can change their values.

SArrays are not mutable. Your example would work if you used an MArray.

I feel that the discussion is going a bit tangential, but what did you mean by this?

I thought it meant performing `A = reinterpret(reshape, Int64, a)`. Then it sounded like we could change the contents of `A` regardless of the type of `a`, because this operation put `a` into an array `A` of bits.