# Could anyone make this run faster?
**URL:** https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145
**Category:** Performance
**Created:** [December 26, 2022, 1:13pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145 "2022-12-26T13:13:01Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 1:13pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/1 "2022-12-26T13:13:01Z")
Hello!
I apologize for using unicode characters before hand, but I really like it π I have a function given as:
```julia
function βα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
# Skip distances outside the support of the kernel:
if q < 0.0 || q > 2.0
return SVector(0.0,0.0,0.0)
end
gradWx = Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (xα΅’β±Ό[1] / (q*h+1e-6))
gradWy = Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (xα΅’β±Ό[2] / (q*h+1e-6))
gradWz = Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (xα΅’β±Ό[3] / (q*h+1e-6))
return SVector(gradWx,gradWy,gradWz)
end
```
Which by using BenchmarkTools has a performance as such:
```julia
Ξ±D = 174
q = 1
h = 0.5659
xα΅’β±Ό = rand(SVector{3,Float64})
@benchmark βα΅’Wα΅’β±Ό($Ξ±D,$q,$xα΅’β±Ό,$h)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min β¦ max): 6.100 ns β¦ 115.700 ns β GC (min β¦ max): 0.00% β¦ 0.00%
Time (median): 6.300 ns β GC (median): 0.00%
Time (mean Β± Ο): 6.382 ns Β± 1.985 ns β GC (mean Β± Ο): 0.00% Β± 0.00%
β
ββββββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββ β
6.1 ns Histogram: frequency by time 7.1 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
And 6.1 ns is really fast, but any improvement can be of use, since this is a function which I have to run extremely often. I made some improvements to it by combining everything into a factor:
```julia
function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
# Skip distances outside the support of the kernel:
if q < 0.0 || q > 2.0
return SVector(0.0,0.0,0.0)
end
Fac = Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (1/(q*h+1e-6))
gradWx = Fac * xα΅’β±Ό[1]
gradWy = Fac * xα΅’β±Ό[2]
gradWz = Fac * xα΅’β±Ό[3]
return SVector(gradWx,gradWy,gradWz)
end
```
Which produces an improvement of ~8%:
```julia
@benchmark Optimβα΅’Wα΅’β±Ό($Ξ±D,$q,$xα΅’β±Ό,$h)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min β¦ max): 5.600 ns β¦ 149.000 ns β GC (min β¦ max): 0.00% β¦ 0.00%
Time (median): 5.700 ns β GC (median): 0.00%
Time (mean Β± Ο): 5.967 ns Β± 2.516 ns β GC (mean Β± Ο): 0.00% Β± 0.00%
ββ
β
ββββ
ββ
β
βββββββββββββββββββ
βββ
βββββββββββββββββββββββββββββ
β β
5.6 ns Histogram: log(frequency) by time 20.1 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
Iβve seen some of you being able to speed things up in an (to me) incredible way, so I just felt it would be stupid of me not to ask, if someone would know how to make this calculation faster π
Kind regards
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [December 26, 2022, 1:32pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/2 "2022-12-26T13:32:11Z")
> [@Ahmed\_Salih](#):
>
> `@benchmark Optimβα΅’Wα΅’β±Ό($Ξ±D,$q,$xα΅’β±Ό,$h)`
You can add `@fastmath`. Also note that you can just multiply `Fac` by `x`, to simplify the code:
```julia
julia> function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
# Skip distances outside the support of the kernel:
if q < 0.0 || q > 2.0
return SVector(0.0,0.0,0.0)
end
Fac = @fastmath Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (1/(q*h+1e-6))
return Fac * xα΅’β±Ό
end
Optimβα΅’Wα΅’β±Ό (generic function with 1 method)
```
Thatβs sligthly faster here.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 26, 2022, 1:32pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/3 "2022-12-26T13:32:45Z")
> [@Ahmed\_Salih](#):
>
> ```julia
> function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
> # Skip distances outside the support of the kernel:
> if q < 0.0 || q > 2.0
> return SVector(0.0,0.0,0.0)
> end
>
> Fac = Ξ±D * 1/h * (5*(q-2)^3*q)/8 * (1/(q*h+1e-6))
>
> gradWx = Fac * xα΅’β±Ό[1]
> gradWy = Fac * xα΅’β±Ό[2]
> gradWz = Fac * xα΅’β±Ό[3]
>
> return SVector(gradWx,gradWy,gradWz)
> end
>
> ```
Probably not faster, but more elegant and general:
```julia
function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
# Skip distances outside the support of the kernel:
if 0 <= q <= 2
Fac = Ξ±D*5*(q-2)^3*q / (8h*(q*h+1e-6))
else
Fac = 0.0 # or return zero(xα΅’β±Ό)
end
return Fac .* xα΅’β±Ό
end
```
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [December 26, 2022, 2:04pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/4 "2022-12-26T14:04:49Z")
I would guess you are approaching the theoretical maximum performance for running this one function, especially after `@fastmath`. Itβs only 2.5ns on my machine.
So most of the remaining gains will be in how you can parallelize it in the calling function.
For example running it on a GPU can be easily 30x faster for a large array:
```julia
julia> using CUDA
julia> xα΅’β±Ό = [rand(SVector{3,Float64}) for i in 1:100_000_000];
julia> cu_xα΅’β±Ό = CuArray(xα΅’β±Ό);
julia> @benchmark Optimβα΅’Wα΅’β±Ό.($Ξ±D,$q,$xα΅’β±Ό,$h)
BenchmarkTools.Trial: 6 samples with 1 evaluation.
Range (min β¦ max): 808.936 ms β¦ 1.161 s β GC (min β¦ max): 0.17% β¦ 17.53%
Time (median): 861.737 ms β GC (median): 9.89%
Time (mean Β± Ο): 910.588 ms Β± 131.796 ms β GC (mean Β± Ο): 10.66% Β± 8.91%
β β ββ β β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
809 ms Histogram: frequency by time 1.16 s <
Memory estimate: 2.24 GiB, allocs estimate: 2.
julia> @benchmark CUDA.@sync Optimβα΅’Wα΅’β±Ό.($Ξ±D,$q,$cu_xα΅’β±Ό,$h)
BenchmarkTools.Trial: 160 samples with 1 evaluation.
Range (min β¦ max): 21.589 ms β¦ 33.286 ms β GC (min β¦ max): 0.00% β¦ 8.57%
Time (median): 31.562 ms β GC (median): 8.11%
Time (mean Β± Ο): 31.223 ms Β± 1.467 ms β GC (mean Β± Ο): 8.03% Β± 1.18%
ββ ββββ
ββββββββββββββββββββββββββββββββββββββββββββββββ
ββ
βββββββ
ββ β
21.6 ms Histogram: frequency by time 33.1 ms <
Memory estimate: 4.80 KiB, allocs estimate: 85.
```
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 2:05pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/5 "2022-12-26T14:05:22Z")
Thank you both @lmiq and @DNF !
I tried the @fastmath macro, but it did not seem to have any effect on my system, so I removed it again.
I gave your version, DNF, a try, and actually yours is much faster π !
```julia
function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
# Skip distances outside the support of the kernel:
if 0 < q < 2
Fac = Ξ±D*5*(q-2)^3*q / (8h*(q*h+1e-6))
else
Fac = 0.0 # or return zero(xα΅’β±Ό)
end
return Fac .* xα΅’β±Ό
end
@benchmark Optimβα΅’Wα΅’β±Ό($Ξ±D,$q,$xα΅’β±Ό,$h)
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
Range (min β¦ max): 2.500 ns β¦ 31.300 ns β GC (min β¦ max): 0.00% β¦ 0.00%
Time (median): 2.700 ns β GC (median): 0.00%
Time (mean Β± Ο): 2.863 ns Β± 1.398 ns β GC (mean Β± Ο): 0.00% Β± 0.00%
β
ββββββββββββββββββββ
ββββββββββββββββββββββββββββββββββββββ β
2.5 ns Histogram: frequency by time 3.4 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
```
This is about a 50% speed up compared to base-line. I am not sure why this difference is this big, but I assume the way you multiply constants together makes it easier for the compiler perhaps?
Kind regards
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 2:06pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/6 "2022-12-26T14:06:49Z")
Thank you for this! π
One of my aims is also running my particle simulation code (which this function above is a part of) on the GPU, but part of that journey is making sure I understand fully how a CPU versions works.
I have already tried GPU, but I wrote it horribly so it ended up 10x slower than CPU. Will make sure to ask when I am ready for that again π
Kind regards
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [December 26, 2022, 2:24pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/7 "2022-12-26T14:24:39Z")
This version will still be much slower for small loads on a GPU - you need a large task to see the kind of gains I posted.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [December 26, 2022, 2:25pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/8 "2022-12-26T14:25:26Z")
> [@Ahmed\_Salih](#):
>
> but I assume the way you multiply constants together makes it easier for the compiler perhaps?
Probably thatβs because less divisions are performed. Divisions are expensive.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 26, 2022, 2:26pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/9 "2022-12-26T14:26:14Z")
> [@Ahmed\_Salih](#):
>
> but I assume the way you multiply constants together makes it easier for the compiler perhaps?
Perhaps itβs that mine has just one division instead of three. And the `@fastmath` macro perhaps gives the compiler permission to reorder the expression similarly.
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 2:49pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/10 "2022-12-26T14:49:20Z")
Great, thank you all π
As a final note, @DNF, using the form provided by you, my test case went from 30 to ~21-22 minutes, so roughly a 33% improvement, which was really nice!
Kind regards
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 26, 2022, 3:13pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/11 "2022-12-26T15:13:41Z")
In an effort to make calculation even faster, what can you tell us about the values of `h`? Are they different in every call, or do they represent a stable parameter or hyperparameter?
As for `q`, what is its meaning? When does it get very small i.e. \< 1e-5 ?
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 3:45pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/12 "2022-12-26T15:45:41Z")
h is in this case a constant parameter called the βSmoothing lengthβ which determines the βradius of influenceβ material points have on each other in Smoothed Particle Hydrodynamics.
The normalized distance is then given as q which is defined as: r/h, where r is the Euclidean distance between two points.
q is a parameter which goes into a kernel function to estimate βamount of influenceβ - the kernel is defined in such a way that it only has an influence on another material point if the distance is between 0 to 2, i.e. 0 \< q \< 2.
This condition is in a sense fulfilled by the neighbour searching algorithm, but as a safety precaution I like keeping it in.
q gets extremely small when the same material point is considered twice or two points are exceptionally close, and thereby have a large influence on each other.
Kind regards
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 26, 2022, 4:23pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/13 "2022-12-26T16:23:22Z")
How common is the case when `q` is _not_ between 0 and 2? If it is quite rare, you may be able to get further speedups by removing the branch and use Loopvectorization.jl
```julia
function Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
Fac = ifelse(0 < q < 2, Ξ±D*5*(q-2)^3*q / (8h*(q*h+1e-6)), 0.0)
return Fac .* xα΅’β±Ό
end
```
And then use something like (donβt know exactly what you want here)
```julia
@tturbo for k in ...
xα΅’β±Ό[k] = Optimβα΅’Wα΅’β±Ό(Ξ±D,q,xα΅’β±Ό,h)
end
```
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 26, 2022, 6:16pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/14 "2022-12-26T18:16:44Z")
> [@DNF](#):
>
> `0 < q < 2`
Actually, this may include a branch. Not sure.
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 6:19pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/15 "2022-12-26T18:19:23Z")
In theory it should βneverβ happen, since the neighbor list would throw it out if it was above 2, so only the condition where it is close to zero is realistically possible. Still I like to keep it in π
I tried to remove it and didnβt see any performance gain.
Kind regards
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [December 26, 2022, 6:21pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/16 "2022-12-26T18:21:35Z")
The gain would be when iterating over this in a hot loop. When all branches are removed it could be simd vectorized, with potentially huge performance improvements.
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 6:23pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/17 "2022-12-26T18:23:45Z")
Okay got it - will take note of that.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 26, 2022, 6:34pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/18 "2022-12-26T18:34:52Z")
There is another gain to squeeze and that is to remove the remaining division. The `q` factor would cancel if it wasnβt for the `1e-6` epsilon making sure the interaction zeroes at `q = 0`. But breaking the function at another point close to zero, allows to approximate with no division.
I didnβt really find a best approximation, but the final function should look something like this:
```julia
function Optimβ2(Ξ±D,q,invh)
Fac = ifelse(q < 2, ifelse(q < 2e-6, -Ξ±D*5*q*invh*1e6, Ξ±D*5*(q-2)^3*0.125*invh^2), 0.0)
return Fac
end
```
This is for testing, the real `Optim` scalar multiplies with `x_ij`. Also, I suppose `q` never really goes below zero (Euclidean distance). In any case, this version uses `invh` i.e. `1/h` which can be calculated outside the hot loop.
```julia
julia> lineplot(q->Optimβ2(Ξ±D,q,ih), 0.0,0.0001)
ββββββββββββββββββββββββββββββββββββββββββ
0 ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β #87(x)
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
f(x) ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’£β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’Έβ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’Έβ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’Έβ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ Έβ €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β €β
-3 000 ββ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββββββββββββββββββββββββββββββββββββββββββ
β 0β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β 0.0001β
β β β β β β β β β β β β β β β β β β β β β xβ β β β β β β β β β β β β β β β β β β β
julia> lineplot(q->Optimβ1(Ξ±D,q,h), 0.0,0.0001)
ββββββββββββββββββββββββββββββββββββββββββ
0 ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β #89(x)
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’£β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ’Έβ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
f(x) ββ’Έβ‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β’±β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β β’£β‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β β β β’β‘β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β β β β β β β β ¦β €β €β£β£β£β£β£β β β β β β β β β β β β β β β β β β β β β β β β β
ββ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
-3 000 ββ β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β
ββββββββββββββββββββββββββββββββββββββββββ
β 0β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β β 0.0001β
β β β β β β β β β β β β β β β β β β β β β xβ β β β β β β β β β β β β β β β β β β β
```
This is a comparison chart. Note the harsher knee with the divisionless approximation. The softer knee can be recovered with a quadratic approximation for small `q`s (still no division). Anyway, this will not get a 10x increase for sure.
BTW if you are running on an M1 CPU Mac, the difference might be even smaller, as division seems to be quite fast on it.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [December 26, 2022, 7:12pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/19 "2022-12-26T19:12:23Z")
> [@Ahmed\_Salih](#):
>
> In theory it should βneverβ happen, since the neighbor list would throw it out if it was above 2, so only the condition where it is close to zero is realistically possible. Still I like to keep it in π
In terms of the actual particle simulation, there are two options: update, or not, the pair list at every step.
It is common to compute the list up to a distance greater than the actual cutoff, and do not upate it at every step. In that case you need to compute the distance again in that loop and check for the cutoff. This is the easiest strategy if you later want to port the computation of the forces to the GPU, for example. This is what Molly does.
Alternatively, you can update the list at every step. In that case you donβt need to check for the distance again.
The first is usually faster, but not always the particles move in a pace that can justify not updating the lists at a reasonable interval.
---
### Author:  [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)
#### Post date: [December 26, 2022, 7:33pm UTC](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145/20 "2022-12-26T19:33:00Z")
> [@lmiq](#):
>
> Molly
Yep, that is also a very valid approach! I cannot implement it yet though, since I do not know how to estimate when particles βroughlyβ have left a box.
I am not entirely sure how you handle it in the CellListMap package, but I have not noticed any significant delay in updating the neighbourlist at each time step (even at 1e-5 increments).
If you know of some material where to read about βwhen to updateβ the list, I would love to see it. Interesting to see how much there is in common with Molecular simulations π
Kind regards
[Next page](https://discourse.julialang.org/t/could-anyone-make-this-run-faster/92145.md?page=2)