# @simd with array views
**URL:** https://discourse.julialang.org/t/simd-with-array-views/19207
**Category:** General Usage
**Tags:** performance, simd
**Created:** [January 2, 2019, 9:04pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207 "2019-01-02T21:04:11Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Duls91](https://discourse.julialang.org/u/Duls91)
#### Post date: [January 2, 2019, 9:04pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/1 "2019-01-02T21:04:11Z")
Hello everyone and Happy New Year!
I have a decrease in performance in my cross-correlation function when I use array views. I am comparing functions without views, with views and one with the classical imfilter from ImageFiltering.jl:
```julia
using ImageFiltering, Statistics, BenchmarkTools
# basically, slide the template over the padded image and calculate the correlation between template and a slice of the image under it.
function crossCorrelation_raw!(out::Array{T}, image::Array{T,},
template::Array{T, N}, outsize::Tuple{Int64, Int64}) where {T,}
size_template = size(template)
@inbounds for i2=1:outsize[2], i1=1:outsize[1]
# this is the slice of the image under the template for the current iteration
subImage = image[i1:i1+size_template[1]-1, i2:i2+size_template-1]
# caulculate the product for each position and sum them all up, then write to out
s = dot_simd(template, subImage)
@inbounds out[i1, i2] += s
end
return out
end
# the same thing as above, only that subImage is now a view, template is also view passed as an argument in the wrapper function func_view
function crossCorrelation_raw_view!(out::Array{T}, image::Array{T},
template::AbstractArray{T}, outsize::Tuple{Int64, Int64}) where {T}
template_size = size(template)
@inbounds for i2=1:outsize[2], i1=1:outsize[1]
subImage = view(image, i1:i1+template_size[1]-1, i2:i2+template_size[2]-1)
@inbounds out[i1, i2] += dot_simd(template, subImage)
end
return out
end
@inline function dot_simd(x::AbstractArray{T}, y::AbstractArray{T}) where {T}
s = 0.0::T
@simd for i in eachindex(x)
@inbounds s += x[i] * y[i]
end
s
end
# just a wrapper function for the non-views version. Takes the images, prepares image by padding and allocates the output
function func(image::Array{T}, template::Array{T}) where {T}
size_template = size(template)
image_pad = padarray(image, Fill(0.0::T, (size_template[1]-1, size_template[2]-1)))
image_pad = Array{T, 2}(image_pad[Base.axes(image_pad)[1], Base.axes(image_pad)[2]]);
outsize = size(image) .+ size_template.-1
out = zeros(T, outsize)
crossCorrelation_raw!(out, image_pad, template, outsize);
end
# just a wrapper function for the views version. Takes the input images, prepares image by padding and allocates the output
function func_view(image::Array{T, N}, template::Array{T, N}) where {T, N}
size_template = size(template)
template_view = view(template, :, :)
image_pad = padarray(image, Fill(0.0::T, (size_template[1]-1, size_template[2]-1)))
image_pad = Array{T, 2}(image_pad[Base.axes(image_pad)[1], Base.axes(image_pad)[2]]);
outsize = size(image) .+ size_template.-1
out = zeros(T, outsize)
crossCorrelation_raw_view!(out, image_pad, template_view, outsize);
end
# the imfilter version from ImageFiltering.jl
function correlate_image(img::Array{T, N}, template::Array{T, N}) where {T, N}
img = padarray(img, Fill(0, (Int64(round((size(template)[1]-1)/2)), Int64(round((size(template)[2]-1)/2)))))
z = imfilter(T, img, centered(template), Algorithm.FIR())
end
```
Benchmark:
```julia
#creating two images to correlate and a
n = 100
im1 = rand(1.0:2.0, n, n)
im2 = copy(im1);
@btime func(im1, im2)
@btime func_view(im1, im2)
@btime correlate_image(im1, im2)
```
```julia
w/o views: 702.911 ms (79212 allocations: 2.96 GiB)
w/ views: 1.882 s (12 allocations: 2.34 MiB)
imfilter: 616.570 ms (46 allocations: 1.85 MiB)
```
Although there are only 12 allocs for the views version, it takes almost 3 times longer to run. Why is that? w/o views allocates almost 3 GiB of data and is still so much faster than the views function…
Some notes:
- adding @inline to dot\_simd enormously decreased the allocations. Which is good I guess? Only after doing this and reading more about @inline, I got that @simd goes with @inline.
- I benchmarked dot\_simd with both arrays and views and the performance was the same. So dot\_simd should not be the problem here.
- If someone would be kind to check whether my type declarations are “the Julia way” and if they are properly set. Maybe some of them are unnecessary?
Thanks!
Edit: forgot to include the definition of correlate\_image().
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [January 2, 2019, 9:29pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/2 "2019-01-02T21:29:46Z")
Hi,
Note that to perform `SIMD` instructions (single instruction, multiple data), the CPU likes to load contiguous vectors of numbers from the underlying array, and perform the same operation on them.
When you copy an array to make a new one, then all the underlying data is contiguous, so it can be loaded quickly.
However, when you use views, there is no such guarantee.
If you’re comfortable looking at assembly, check this out
```julia
julia> @inline function dot_simd(x::AbstractArray{T}, y::AbstractArray{T}) where {T}
s = 0.0::T
@simd for i in eachindex(x)
@inbounds s += x[i] * y[i]
end
s
end
dot_simd (generic function with 1 method)
julia> X = rand(100,100);
julia> Y = rand(100,100);
julia> Xv = @view X[1:100,1:100];
julia> Yv = @view Y[1:100,1:100];
julia> @code_native dot_simd(X, Y)
```
Here, we see this as the primary loop body:
```julia
L80:
vmovupd (%rcx,%rdi,8), %ymm4
vmovupd 32(%rcx,%rdi,8), %ymm5
vmovupd 64(%rcx,%rdi,8), %ymm6
vmovupd 96(%rcx,%rdi,8), %ymm7
; ││└
; ││┌ @ float.jl:395 within `+'
vfmadd231pd (%rdx,%rdi,8), %ymm4, %ymm0
vfmadd231pd 32(%rdx,%rdi,8), %ymm5, %ymm1
vfmadd231pd 64(%rdx,%rdi,8), %ymm6, %ymm2
vfmadd231pd 96(%rdx,%rdi,8), %ymm7, %ymm3
; ││└
; ││ @ simdloop.jl:74 within `macro expansion'
; ││┌ @ int.jl:53 within `+'
addq $16, %rdi
cmpq %rdi, %rsi
jne L80
```
It says move 4 vectors into `ymm` registers (these registers are 256 bits, meaning move 4 \* 4, or 16 doubles total), and then perform fused multiply add instructions (`vfmadd`) on another region of memory (the other input array) accumulate in other registers. Then perform a comparison, and conditionally jump back to location L80 (start of the loop).
This is what it means for the compiler to “vectorize” the loop.
Unfortunately, this vectorization requires loading contiguous blocks of memory with those `vmovupd` instructions. This is guaranteed with a regular old array, where `A[i+1]` always follows `A[i]`, but not in a view, where who-knows how far apart they can be.
We therefore don’t see any such vectorized portion when we do `@code_native dot_simd(Xv, Yv)`. You can also use BenchmarkTools to confirm that, even though both arrays are actually contiguous in this case, the function on the guaranteed-contiguous arrays is much faster.
EDIT:
If you’d like to recover some of the performance, here is a simple way:
```julia
julia> @inline function dot_simd(x::SubArray{T,2}, y::SubArray{T,2}) where {T}
s = 0.0::T
for i in 1:size(x,2)
s += @views dot_simd(x[:,i], y[:,i])
end
s
end
dot_simd (generic function with 2 methods)
```
These views, because they are not across columns and because the underlying arrays are dense, are guaranteed to be contiguous.
Benchmarks on Ryzen:
```julia
julia> @benchmark dot_simd($X, $Y)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 1.234 μs (0.00% GC)
median time: 1.236 μs (0.00% GC)
mean time: 1.242 μs (0.00% GC)
maximum time: 3.361 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 10
julia> @benchmark dot_simd($Xv, $Yv)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 2.314 μs (0.00% GC)
median time: 2.316 μs (0.00% GC)
mean time: 2.337 μs (0.00% GC)
maximum time: 6.536 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 9
```
and with the new specialized method for matrix-views:
```julia
julia> @benchmark dot_simd($Xv, $Yv)
BenchmarkTools.Trial:
memory estimate: 0 bytes
allocs estimate: 0
--------------
minimum time: 1.630 μs (0.00% GC)
median time: 1.659 μs (0.00% GC)
mean time: 1.682 μs (0.00% GC)
maximum time: 2.511 μs (0.00% GC)
--------------
samples: 10000
evals/sample: 10
```
EDIT:
Yikes! In `@code_native dot_simd(im2_view, im2_view)`
There is not a 4x unrolled block, but there is a vectorized 2x unrolled block, which I’d think really ought to be fine.
```julia
L224:
vmovupd -32(%rdx,%rbp,8), %ymm3
vmovupd (%rdx,%rbp,8), %ymm4
; ││└└└
; ││ @ simdloop.jl:73 within `macro expansion' @ float.jl:395
vfmadd231pd -32(%rdi,%rbp,8), %ymm3, %ymm0
vfmadd231pd (%rdi,%rbp,8), %ymm4, %ymm2
; ││ @ simdloop.jl:74 within `macro expansion'
; ││┌ @ int.jl:53 within `+'
addq $8, %rbp
cmpq %rbp, %rbx
jne L224
```
---
### Author:  [@Duls91](https://discourse.julialang.org/u/Duls91)
#### Post date: [January 2, 2019, 9:58pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/3 "2019-01-02T21:58:21Z")
Thank you @Elrod for your reply and insights. Especially for the detailed explanation of the code\_native output.
I just benchmarked it on my PC:
```julia
n = 10000
im1 = rand(1.0:2.0, n, n)
im1_view = view(im1, :,:);
@btime dot_simd(im1, im1)
@btime dot_simd(im1_view, im1_view)
@btime dot_simd_new(im1_view, im1_view) # this is wwith your added dot_simd version
43.473 ms (1 allocation: 16 bytes)
42.917 ms (1 allocation: 16 bytes)
42.961 ms (1 allocation: 16 bytes)
```
In all cases, it is comparable and there is no trouble with vectorization. Is this is what you meant by
> even though both arrays are actually contiguous in this case
So for some reason, dot\_simd is also vectorizing the views (I also checked and compared the code\_native output). But why do I see such a discrepancy when I use dot\_simd in the cross-correlation.
Maybe in this specific case, dot\_simd is not the problem but something else? Sorry if I am being a hassle…
---
### Author:  [@Duls91](https://discourse.julialang.org/u/Duls91)
#### Post date: [January 2, 2019, 10:17pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/4 "2019-01-02T22:17:07Z")
OK, I added your dot\_simd to crossCorrelation\_raw\_view():
Benchmark:
```julia
n = 100
im1 = rand(1.0:2.0, n, n)
im2 = copy(im1)
@btime func(im1, im2)
@btime func_view(im1, im2)
@btime correlate_image(im1, im2)
624.448 ms (79212 allocations: 2.96 GiB)
130.391 ms (39612 allocations: 4.75 MiB)
540.183 ms (46 allocations: 1.85 MiB)
```
It did solve the issue. Simd and views are compatible. But I am just wondering the same as above: Testing dot\_simd with arrays or views gives the same performance. So why does it not keep that performance in the parent crossCorrelation\_raw/\_view functions?
Either way, This seems to work now and I am glad that it does. Hope I can translate it to some other simd functions I have. Thanks a lot again!
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [January 2, 2019, 10:57pm UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/5 "2019-01-02T22:57:45Z")
There are two factors at play here.
1. Memory bandwidth
2. SubArrays are smart.
`:` has a different type from `a:b`. Therefore, when you create a view with `:` for the rows, the `SubArray` dispatch knows the memory is actually contiguous:
```julia
julia> n = 224;
julia> im1 = rand(1.0:2.0, n, n);
julia> im1_view = @view im1[:,:];
julia> im2_view = @view im1[1:end,:];
julia> IndexStyle(typeof(im1_view))
IndexLinear()
julia> IndexStyle(typeof(im2_view))
IndexCartesian()
```
In fact, when we use `@code_native`, we can see that calling on `im1_view` is vectorized, but that `im2_view` is not. `im2_view` is like the view in your function.
However, that change alone does not solve the issue. That is because, when you choose `n=10000`, the speed is constrained by memory bandwidth, not by CPU speed. Both are using the same amount of data, so it takes the same time to move memory around.
Shrinking `n` lets CPU time become the more important factor. I chose `n=224` because 224 is the second largest number divisible by 32 and smaller than `sqrt(512e3/8)`. The L2 cache on the computer I’m benchmarking on is 512e3, so this should fit into the cache of a single core. Second largest for a little head room (I think the 512 KiB is for data and instruction cache? Or does the instruction cache go straight to L1?).
This yields:
```julia
julia> @btime dot_simd($im1, $im1)
6.124 μs (0 allocations: 0 bytes)
125230.0
julia> @btime dot_simd($im1_view, $im1_view)
6.115 μs (0 allocations: 0 bytes)
125230.0
julia> @btime dot_simd_new($im1_view, $im1_view)
6.811 μs (0 allocations: 0 bytes)
125230.0
julia> @btime dot_simd($im2_view, $im2_view)
8.493 μs (0 allocations: 0 bytes)
125230.0
julia> @btime dot_simd_new($im2_view, $im2_view)
6.795 μs (0 allocations: 0 bytes)
125230.0
```
Note the difference in types of the two views:
```julia
julia> typeof(im1_view)
SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
julia> typeof(im2_view)
SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}
```
You could dispatch to the original vs the new `dot_simd` based on whether the last parameter is true or false.
---
### Author:  [@Duls91](https://discourse.julialang.org/u/Duls91)
#### Post date: [January 4, 2019, 1:17am UTC](https://discourse.julialang.org/t/simd-with-array-views/19207/6 "2019-01-04T01:17:19Z")
Thanks again for all the details. Turns out even more complicated then I anticipated. Even the small, but important, detail of the different types from : and a:b.