# Most efficient way of adding elements within matrices in loops

**URL:** <https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789>\
**Category:** Performance\
**Tags:** question, matrices, loopvectorization\
**Created:** [November 18, 2024, 9:14pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789 "2024-11-18T21:14:55Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![BMI\_OR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bmi_or/32/213496_2.png) [@BMI\_OR](https://discourse.julialang.org/u/BMI_OR)\
**Post date:** [November 18, 2024, 9:14pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789/1 "2024-11-18T21:14:55Z")

</div>

I have a vector `A` and I need to update its values with the ones of matrix `B`. Each element of `B` needs to be added to some values in `A`. The mapping of indeces from matrix to vector is given by a matrix of vectors `C`. I have a piece of code that currently does the work but it is inefficient:

```julia
for i in 1:I, j in 1:J, k in 1:K
    vector_indices_to_update = C[i,j,k];
    A[vector_indices_to_update] .+= B[i,j,k]
end

```

What is the most efficient way to implement this operation?  
As additional info, the length of vector A is significantly larger that the size of the dimensions of the matrices,  
Thanks!

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [November 18, 2024, 9:41pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789/2 "2024-11-18T21:41:09Z")

</div>

Welcome to this forum!

We can better engage with your question if you [follow the PSA](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757/1) on questions.

The thing that your example is missing is “self contained runnable + mimics your true usecase”. To give an example in style how this could look like:

```julia
using BenchmarkTools

I=10;J=11;K=12; N=500; A = rand(N); C=rand(1:N, I, J, K); B=rand(I,J,K);

function update_indices!(A, B, C)
  for i in 1:I, j in 1:J, k in 1:K
    vector_indices_to_update = C[i,j,k];
    A[vector_indices_to_update] += B[i,j,k]
  end
end

@btime update_indices!($A, $B, $C)

```

What you need to do in order to get useful help is write code snippets that initialize A,B,C. “the length of A is significantly larger…” is all fine and dandy, but it is better to simply incorporate that into the sample problem you show us.

Please also think about all the details, and make sure that they all mirror your true usecase. If your real data is not sorted, then you SHOULD NOT give sorted samples; if your real data is Float32 then you give us Float32, etc etc.

Don’t be go :surprised\_pikachu: if you post `collect(1:10_000)` as a test vector and proposed solutions delete the `collect`, a la

```julia
julia> @btime sum($(1:1000_000))
  2.016 ns (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [November 18, 2024, 10:38pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789/3 "2024-11-18T22:38:59Z")

</div>

> [@BMI\_OR](#):
>
> I have a piece of code that currently does the work but it is inefficient:
> 
> ```julia
> for i in 1:I, j in 1:J, k in 1:K
> vector_indices_to_update = C[i,j,k];
> A[vector_indices_to_update] .+= B[i,j,k]
> end
> 
> ```
> 
> What is the most efficient way to implement this operation?

This looks very efficient to me, provided `I`, `J`, `K`, `A`, `B`, and `C` are either local variables, or `const` global variables, and `A` and `B` have the same element type. Otherwise, there may be considerable performance penalties. Preferably, the vectors in `C` should be sorted, but I’m uncertain how much that matters. If you are sure no index is out of bounds you can add `@inbounds` before the `for`.

Just one thing. The `A[v] .+= B[i,j,k]` translates to `A[v] .= A[v] .+ B[i,j,k]`, so when `v` is a vector, there will be a copy on the right hand side. You can do `@view(A[v]) .+= B[i,j,k]` to avoid that.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 18, 2024, 11:27pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789/4 "2024-11-18T23:27:17Z")

</div>

> [@BMI\_OR](#):
>
> ```julia
> for i in 1:I, j in 1:J, k in 1:K
> vector_indices_to_update = C[i,j,k];
> A[vector_indices_to_update] .+= B[i,j,k]
> end
> 
> ```

If each `vector_indices_to_update` is a vector, then defitely add `@views` or iterate over it in a loop.

Also, your iteration order is suboptimal, since Julia arrays are column major. Either change the loop order, or, even better, get the correct order automatically:

```julia
for i in eachindex(B, C) 
    vector_indices_to_update = C[i]
    @views A[vector_indices_to_update] .+= B[i]
end

```

or

```julia
for i in (val, inds) in zip(B, C)
    @views A[inds] .+= val
end

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 18, 2024, 11:31pm UTC](https://discourse.julialang.org/t/most-efficient-way-of-adding-elements-within-matrices-in-loops/122789/5 "2024-11-18T23:31:41Z")

</div>

> [@sgaure](#):
>
> If you are sure no index is out of bounds you can add `@inbounds` before the `for`.

I don’t think that is going to be safe in this case, since loaded indices don’t have the same guarantees as the ones coming from `eachindex`.