Why is accessing an SMatrix from a struct allocating?

Considering the following simple MWE:

using StaticArrays, BenchmarkTools
struct SomeStruct
    matrix::SMatrix{3,3,Float64}
end

function testing_allocation(g::MMatrix{3,3,Float64}, struct_matrix::SomeStruct)
  fill!(g, zero(Float64))
  g .+= struct_matrix.matrix
end

@btime testing_allocation(g, struct_matrix) setup = (g = MMatrix{3,3,Float64}(undef); struct_matrix = SomeStruct(rand(SMatrix{3,3,Float64})))

Output:

  458.675 ns (4 allocations: 304 bytes)
3×3 MMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 0.981097  0.90833   0.795176
 0.915934  0.482159  0.0795738
 0.250249  0.906489  0.864033

Removing the intermediate struct of using vectors instead of matrices fixes this issue. Am I missing something obvious?

1 Like

This is an abstract type. The full concrete type is matrix::SMatrix{3,3,Float64,9}.

See also many previous discussions, such as SMatrix{2, 2, Float64} is not a concrete type? and Computing Inverse of a stack of matrices - #19 by tchr

3 Likes

Oh my god, thank you, of course. I assume there is no way to define Some_Struct{n} in such a way that matrix is always a SMatrix{n,n,Float64,n*n}? I just need to add a second type parameter and propagate that up through any structures that include this one?

Unfortunately yes, because Julia’s type system cannot do arithmetic in type parameters.

1 Like

Depending on what you want to do, it might be possible to make the struct a bit more convenient for your purposes. For example, defining

using StaticArrays
struct SomeStruct{N, T<:SMatrix{N,N,Float64}}
        matrix::T
end

allows you to do the following:

a = @SMatrix rand(3,3)
b = @SMatrix rand(4,4)
c = @SMatrix rand(4,3)
sa = SomeStruct(a) # Works fine
sb = SomeStruct(b) # Works fine
sc = SomeStruct(c) # Errors

This may or may not be easier to propagate up your structures depending on your use case.

1 Like