Generate function to build a static array from vector of vector

I have a Vector of Vectors (actually Vector of SVectors, but I do not believe this is relevant to the question.) All vectors, inner and outer, have length N. I wish to generate an SMatrix from it

  • in a typestable manner
  • without allocation
    ('tis better t’was done quickly)

So I create a generated function, here’s the MWE

using StaticArrays, Printf
@generated function foo(a,::Val{N}) where{N} 
    els =  :()
    for j = 1:N,i=1:N
        els = :($els...,a[$i][$j])
    end
    return quote
        SMatrix{N,N}($els) 
    end
end
@code_warntype foo([randn(2),randn(2)],Val(2))

which returns values as expected and is typestable, but performance is atrocious.

On the other hand if I hardcode

foo(a,::Val{2}) = @SMatrix{2,2}(a[1][1], ... ,a[2][2])

(in the real world, N=8, fun!) I get great performance. Hence
a) the code I want to generate is the right code
b) but that is not what I do generate.
c) afterthought: or generation is triggered more often that I though - this would then not be reflected in the MWE.

I suspect I am generating at each call, a code containing the actual numerical values (The Julia manual does not guarantee that a generated function is compiled only when called with new signature)

@macroexpand does not work here - I there a way I can visualize my generated code? And where did I go wrong in my metacode?

I would just use SizedArray to hint this:

foo(a, ::Val{N}) where N = reduce(hcat, SizedVector{N}(a))
1 Like

Not possible, since the length of a Vector is not part of its type, while the size of an SMatrix is part of its type. You need a different type of container to achieve this (such as the SizedArray mentioned by @Tamas_Papp above), but of course that is only possible of the length of your Vector is known statically.

I forgot to mention - the size is known staticaly - actualy my subvector are Static, and I know the outer Vector to have the same length.

1 Like

foo(a, ::Val{N}) where N = reduce(hcat, SizedVector{N}(a))

Works like a charm, very fast, and elegant. Brilliant, thanks!

2 Likes

Just in case someone reads this in the future: As coded above, the border case N==1 fails: hcat is never called, and the function returns a SVector{1} instead of the wanted SMatrix{1,1}.

This is easily fixed:
foo(a, ::Val{N}) where N = SMatrix{N,N}(reduce(hcat, SizedVector{N}(a)))