# Speed up multidimensional code without Base.Cartesian
**URL:** https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328
**Category:** Performance
**Tags:** question
**Created:** [September 17, 2021, 11:13am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328 "2021-09-17T11:13:27Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 17, 2021, 11:13am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/1 "2021-09-17T11:13:28Z")
Hello, I wonder if someone could help me speeding up my arbitrary dimension code to run as fast as the explicit version but without using `Base.Cartesian`
Consider the following examples
```julia
using Random
using BenchmarkTools
@inline bc(x, L) = x < L ? x + 1 : 1
function energy_explicit(Λ::AbstractArray{T,3}) where T
L = size(Λ, 1)
E = 0
@simd for r ∈ CartesianIndices(Λ)
x, y, z = Tuple(r)
@inbounds E += Λ[x, y, z] * (Λ[bc(x, L), y, z] + Λ[x, bc(y, L), z] + Λ[x, y, bc(z, L)])
end
return -E
end
function energy_arbitrary(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@simd for r ∈ CartesianIndices(Λ)
@inbounds E += Λ[r] * sum(Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...] for d ∈ 1:N)
end
return -E
end
```
On my machine I get the following results
```julia
julia> Λ = rand(Int8[-1,1], 64, 64, 64);
julia> @btime energy_explicit($Λ)
252.054 μs (0 allocations: 0 bytes)
-192
julia> @btime energy_arbitrary($Λ)
498.445 μs (0 allocations: 0 bytes)
-192
```
I know that I could use `Base.Cartesian` and explicitly generate the expression from `energy_explicit` but can someone come up with a more simple solution?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [September 17, 2021, 12:09pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/2 "2021-09-17T12:09:04Z")
I think what you want to do is make `Λ` be an `Array{StaticVector{T,L},N-1}` which will make the `ntuple` construction type stable.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 17, 2021, 12:11pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/3 "2021-09-17T12:11:04Z")
Could you explain how to use it exactly and the second question is: why is ntuple not type stable in the first place? I mean, the `N` is known from dispatch…
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 17, 2021, 12:11pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/4 "2021-09-17T12:11:43Z")
For comparison reasons, here the generated version:
```julia
@generated function energy_generated(Λ::AbstractArray{T,N}) where {T,N}
quote
E = 0
@simd for r ∈ CartesianIndices(Λ)
E += @inbounds Λ[r] * Base.Cartesian.@ncall($N, +, d -> Base.Cartesian.@nref($N, Λ, x -> (x == d ? bc(r[x], size(Λ, x)) : r[x])))
end
return -E
end
end
```
```julia
julia> @btime energy_generated($Λ)
252.085 μs (0 allocations: 0 bytes)
564
```
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [September 17, 2021, 12:23pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/5 "2021-09-17T12:23:38Z")
> [@stakaz](#):
>
> `for d ∈ 1:N)`
I think you want to replace this with an `ntuple`.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 17, 2021, 12:27pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/6 "2021-09-17T12:27:16Z")
> [@mcabbott](#):
>
> I think you want to replace this with an `ntuple` .
```julia
function energy_arbitrary_ntuple(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@simd for r ∈ CartesianIndices(Λ)
@inbounds E += @inbounds Λ[r] * (+)(ntuple(d -> Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...], N)...)
end
return -E
end
```
```julia
julia> @btime energy_arbitrary_ntuple($Λ)
498.140 μs (0 allocations: 0 bytes)
564
```
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [September 17, 2021, 12:55pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/7 "2021-09-17T12:55:47Z")
What if you write a custom (simple) sum function? Base sum is not the fastest, and it does not specialize for N either (AFAIK).
---
### Author:  [@jipolanco](https://discourse.julialang.org/u/jipolanco)
#### Post date: [September 17, 2021, 1:43pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/8 "2021-09-17T13:43:55Z")
Following @lmiq’s suggestion, I was able to recover the same performance of the explicit version by defining a sum function that specialises on the dimension `N`:
```julia
@inline nsum(f::F, ::Val{N}) where {F,N} = sum(ntuple(f, Val(N)))
function energy_nsum(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@simd for r ∈ CartesianIndices(Λ)
@inbounds E += Λ[r] *
nsum(Val(N)) do d
@inbounds Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...]
end
end
return -E
end
```
```julia
ulia> @btime energy_explicit($Λ)
263.647 μs (0 allocations: 0 bytes)
-832
julia> @btime energy_nsum($Λ)
263.547 μs (0 allocations: 0 bytes)
-832
```
---
### Author:  [@jipolanco](https://discourse.julialang.org/u/jipolanco)
#### Post date: [September 17, 2021, 1:58pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/9 "2021-09-17T13:58:02Z")
> [@stakaz](#):
>
> ```julia
> function energy_arbitrary_ntuple(Λ::AbstractArray{T,N}) where {T,N}
> E = 0
> @simd for r ∈ CartesianIndices(Λ)
> @inbounds E += @inbounds Λ[r] * (+)(ntuple(d -> Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...], N)...)
> end
> return -E
> end
>
> ```
I think the issue with this solution is that the `@inbounds` is not propagated to the interior of the `ntuple`s. Otherwise, it should be exactly equivalent to my solution above.
If I throw in some additional `@inbounds` (and rewrite it to make it a little easier to read), I can also achieve the performance of the explicit version:
```julia
function energy_arbitrary_ntuple(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@simd for r ∈ CartesianIndices(Λ)
E += @inbounds Λ[r] * (+)(ntuple(N) do d
inds = ntuple(N) do i
@inbounds i == d ? bc(r[i], size(Λ, i)) : r[i]
end
@inbounds Λ[inds...]
end...)
end
return -E
end
```
```julia
julia> @btime energy_arbitrary_ntuple($Λ)
262.834 μs (0 allocations: 0 bytes)
-832
```
---
### Author:  [@N5N3](https://discourse.julialang.org/u/N5N3)
#### Post date: [September 17, 2021, 2:42pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/10 "2021-09-17T14:42:13Z")
The performance difference of the original version comes from:
1. boundcheck (@inbounds won’t work if used outside the generator)
2. the length of generator is unknow.
The following modified version has a similar performance.
```julia
function energy_arbitrary_faster(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@inbounds @simd for r ∈ CartesianIndices(Λ)
E += Λ[r] * sum((@inbounds Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...]) for d in ntuple(identity, N))
end
return -E
end
```
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [September 17, 2021, 3:21pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/11 "2021-09-17T15:21:07Z")
The fact that `bc` makes the indexing unpredictable is also about a factor of 2 effect, for me. Can you duplicate the first row/col/etc of `Λ` as the last+1? Otherwise, you can deal with the edge afterwards. Not hard to write for `N=2`, more complicated for N=3… in general you could try [TiledIteration.EdgeIterator](https://github.com/JuliaArrays/TiledIteration.jl#edgeiterator) although I haven’t timed it here.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 17, 2021, 4:48pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/12 "2021-09-17T16:48:45Z")
> [@lmiq](#):
>
> What if you write a custom (simple) sum function? Base sum is not the fastest, and it does not specialize for N either (AFAIK).
It specializes/inlines if you pass a `tuple` (whose length is known at compile-time) instead of a generator.
---
### Author:  [@Mason](https://discourse.julialang.org/u/Mason)
#### Post date: [September 17, 2021, 5:13pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/13 "2021-09-17T17:13:22Z")
> [@Oscar\_Smith](#):
>
> I think what you want to do is make `Λ` be an `Array{StaticVector{T,L},N-1}` which will make the `ntuple` construction type stable.
Another nice option is [HybridArrays.jl](https://github.com/mateuszbaran/HybridArrays.jl) which are basically what you suggest, but pretend to be `N` dimensional, with some of the dimensions statically sized.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [September 17, 2021, 5:43pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/14 "2021-09-17T17:43:42Z")
> [@stakaz](#):
>
> `@inline bc(x, L) = x < L ? x + 1 : 1 `
Note that this is a perfect case for using [ghost cells](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/2). Just make your array `Λ` bigger by one element in each direction, and copy the periodic boundaries before doing any “stencil” calculations.
Not only is this much more flexible and tends to lead to cleaner code (it separates your boundary conditions from your other looping code), but it is also significantly faster because you eliminate all of the `?` conditionals in your inner loops.
In particular, the following code gives the same answer as `energy_explicit` but uses ghost cells:
```julia
Λghost = similar(Λ, size(Λ) .+ 1);
Λghost[CartesianIndices(Λ)] = Λ;
function periodicboundaries!(Λghost::AbstractArray{<:Any,N}) where {N}
sz = size(Λghost)
for d = 1:N
from = CartesianIndices(ntuple(i -> i == d ? (1:1) : (1:sz[i]), N))
to = CartesianIndices(ntuple(i -> i == d ? (sz[i]:sz[i]) : (1:sz[i]), N))
for (i,j) in zip(from,to)
Λghost[j] = Λghost[i]
end
end
return Λghost
end
function energy_arbitrary_faster_ghost(Λ::AbstractArray{T,N}) where {T,N}
E = 0 + zero(T) # promote to at least Int
@inbounds @simd for r ∈ CartesianIndices(ntuple(i -> 1:size(Λ,i)-1, N))
E += Λ[r] * sum(ntuple(d -> @inbounds(Λ[ntuple(i -> i == d ? r[i]+1 : r[i], N)...]), N))
end
return -E
end
@btime periodicboundaries!($Λghost)
@btime energy_arbitrary_faster_ghost($Λghost)
```
giving
```julia
57.329 μs (0 allocations: 0 bytes)
91.855 μs (0 allocations: 0 bytes)
```
vs 337µm for `energy_explicit` on my machine. So, even if you count the cost of copying over the periodic boundaries (which only needs to be done once whenever you update the array) it is more than 2× faster. If you don’t count the time for `periodicboundaries!`, then it is almost 4× faster.
PS. Note that I eliminated the generator completely, compared to @N5N3’s solution.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 17, 2021, 6:00pm UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/15 "2021-09-17T18:00:20Z")
> [@jipolanco](#):
>
> I think the issue with this solution is that the `@inbounds` is not propagated to the interior of the `ntuple` s. Otherwise, it should be exactly equivalent to my solution above.
Yeah, that is a good point. Thank you. I have thought about `@inbounds` but about it inside `ntuple`.
> [@N5N3](#):
>
> `for d in ntuple(identity, N)`
This made my day. I was unaware that `1:N` and `ntuple(identity, N) can have different performance in a for-loop. Is it always the case, because that would make it very sensitive with respect to for loops in all such similar cases.
> [@mcabbott](#):
>
> The fact that `bc` makes the indexing unpredictable is also about a factor of 2 effect, for me. Can you duplicate the first row/col/etc of `Λ` as the last+1? Otherwise, you can deal with the edge afterwards. Not hard to write for `N=2` , more complicated for N=3… in general you could try [TiledIteration.EdgeIterator](https://github.com/JuliaArrays/TiledIteration.jl#edgeiterator) although I haven’t timed it here.
> [@stevengj](#):
>
> Note that this is a perfect case for using [ghost cells](https://discourse.julialang.org/t/finite-difference-laplacian-with-five-point-stencil/25014/2). Just make your array `Λ` bigger by one element in each direction, and copy the periodic boundaries before doing any “stencil” calculations.
Yes, there are even more options like precalculate the boundary indices and so on. But here I would like to have the most straight-forward but fast implementation.
> [@stevengj](#):
>
> So, even if you count the cost of copying over the periodic boundaries (which only needs to be done once whenever you update the array) it is more than 2× faster. If you don’t count the time for `periodicboundaries!` , then it is almost 4× faster.
>
> PS. Note that I eliminated the generator completely, compared to @N5N3’s solution.
Thanks, that sounds great. I will inspect it more soon.
Thank all the answers. I think it answers all my questions by now.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 20, 2021, 8:49am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/16 "2021-09-20T08:49:16Z")
Im principle it is a great thing. However, for a LxLxL = V lattice I usualy have V updates of individul sites before calculating the energy. So, it is still advisable to not copy the updates. But thanks for the great solution at this point.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 20, 2021, 8:50am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/17 "2021-09-20T08:50:17Z")
Unfortunately, it is still slower on my machine… ~ 370 μs
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 20, 2021, 9:12am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/18 "2021-09-20T09:12:42Z")
Ok, so the culpit here was the`sum` function. I achieve the same performance of 250 μs with this one, which uses a normal `for` loop
```julia
function energy_arbitrary(Λ::AbstractArray{T,N}) where {T,N}
E = 0
@simd for r ∈ CartesianIndices(Λ)
E += @inbounds Λ[r] * (+)((@inbounds Λ[ntuple(i -> i == d ? bc(r[i], size(Λ, i)) : r[i], N)...] for d ∈ 1:N)...)
end
return -E
end
```
---
### Author:  [@N5N3](https://discourse.julialang.org/u/N5N3)
#### Post date: [September 20, 2021, 9:43am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/19 "2021-09-20T09:43:19Z")
Well, I looked into the generated LLVM IR.
Both `sum(x)` (my version) and `+(x...)` (your version) is unrolled by the optimizer.
The difference is: `sum(x::Generator)` call `mapreduce(identity, add_sum, x)` to avoid overflow, while `+(x...)` not.
If you replace `sum(x)` with `reduce(+,x)`, the performance should be exactly same.
---
### Author:  [@stakaz](https://discourse.julialang.org/u/stakaz)
#### Post date: [September 20, 2021, 9:54am UTC](https://discourse.julialang.org/t/speed-up-multidimensional-code-without-base-cartesian/68328/20 "2021-09-20T09:54:20Z")
Ah, thanks. That sounds intuitive 😉