Convert Vector{SVector} to Matrix

General Usage Performance
#1

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)

#2

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.

4 Likes
#3

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…

#4

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]
#5

It is
https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api.html#Arrays-of-static-arrays-1

2 Likes
#6

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!

#7

reinterpret does not give you a Matrix though.

#8

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
#9

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