# Usage of vector of arrays

Hi everyone,
I am trying to manipulate data using a vector of arrays. The aim is to store the connectivity of a Poisson finite difference stencil using vector of arrays, rather than storing them into independent vectors.
To this end, I’ve initialised a data structure `Tsten` that contains at entry `[1]` indices of center nodes, entries `[2]` and `[3]` are west and east nodes, respectively.
In the “independent vectors” version I use correspondingly `TC`, `TW` and `TE`.

1. I want to populate the first entry of `Tsten` (Tsten[1] .= T) and this, introduce entries in `Tsten[2]` and `Tsten[3]`, which is not what I would have expected. Any idea of what is wrong there?

2. I want to populate `Tsten[2]` and `Tsten[3]` using the function `PoissonSparsity1D_v1!`
and I expect same results than using an “independent vectors” approach using `PoissonSparsity1D!`. But it is not the case. It is likely related to the previous point. Would anyone have some hints?

3. Finally, it would be nice to ‘merge’ `PoissonSparsity1D_v1` and `PoissonSparsity1D` using multiple dispatch (such that a single function apply to either a vector of arrays or independent vectors). But this is a next step…

Cheers!

``````function PoissonSparsity1D!(TC, TW, TE)
TW[2:end-0] .= TC[1:end-1];
TE[1:end-1] .= TC[2:end-0];
end

function PoissonSparsity1D_v1!(x)
TC, TW, TE = x
TW[2:end-0] .= TC[1:end-1];
TE[1:end-1] .= TC[2:end-0];
end

function PrintAll(Tsten, TC, TW, TE)
println( "Center:  ", TC )
println( "Center:  ", Tsten[1] )
println( "West:    ", TW )
println( "West:    ", Tsten[2] )
println( "East:    ", TE )
println( "East:    ", Tsten[3] )
end

function main()

# Initialise
ncx = 5
TC, TW, TE = fill(0,  ncx), fill(0,   ncx), fill(0,   ncx) # independent arrays
T = collect(1:ncx)

# Define and allocate vector of arrays
Tsten = Array{Vector{Int64}}(undef, 3);
fill!(Tsten, zeros(ncx))

# Initially all entries are zero
println("Check #1")
PrintAll(Tsten, TC, TW, TE)

# Now let's just put non zero values in the first vector
Tsten[1] .= T

# Why are there non-zeros in Tsten[i] with i>1?
println("Check #2 -  Why are there non-zeros in Tsten[i] with i>1?")
PrintAll(Tsten, TC, TW, TE)

# Here, I try to generate the connectivity - similar results are expected
TC .= T
PoissonSparsity1D!(TC,TW,TE)
Tsten[1] .= T
PoissonSparsity1D_v1!(Tsten)

println("Check #3 - Then, results also differ there:")
PrintAll(Tsten, TC, TW, TE)

end

for it=1:1
@time main()
end
``````

Regarding question 1, be careful about initializing a vector of vectors using `fill`.
This is a confusing aspect of Julia at first → link to post where this is explained.

All 3 entries of `Tsten` will point to the same location and the inplace assignment `Tsten[1] .= T` will keep all 3 pointing to the new replaced `T` values.

Do instead: `Tsten[1] = T`, or keep the inplace assignment but then you will have to pre-assign the individual arrays of zeros independently using a loop or comprehension as in the linked post .

1 Like

Wow, thanks!
That was not super intuitive to me. Sorry for not having found the related post beforehand.
Replacing the lines:

``````Tsten = Array{Vector{Int64}}(undef, 3);
fill!(Tsten, zeros(ncx))
``````

by a comprehension as:

``````Tsten= [zeros(Int64,ncx) for _ in 1:3]
``````

solved point #1 as well a #2. I can tackle point #3 now.

Be careful because I assigning the complete vectors will only relabel them, if you are working with mutable arrays. For example, here all vectors are the same:

``````julia> a = [1,2]
2-element Vector{Int64}:
1
2

julia> b = [ a, a ]
2-element Vector{Vector{Int64}}:
[1, 2]
[1, 2]

julia> b[1][1] = 0
0

julia> a
2-element Vector{Int64}:
0
2

julia> b
2-element Vector{Vector{Int64}}:
[0, 2]
[0, 2]

``````

To avoid this kind of thing you have to explicitly copy the arrays, as in:

``````julia> b = [ copy(a), copy(a) ]
2-element Vector{Vector{Int64}}:
[0, 2]
[0, 2]

julia> b[1][1] = 1
1

julia> b
2-element Vector{Vector{Int64}}:
[1, 2]
[0, 2]

julia> a
2-element Vector{Int64}:
0
2

``````

But it may be that in this case the best solution is to use immutable arrays (from StaticArrays):

``````julia> using StaticArrays

julia> a = @SVector [1,2]
2-element SVector{2, Int64} with indices SOneTo(2):
1
2

julia> b = [ a, a ]
2-element Vector{SVector{2, Int64}}:
[1, 2]
[1, 2]

julia> b[1] = @SVector [0,0]
2-element SVector{2, Int64} with indices SOneTo(2):
0
0

julia> a
2-element SVector{2, Int64} with indices SOneTo(2):
1
2

julia> b
2-element Vector{SVector{2, Int64}}:
[0, 0]
[1, 2]

``````
1 Like