`rand(SMatrix{N,N}, M)` produces array of abstractly typed values

[c42f]

Originally spotted in the wild on discourse https://discourse.julialang.org/t/grassmann-jl-a-b-3x-faster-than-julias-staticarrays-jl/41451/1

The former of these two cases you'd naturally expect to work, but it allocates an array of pointers rather than allocating values inline:

julia> typeof(rand(SMatrix{5,5},10000))
Array{SArray{Tuple{5,5},T,2,L} where L where T,1}

julia> typeof(rand(SMatrix{5,5,Float64,25},10000))
Array{SArray{Tuple{5,5},Float64,2,25},1}

I guess we can try to overload rand for a set of our abstract types in an attempt to patch this up. The Random stdlib doesn't provide a way to make the output array concrete, it appears we dispatch to the following

rand(r::AbstractRNG, ::Type{X}, dims::Dims) where {X} = rand!(r, Array{X}(undef, dims), X)

[andyferris]

Yeah I think you just have to put in the appropriate methods so you get a default eltype of Float64.

[c42f]

Oh wait, we actually have these rand overrides already so this is certainly a bug in StaticArrays. I'm unsure why we even have them, as we don't seem to do much different from stdlib/Random... maybe to avoid ambiguities.

[chakravala]

I believe this issue can't be solved within StaticArrays but somehow needs to be solved in the rand interface itself.

I've encountered a similar problem in Grassmann package, not sure how to fix it.

[c42f]

I believe we can work around this issue fairly easily by adding making sure we have appropriate overloads in StaticArrays. It won't be an elegant or satisfying solution, but 🤷‍♂️