Usage of vector of arrays

General Usage
1

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
2

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
3

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.

4

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
5

Thanks a lot for advices.
Ok, so using StaticArray one gets MATLAB behaviour while, without, it is rather similar to Python/numpy.