# Loop vs vectorization
**URL:** https://discourse.julialang.org/t/loop-vs-vectorization/104796
**Category:** Performance
**Created:** [October 10, 2023, 11:21am UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796 "2023-10-10T11:21:02Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Alberto\_Roman](https://discourse.julialang.org/u/Alberto_Roman)
#### Post date: [October 10, 2023, 11:21am UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796/1 "2023-10-10T11:21:02Z")
Hi all,
I have a partial differential equation to solve using finite difference. I am trying to understand if it is better to use a vectorized version or a loop version to calculate the right hand side of the equation, which has the spatial discretization. What I am struggling to understand is that for a given grid size (e.g. number of points at which the RHS is evaluated) the vector version performs more poorly than the loop version.
However if I increase the grid size the difference in performance between the two seems to become less.
In order to improve my understanding if Julia performs better with loop or vectorized code I wrote the following code
```julia
function vect_sum(A,B)
C .= A + B;
return C
end
function loop_sum(A,B,C,m,n)
for i=1:m
for j=1:n
C[i,j] = A[i,j] + B[i,j]
end
end
return C
end
n= 1000
m = 1000
A = rand(m,n)
B = rand(m,n)
C= zeros(m,n);
@btime vect_sum(A,B);
C= zeros(m,n);
@btime loop_sum(A,B,C,m,n);
```
Here, despite the common sense that loop in Julia are faster, I am getting that the vector version of the code perform better than the loop.
Does anybody have an hint on why this is happening?
Thanks
---
### Author:  [@ufechner7](https://discourse.julialang.org/u/ufechner7)
#### Post date: [October 10, 2023, 11:29am UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796/2 "2023-10-10T11:29:52Z")
Can you try to fix your backticks that surround your code?
They should look like this
```
On an English keyboard that is often the top left key. On some keyboards you have to press it twice.
Or copy-and-paste mine…
---
### Author:  [@maxfreu](https://discourse.julialang.org/u/maxfreu)
#### Post date: [October 10, 2023, 12:05pm UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796/3 "2023-10-10T12:05:56Z")
In theory, both should be equally fast. In julia it shouldn’t matter whether you write loops or numpy-style “vectorized” functions, because the loops are compiled. After that, they should do pretty much the same thing as the “vectorized” code, which at some point is also just a bunch of loops. Nevertheless, you can learn some things even on the your simple use case - here are some remarks:
1. The benchmark is not fair, because `vect_sum` allocates the output vector C, which takes time. `loop_sum` only writes to the already allocated output C that you pass in. A good rule is to avoid allocations whenever you can, because this can make your code much faster, depending on the exact problem.
2. The reason why `loop_sum` probably is slower is because it performs bounds checks. So it checks whether the index you use to access the array is within allowed limits. This involves comparison operations that take time. You can tell the compiler that all array accesses are within bound by using the `@inbounds` macro either in front of an expression (`@inbounds C[i,j] = ...`) or in front of a loop / scope (`@inbounds for ... or @inbounds begin ...`). That speeds things up a lot, but if you get things wrong, you get a segmentation fault and your program crashes.
3. Here, A and B have the same size. Therefore you can simply index into them with a linear index like so:
```julia
for i in eachindex(A)
C[i] = A[i] + B[i]
end
```
1. This has the benefit of traversing your computer’s memory in a linear fashion, which speeds up things even more.
2. Speaking about memory access order: julia uses a column major memory layout. So when you perform linear indexing as just described, you traverse the first column of your matrix from top to bottom, row by row, followed by the second column and so on. The rows in the matrix are indexed by the first index in the brackets (here i) and the columns by the second. So to achieve linear memory access with nested loops, you either have to switch the order of the loops or the i and j in the brackets. A good rule is that the index of the innermost loop comes first in the array access, followed by the outer ones, like so:
```julia
for i ...
for j ...
for k ...
C[k,j,i] = ...
```
1. You don’t have to explicitly tell `loop_sum` how big your arrays are by passing m and n. This is not C (: So you can either do something like `for i in 1:size(A,1) ...` to traverse the first dimension of A. But it is much better to use other functions, such as `eachindex`, `axes` or `CartesianIndices`, as they can make your life easier and safer once you understood how to use them. So I highly recommend reading their docs.
2. You already used the `BenchmarkTools` package, which is good. It is good practice to “interpolate” the variables using $ to “pre-compute” them before benchmarking, which keeps your results clean from allocations that otherwise could occur. So write `@btime vect_sum($A, $B)` - however, here it shouldn’t make much of a difference.
Try to implement these hints and then see how it turns out - I’m curious 🙂
Edit: `vect_sum` already implements the right array access order and inbounds accesses, ~~which is probably why it is faster~~ but seemingly this does not outweigh the allocation, which is why it’s slower.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [October 10, 2023, 12:22pm UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796/4 "2023-10-10T12:22:28Z")
As said above, the reason is absolutely the order your are running the loops:
```julia
julia> @btime vect_sum(A,B);
1.799 ms (3 allocations: 7.63 MiB)
julia> C= zeros(m,n);
julia> @btime loop_sum(A,B,C,m,n);
13.313 ms (0 allocations: 0 bytes)
julia> function loop_sum2(A,B,C,m,n)
for j=1:n
for i=1:m
C[i,j] = A[i,j] + B[i,j]
end
end
return C
end
loop_sum2 (generic function with 1 method)
julia> @btime loop_sum2(A,B,C,m,n);
1.160 ms (0 allocations: 0 bytes)
```
But also note more Julian implementations could be:
```julia
julia> function loop3!(C,A,B)
for i in eachindex(C,A,B)
C[i] = A[i] + B[i]
end
return C
end
loop3! (generic function with 1 method)
julia> @btime loop3!($C,$A,$B);
1.172 ms (0 allocations: 0 bytes)
julia> function loop4!(C,A,B)
for j in axes(C,2), i in axes(C,1)
C[i,j] = A[i,j] + B[i,j]
end
return C
end
loop4! (generic function with 1 method)
julia> @btime loop4!($C,$A,$B);
1.185 ms (0 allocations: 0 bytes)
```
But also note that your vectorized version is not exactly comparable, because there you were allocating some intermediates (and the result). You get the same performance with:
```julia
julia> function vect_sum!(C,A,B)
C .= A .+ B
return C
end
vect_sum! (generic function with 1 method)
julia> @btime vect_sum!($C,$A,$B);
1.192 ms (0 allocations: 0 bytes)
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [October 10, 2023, 12:40pm UTC](https://discourse.julialang.org/t/loop-vs-vectorization/104796/5 "2023-10-10T12:40:15Z")
> [@Alberto\_Roman](#):
>
> I have a partial differential equation to solve using finite difference. I am trying to understand if it is better to use a vectorized version or a loop version to calculate the right hand side of the equation, which has the spatial discretization.
For finite differences, I would tend to prefer loops (properly written, as noted by @lmiq and others above). Loops just give you more flexibility for this kind of thing (as well as excellent performance).
For example, see this post for example code of writing a scalar-wave equation finite-difference routine using loops that, combined with Julia’s [`CartesianIndex` facilities](https://julialang.org/blog/2016/02/iteration/), support _arbitrary dimensionality_ with a _single_ loop: [Seemingly unnecessary allocation within for loop - #8 by stevengj](https://discourse.julialang.org/t/seemingly-unnecessary-allocation-within-for-loop/72556/8)