Increment some components of a vector

[Benoit_Gaudeul]

Hi,
I’m trying to increment some components of a vector (for a finite volume code). Here is a MWE:

components=[1,1,2]
vector=zeros(3)
increment=collect(1:3)
vector[components]+=increment

The value of the I’m hoping for is [3,3,0], the one I have is [2,3,0]. This feature is used almost everywhere in the existing code (where the components appeared only once thus not triggering the bug). I could obviously write a function with a loop to get the desired issue, but I’m afraid that will make the code much slower and slightly less clear. Is there a better way to go?

Best,
Benoît

[stillyslalom]

For-loops are fast in Julia! Don’t shy away from using them.

julia> function increment!(vector, components, Δ)
           for (c, δ) in zip(components, Δ)
               vector[c] += δ
           end
       end
increment! (generic function with 1 method)

julia> using BenchmarkTools

julia> @btime increment!($vector, $components, $increment)
  15.515 ns (0 allocations: 0 bytes)

[mauro3]

I think this is pretty idiomatic:

julia> components = [1,1,2]
3-element Array{Int64,1}:
 1
 1
 2

julia> vector = zeros(3)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> increment= 1:3 # no need to collect
1:3

julia> @view(vector[components]) .+= increment
3-element view(::Array{Float64,1}, [1, 1, 2]) with eltype Float64:
 3.0
 3.0
 3.0

julia> vector
3-element Array{Float64,1}:
 3.0
 3.0
 0.0

[Benoit_Gaudeul]

Thanks a lot for your replies, both ways seem great, and I’m glad to discover this @view functionality.

Could you explain what are views?

Try and modify the result, and you’ll see the difference: julia> A = [1 2; 3 4]; julia> B = A[:,:]; # This is actually a copy of A julia> B[1,1]=0; julia> A 2×2 Array{Int64,2}: 1 2 3 4 julia> B 2×2 Array{Int64,2}: 0 2 3 4 But: julia> B = @view A[:,:]; # B is a view of A julia> B[1,1]=0; julia> A # Ooops! 2×2 Array{Int64,2}: 0 2 3 4 julia> B 2×2 view(::Array{Int64,2}, :, with eltype Int64: 0 2 3 4 In this particular case, it’s the same as B = A, so @view does not se…

[stevengj]

@Benoit_Gaudeul

I could obviously write a function with a loop to get the desired issue, but I’m afraid that will make the code much slower and slightly less clear.

Loops aren’t slow in Julia!

It sounds like you may be translating an existing “vectorized” code from some other language like Matlab or Python? You may not get much benefit from Julia unless you write in a completely different style.

For finite-volume/finite-difference algorithms, my feeling is that a simple loop is typically the clearest as well as the fastest method — just loop over your grid, and apply your stencil each point, e.g. u[i,j] += c[i,j] * (u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4u[i,j]) for a diffusion-like equation. Use ghost cells for boundary conditions.

[Benoit_Gaudeul]

indeed I’m writing code for a nonlinear finite volume scheme and came to Julia to benefit from the automatic differentiation (implicit euler+nonlinear = painful jacobian+newton’s method). I used to use Matlab before, I’ll try to rewrite my code using loops to compare the performances (but not today, today I’ll simply enjoy my now working module for some test cases )

Thanks a lot for the input, it is really nice to be welcomed by such a helpful community