Struggling with type stability - is SMatrix type stable?

General Usage Performance
,
1

I have a problem where I want to collect an Array of vectors as I perform an iteration. Something like this:

function foo(x::Vector{Float64}, size::Int8)
  tmat = SMatrix{4,4}{Float64}(eye(4))
  tmatrices = Array{SMatrix{4,4}{Float64}}(size)
  for j=1:size
	tmatrices[j] = bar(x[j])
	tmat *= tmatrices[j]
  end
  return tmat, tmatrices;
end

function bar(x)
    return eye(4)
end

where SMatrix is a Static Arrays type. When I run this using @code_warntype, Julia complains that tmat is of an abstract type:

tmat::Any = SSAValue(4) # line 5:
  tmatrices::Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1} = $(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}, svec(Any, Int64), Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}, 0, 5, 0)) # line 6:
  SSAValue(6) = (Base.select_value)((Base.sle_int)(1, size::Int64)::Bool, size::Int64, (Base.sub_int)(1, 1)::Int64)::Int64
  #temp#@_5::Int64 = 1
  35: 
  unless (Base.not_int)((#temp#@_5::Int64 === (Base.add_int)(SSAValue(6), 1)::Int64)::Bool)::Bool goto 48
  SSAValue(7) = #temp#@_5::Int64
  SSAValue(8) = (Base.add_int)(#temp#@_5::Int64, 1)::Int64
  j::Int64 = SSAValue(7)
  #temp#@_5::Int64 = SSAValue(8) # line 7:
  SSAValue(2) = $(Expr(:invoke, MethodInstance for eye(::Type{Float64}, ::Int64, ::Int64), :(Base.eye), Float64, 4, 4))
  $(Expr(:invoke, MethodInstance for setindex!(::Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}, ::Array{Float64,2}, ::Int64), :(Main.setindex!), :(tmatrices), SSAValue(2), :(j))) # line 8:
  tmat::Any = (tmat::Any * (Base.arrayref)(tmatrices::Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}, j::Int64)::StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L)::Any
  46: 
  goto 35
  48:  # line 10:
  return (Core.tuple)(tmat::Any, tmatrices::Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1})::Tuple{Any,Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}}
end::Tuple{Any,Array{StaticArrays.SArray{Tuple{4,4},Float64,2,L} where L,1}}

Why?

2

tmat starts out as an SArray, but multiplying by an Array changes its type to the latter, hence the warning.

Not sure what your code does, but consider I instead of eye.

3

I’ve modified the code thus:

using StaticArrays

function foo(x::Vector{Float64}, size)
  tmat::SMatrix{4,4}{Float64} = SMatrix{4,4}{Float64}(@SMatrix eye(4))
  tmatrices = Array{SMatrix{4,4}{Float64}}(size)
  for j=1:size
 	tmatrices[j] = bar(x[j])
tmat *= tmatrices[j]
  end
  return tmat, tmatrices;
end

function bar(x)
    return @SMatrix eye(4)
end

but there is no change.

4

Changing it to:

function bar(x)
    return I
end

…also makes no difference:

tmat::Any = SSAValue(7) # line 5:
5

Yes, SMatrix{4,4,Float64} is an abstract type (alias).

julia> x = @SMatrix randn(4,4)
4Ă—4 SArray{Tuple{4,4},Float64,2,16}:
 -1.55576    0.571037   0.927532   1.27703 
  0.162367  -1.78089   -0.832958   0.760511
  0.417475  -0.622358  -1.68881   -1.12098 
 -0.973613   0.941918   0.337696   1.63527 

julia> typeof(x)
SArray{Tuple{4,4},Float64,2,16}

For your array, use SArray{Tuple{4,4},Float64,2,16}, SMatrix{4,4,Float64,16}, or typeof(tmat).
That is use, Array{SArray{Tuple{4,4},Float64,2,16}}(undef, size) or Array{typeof(tmat)}(undef, size).

EDIT:
Your array as written would accept any of the following subtypes of SMatrix{4,4,Float64}:

julia> SMatrix{4,4,Float64,10} <: SMatrix{4,4,Float64}
true

julia> SMatrix{4,4,Float64,100} <: SMatrix{4,4,Float64}
true

julia> SMatrix{4,4,Float64,16} <: SMatrix{4,4,Float64}
true

However, all but the last are impossible, so you’d only want the last anyway.

1 Like
6

OK, attempting to convert to I from eye(). There are places in my code where I really need to store an “initial” matrix that is the identity matrix. Like so:

struct M
    transform::Array{Float64}
end

function foo()
    return M(I)
end

which gives the error:

ERROR: MethodError: Cannot `convert` an object of type UniformScaling{Int64} to an object of type Array{Float64,N} where N

I seems to work fine “inline” when doing matrix algebra but it can’t be stored in a place where you want an array.

1 Like
7

Your array as written would accept any of the following subtypes of SMatrix{4,4,Float64} :

julia> SMatrix{4,4,Float64,10} <: SMatrix{4,4,Float64}
true

julia> SMatrix{4,4,Float64,100} <: SMatrix{4,4,Float64}
true

julia> SMatrix{4,4,Float64,16} <: SMatrix{4,4,Float64}
true

…huh. Why aren’t the “{4,4}” and the “16” parts of that declaration redundant?

8

Yeah, that is unfortunate, because we can derive the “16” deterministic ally from the “4”.
We do need a 16, because the static arrays wrap NTuples.

I think there is an issue meant to allow that feature by allowing some calculations in struct declarations, so we could have something like:

struct SArray{S,T,N} <: AbstractArray{T,N}
    data::NTuple{prod(S),T}
end

but I couldn’t find it. Still, I think this will get resolved someday! Anyone have a link to that issue?

As an alternative, maybe because SMatrices are so heavily used, we could define something like:

struct SMatrix{M,N,T} <: AbstractMatrix{T}
    data::NTuple{N,NTuple{M,T}}
end

to be slightly more convenient on the user’s side of things.

EDIT:
To get an array from I,

using LinearAlgebra
Matrix{Float64}(I,m,n)
2 Likes
9

Not sure what you’re trying to do, but I have some remarks.

This line is messy:

tmat::SMatrix{4,4}{Float64} = SMatrix{4,4}{Float64}(@SMatrix eye(4))

Remove the type assertion, it only causes problems. Also, you are putting an SMatrix in an SMatrix constructor, which is redundant, and may cause more trouble.

Finally, the size input argument in your function is not a good idea, size is the name of a very important function, and also I suspect that it is unneeded. This code is type-stable, I think:

function foo(x::AbstractVector)
    tmat = @SMatrix eye(4)
    tmatrices = Vector{typeof(tmat)}(length(x))
    for j in eachindex(x)
        tmatrices[j] = bar(x[j])
        tmat *= tmatrices[j]
    end
    return tmat, tmatrices
end

function bar(x)
    return @SMatrix eye(4)
end

In general, try to remove type annotations, and see if that helps – some times they just get in your way.

1 Like
10

OK, in my real code I have created a type alias:

    const TransformationMatrix = SMatrix{4,4,Float64,16}

and I have replaced all instances of eye(4) with either I or

one(TransformationMatrix)

and with a bit of other fiddling I do not have any @code_warntype issues, and my code is now 2x faster. Thanks!

This line is messy:

Yeah, I have some cleanup to do now.

1 Like