SVector from concat of subvectors

General Usage
1

Greetings.

I am seeking to replace

# generates local stiffness matrix 
function genLocStiffMat(element::Element)
    h     = element.area
    e1    = element.e1
    e2    = element.e2
    Iloc  = SVector(e1, e1, e2, e2) 
    Jloc  = SVector(e1, e2, e1, e2) 
    Aloc  = SVector(1/h, -1/h, -1/h, 1/h) 
    return Iloc, Jloc, Aloc
end

by a more elegant version that works for larger sets of e_i values.

The version

function genLocStiffMat2(element::Element)
    h     = element.area
    v     = SVector(element.e1, element.e2)
    Iloc  = SVector{4}([v[i] for j=1:2, i=1:2]) 
    Jloc  = SVector{4}([v[i] for i=1:2, j=1:2])
    Aloc  = SVector(1/h, -1/h, -1/h, 1/h) 
    return Iloc, Jloc, Aloc
end

allocates memory. Advice is appreciated.

2

try using generator instead […]

function genLocStiffMat2(element)
    h     = element.area
    v     = SVector(element.e1, element.e2)
    Iloc  = SVector{4}(v[i] for j=1:2, i=1:2) 
    Jloc  = SVector{4}(v[i] for i=1:2, j=1:2)
    Aloc  = SVector(1/h, -1/h, -1/h, 1/h) 
    return Iloc, Jloc, Aloc
end
1 Like
3

I don’t have access to running Julia at the moment, but the key is that the statements [v[i] for j=1:2, i=1:2] and [v[i] for i=1:2, j=1:2] both allocate a 2x2 matrix. I think what you want is

Iloc  = SVector{4}(v[i] for j=1:2 for i=1:2)

though I may have gotten my is and js backwards.

1 Like
4

Hurray! What is “generator” is this context?

5

Yup! That does the trick. I would like to understand the underlying principle. Why does [ … ] allocate memory? How is it avoided by the suggestion you provide? Sincere thanks.

6

Comprehensions can also be written without the enclosing square brackets, producing an object known as a Generator

1 Like
7

Hurray again! Thx for bearing with me.

8

[...] is julias syntax for creating an Array, which is always (or close enough to always, there may be some exceptions, I’m not an expert) heap allocated (reported as allocations by @time). Thus, your original expression first created an Array and then SVector used this allocated Array to construct the SVector.

1 Like
9

Playing in the REPL should be helpful:

julia> [i for j=1:2, i=1:2]
2Ă—2 Matrix{Int64}:
 1  2
 1  2

julia> isbits(ans)
false

julia> (i for j=1:2, i=1:2)
Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, var"#3#4"}(var"#3#4"(), Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((1:2, 1:2)))

julia> isbits(ans)
true

I’d also note that Iterators.map is an alternative to using the generator syntax sugar:

julia> a = (i for j=1:2, i=1:2)
Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, var"#1#2"}(var"#1#2"(), Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((1:2, 1:2)))

julia> b = Iterators.map(last, Iterators.product(1:2, 1:2))
Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, typeof(last)}(last, Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((1:2, 1:2)))

julia> a == b
false

julia> collect(a) == collect(b)
true

(TBH I’m not sure why a == b returns false? EDIT: Base.Generator doesn’t overload ==, so they get compared by ===. Bug? Reported an issue on the Github: equality for iterators · Issue #53046 · JuliaLang/julia · GitHub)

10

a shape ALMOST the same as a=...

julia> i=1:2
1:2

julia> j=1:2
1:2

julia> c=Base.Generator((j,i)->i,Base.product(j,i))
Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, var"#41#42"}(var"#41#42"(), 
Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((1:2, 1:2)))

julia> a = (i for j=1:2, i=1:2)
Base.Generator{Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}, var"#43#44"}(var"#43#44"(), 
Base.Iterators.ProductIterator{Tuple{UnitRange{Int64}, UnitRange{Int64}}}((1:2, 1:2)))
11

There is also the @SVector convenience macro:

Iloc  = @SVector [e1, e1, e2, e2]
Jloc  = @SVector [e1, e2, e1, e2]
Aloc  = @SVector [1/h, -1/h, -1/h, 1/h]

which doesn’t allocate and looks quite clear.