# Can this be written even faster? (CPU)

**URL:** <https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924>\
**Category:** Performance\
**Created:** [February 8, 2024, 12:05pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924 "2024-02-08T12:05:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 12:05pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/1 "2024-02-08T12:05:21Z")

</div>

Hello!

My usual experience has been that when I think something is as fast as I think it can be, then someone shows me how to do it even faster 🙂

I’ve made a script attached below with a serial version and different ways of multi-threading. The results I get from **s** erial, **t** hreaded etc. are as follows:

```julia
serial: Trial(953.400 μs)
@threads: Trial(254.000 μs)
@threads :static :Trial(248.900 μs)
Polyester.@batch: Trial(239.600 μs)

```

I know Polyester is recommended when it works to reduce allocations, but I was wondering if there was an even faster way to calculate this. For me 240 μs is slow. I need to calculate this function a lot of times, so any reduction in time is worthwhile.

Kind regards

```julia
using BenchmarkTools
using Polyester
using StaticArrays
using LoopVectorization

# Generate the needed constants
N = 6000
Cb = 100000
xᵢⱼ = rand(SVector{3,Float64},N*4)
list = [(rand(1:N), rand(1:N), rand()) for _ in 1:N*4]
ρ₀ = 1000
γ⁻¹ = 1/7
g = 9.81
drhopLp = zeros(length(list))
drhopLn = zeros(length(list))

# In this code, use of multi-threading to calculate the heavy part for density diffusion to work
function serial_approach(Cb,xᵢⱼ,list,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    for iter in eachindex(list)
        xⱼᵢ = -xᵢⱼ[iter]
        Pᵢⱼᴴ = ρ₀ * (-g) * xⱼᵢ[2]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

# In this code, use of multi-threading to calculate the heavy part for density diffusion to work
function threads_approach(Cb,xᵢⱼ,list,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    Base.Threads.@threads for iter in eachindex(list)
        xⱼᵢ = -xᵢⱼ[iter]
        Pᵢⱼᴴ = ρ₀ * (-g) * xⱼᵢ[2]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

function threads_static_approach(Cb,xᵢⱼ,list,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    Base.Threads.@threads :static for iter in eachindex(list)
        xⱼᵢ = -xᵢⱼ[iter]
        Pᵢⱼᴴ = ρ₀ * (-g) * xⱼᵢ[2]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

function polyester_approach(Cb,xᵢⱼ,list,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @batch for iter in eachindex(list)
        xⱼᵢ = -xᵢⱼ[iter]
        Pᵢⱼᴴ = ρ₀ * (-g) * xⱼᵢ[2]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

function loopvectorization_approach(Cb,xᵢⱼ,list,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @tturbo for iter in eachindex(list)
        xⱼᵢ = -xᵢⱼ[iter]
        Pᵢⱼᴴ = ρ₀ * (-g) * xⱼᵢ[2]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

s = @benchmark serial_approach($Cb,$xᵢⱼ,$list,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);
t = @benchmark threads_approach($Cb,$xᵢⱼ,$list,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);
ts = @benchmark threads_static_approach($Cb,$xᵢⱼ,$list,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);
p = @benchmark polyester_approach($Cb,$xᵢⱼ,$list,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);
# why @tturbo no work?
#l = @benchmark loopvectorization_approach($Cb,$xᵢⱼ,$list,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);

println(s,t,ts,p)

```

---

<div class="post-metadata">

**Author:** ![linwaytin](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@linwaytin](https://discourse.julialang.org/u/linwaytin)\
**Post date:** [February 8, 2024, 2:32pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/2 "2024-02-08T14:32:34Z")

</div>

Maybe try to put `const` in front of those global variables?

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [February 8, 2024, 2:51pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/3 "2024-02-08T14:51:32Z")

</div>

It appears that he is passing all the globals to the functions, so the functions do not access the globals and wouldn’t incur the performance degradation due to such access, but your suggestion is good practice anyway.

---

<div class="post-metadata">

**Author:** ![linwaytin](https://avatars.discourse-cdn.com/v4/letter/l/898d66/32.png) [@linwaytin](https://discourse.julialang.org/u/linwaytin)\
**Post date:** [February 8, 2024, 2:53pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/4 "2024-02-08T14:53:04Z")

</div>

You’re right. Sorry I didn’t notice that.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 8, 2024, 2:56pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/5 "2024-02-08T14:56:41Z")

</div>

Is Float32 precision sufficient? If so, that will be a free 2-3x boost.

---

<div class="post-metadata">

**Author:** ![Jeff\_Emanuel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeff_emanuel/32/15440_2.png) [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)\
**Post date:** [February 8, 2024, 2:57pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/6 "2024-02-08T14:57:45Z")

</div>

Your seventh roots seem to be the only operation that could take long. Do you know anything a priori about the arguments? Could you tabulate the values? Could you use an approximation?

Maybe you are just hitting memory bandwidth limits and you need to rethink approaching the problem at a higher level so that you can do more calculation on smaller parts of the array. Oscar’s suggestion addresses memory bandwidth by decreasing the size of the elements.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 3:09pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/7 "2024-02-08T15:09:36Z")

</div>

> [@Oscar\_Smith](#):
>
> Is Float32 precision sufficient? If so, that will be a free 2-3x boost.

I would prefer to stay in Float64 for now 😊

> [@Jeff\_Emanuel](#):
>
> Your seventh roots seem to be the only operation that could take long. Do you know anything a priori about the arguments? Could you tabulate the values? Could you use an approximation?

The constants Cb, rho0 and gamma I know before hand. Unfortunately it is not feasibile with an approximation afaik - atleast want to avoid that for now.

I’ve been trying to get a `@tturbo` version to work but unfortunately it allocates for some reason:

```julia
using BenchmarkTools
using Polyester
using StaticArrays
using LoopVectorization

# Generate the needed constants
N = 6000
Cb = 100000
xᵢⱼᶻ = rand(Float64,N*4)
xᵢⱼ = rand(SVector{3,Float64},N*4)
list = [(rand(1:N), rand(1:N), rand()) for _ in 1:N*4]
ρ₀ = 1000
γ⁻¹ = 1/7
g = 9.81
drhopLp = zeros(length(list))
drhopLn = zeros(length(list))

function loopvectorization_approach(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @tturbo for iter in eachindex(xᵢⱼᶻ)
        # IF I COMMENT '* -xᵢⱼᶻ[iter]' NO ALLOCATIONS, WHY?
        Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼᶻ[iter]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

l = @benchmark loopvectorization_approach($Cb,$xᵢⱼᶻ,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);

```

I’ve introduced a new vector since I couldn’t double index into StaticArray using LoopVectorization. When I comment out the part mentioned in green then no allocs, but if I keep it in a lot of allocs:

```julia
BenchmarkTools.Trial: 1371 samples with 1 evaluation.
 Range (min … max): 1.542 ms … 64.285 ms ┊ GC (min … max): 0.00% … 41.45%
 Time (median): 2.444 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 3.648 ms ± 5.515 ms ┊ GC (mean ± σ): 29.74% ± 18.84%

  ▆██▅▃▃
  ██████▆▅▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▆██▅▆▆▆▆▅ █
  1.54 ms Histogram: log(frequency) by time 26.8 ms <

 Memory estimate: 8.24 MiB, allocs estimate: 132000.

BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max): 2.467 μs … 11.067 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 2.600 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 2.986 μs ± 1.229 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▆█▆▂ ▁▁▂▁ ▄ ▂
  ██████████▇▇█▇▇▆▇▄▅▅▄▄▄▄▁▃▁▁▄▃▃▄▁▃▄▁▃▃▁▄▁▁▃▁▁▁▁▁▁▄▃▁▁▁▁▇▄█ █
  2.47 μs Histogram: log(frequency) by time 7.18 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

Would anyone happen to know how to deal with that?

Kind regards

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [February 8, 2024, 3:47pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/8 "2024-02-08T15:47:09Z")

</div>

Could the allocations you’re seeing be related to [this?](https://github.com/JuliaSIMD/LoopVectorization.jl/issues/526#issuecomment-1898969405)

you can try converting `ρ₀ = 1000` to a float before you compute the inner sum. I think the Julia compiler seems to give up more easily on inferring types in Julia 1.10

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [February 8, 2024, 3:50pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/9 "2024-02-08T15:50:31Z")

</div>

That will still infer fine. That said, it will be faster to work with `inv(ρ₀)` since multiplication is faster than division.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 3:52pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/10 "2024-02-08T15:52:37Z")

</div>

> [@Salmon](#):
>
> Could the allocations you’re seeing be related to [this?](https://github.com/JuliaSIMD/LoopVectorization.jl/issues/526#issuecomment-1898969405)
> 
> you can try converting `ρ₀ = 1000` to a float before you compute the inner sum. I think the Julia compiler seems to give up more easily on inferring types in Julia 1.10

Updated the code to strictly use Float64:

```julia
using BenchmarkTools
using StaticArrays
using LoopVectorization

# Generate the needed constants
N = 6000
Cb = Float64(100000)
xᵢⱼᶻ = rand(Float64,N*4)
xᵢⱼ = rand(SVector{3,Float64},N*4)
ρ₀ = Float64(1000)
γ⁻¹ = Float64(1/7)
g = Float64(9.81)
drhopLp = zeros(Float64, length(list))
drhopLn = zeros(Float64, length(list))

function loopvectorization_approach(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @tturbo for iter in eachindex(xᵢⱼᶻ)
        # IF I COMMENT '* -xᵢⱼᶻ[iter]' NO ALLOCATIONS, WHY?
        Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼᶻ[iter]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ/Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ/Cb), γ⁻¹) - 1)
        
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

l = @benchmark loopvectorization_approach($Cb,$xᵢⱼᶻ,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);

display(l)

BenchmarkTools.Trial: 1246 samples with 1 evaluation.
 Range (min … max): 1.339 ms … 59.475 ms ┊ GC (min … max): 0.00% … 36.11%
 Time (median): 2.644 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.011 ms ± 5.803 ms ┊ GC (mean ± σ): 28.21% ± 18.78%

  ▂█▇▇▄▂
  ██████▁▄▁▄▁▄▁▁▁▁▅▁▁▄▁▁▁▁▄▁▁▁▁▁▁▁▁▁▄▁▁▄▅▇▇▆▆▇▄▅▄▄▅▁▄▁▁▄▄▄▄▅ █
  1.34 ms Histogram: log(frequency) by time 32.9 ms <

 Memory estimate: 8.24 MiB, allocs estimate: 132000.

```

Issue with allocations persist.

> [@Oscar\_Smith](#):
>
> That will still infer fine. That said, it will be faster to work with `inv(ρ₀)` since multiplication is faster than division.

Sorry I didn’t get this, you meant Cb right?

Kind regards

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 8, 2024, 3:53pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/11 "2024-02-08T15:53:38Z")

</div>

LoopVectorization often allocates on Julia 1.10.  
Can you try 1.9?

I haven’t looked at this.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 3:54pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/12 "2024-02-08T15:54:14Z")

</div>

The difference between inv(Cb) and Cb is 0.040 ms, nice!

Thank you 😊 In my case this does add up.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 4:01pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/13 "2024-02-08T16:01:34Z")

</div>

@Elrod, I would never have thought of this myself 😅

Using Julia 1.9.4:

```julia
using BenchmarkTools
using StaticArrays
using LoopVectorization

# Generate the needed constants
N = 6000
Cb = inv(Float64(100000))
xᵢⱼᶻ = rand(Float64,N*4)
xᵢⱼ = rand(SVector{3,Float64},N*4)
ρ₀ = Float64(1000)
γ⁻¹ = Float64(1/7)
g = Float64(9.81)
drhopLp = zeros(Float64, length(xᵢⱼ))
drhopLn = zeros(Float64, length(xᵢⱼ))

function loopvectorization_approach(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @tturbo for iter in eachindex(xᵢⱼᶻ)
        # IF I COMMENT '* -xᵢⱼᶻ[iter]' NO ALLOCATIONS, WHY?
        Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼᶻ[iter]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ * Cb), γ⁻¹) - 1)
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ * Cb), γ⁻¹) - 1)
        
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

l = @benchmark loopvectorization_approach($Cb,$xᵢⱼᶻ,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn);

display(l)

```

```julia
julia> display(l)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 163.000 μs … 556.700 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 172.900 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 178.892 μs ± 27.539 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

    ▆▆▂█▅▄▇▄▃▂▂▃▁▃▁ ▁▁▃ ▂▂ ▂
  ▆▁███████████████▇███████▆▆▆▆▆▇▆▅▆▅▅▅▆▅▆▃▄▄▄▅▃▄▃▆▄▄▅▇▆▆▆▅▆▅▆▅ █
  163 μs Histogram: log(frequency) by time 266 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

Would love if this could work on version 1.10, but I suppose my best option for now is to go back to 1.9? 😕

EDIT: The allocations kill my performance since I am doing particle simulation stuff.

Kind regards

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [February 8, 2024, 4:14pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/14 "2024-02-08T16:14:13Z")

</div>

> [@Ahmed\_Salih](#):
>
> Would love if this could work on version 1.10, but I suppose my best option for now is to go back to 1.9? 😕
> 
> EDIT: The allocations kill my performance since I am doing particle simulation stuff.

Can you verify with code\_warntype whether the allocations come from a type instability in the LoopVectorization version?  
If yes, you might still make it work in 1.10 by converting some the other integers to floats explicitly.  
It’s a bit unjulian, but if it works…

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 4:30pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/15 "2024-02-08T16:30:55Z")

</div>

code\_warntype

```julia
MethodInstance for loopvectorization_approach(::Float64, ::Vector{Float64}, ::Float64, ::Float64, ::Float64, ::Vector{Float64}, ::Vector{Float64})
  from loopvectorization_approach(Cb, xᵢⱼᶻ, ρ₀, γ⁻¹, g, drhopLp, drhopLn) @ Main REPL[33]:1
Arguments
  #self#::Core.Const(loopvectorization_approach)
  Cb::Float64
  xᵢⱼᶻ::Vector{Float64}
  ρ₀::Float64
  γ⁻¹::Float64
  g::Float64
  drhopLp::Vector{Float64}
  drhopLn::Vector{Float64}
Locals
  @_9::Union{}
  val@_10::Union{}
  @_11::Int64
  @_12::Int64
  @_13::Int64
  val@_14::Base.OneTo{Int64}
  Wvecwidth##::Static.StaticInt{4}
  Tloopeltype##::Type{Float64}
  vargsym#233::Tuple{Tuple{CloseOpenIntervals.CloseOpen{Static.StaticInt{0}, Int64}}, Tuple{LayoutPointers.GroupedStridedPointers{Tuple{Ptr{Float64}, Ptr{Float64}, Ptr{Float64}}, (1, 1, 1), (0, 0, 0), ((1,), (1,), (1,)), ((1,), (2,), (3,)), Tuple{Static.StaticInt{8}, Static.StaticInt{8}, Static.StaticInt{8}}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}, Static.StaticInt{0}}}, Vararg{Float64, 4}}}
  ##grouped#strided#pointer####19###::LayoutPointers.GroupedStridedPointers{Tuple{Ptr{Float64}, Ptr{Float64}, Ptr{Float64}}, (1, 1, 1), (0, 0, 0), ((1,), (1,), (1,)), ((1,), (2,), (3,)), Tuple{Static.StaticInt{8}, Static.StaticInt{8}, Static.StaticInt{8}}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}, Static.StaticInt{0}}}
  preserve#buffer#@_19::Vector{Float64}
  vptr##_drhopLn::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
  preserve#buffer#@_21::Vector{Float64}
  vptr##_drhopLp::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
  #loopconstnumber###12###::Int64
  preserve#buffer#@_24::Vector{Float64}
  vptr##_xᵢⱼᶻ::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
  #temp###5###::Float64
  #temp###6###::Float64
  #iter_loop_upper_bound###4###::Int64
  #iter_loop_lower_bound###3###::Static.StaticInt{1}
  #loopleniter###2###::Int64
  #looprangeiter###1###::Static.OptionallyStaticUnitRange{Static.StaticInt{1}, Int64}
  msg::Union{}
  kwargs::Union{}
  line::Union{}
  file::Union{}
  id::Union{}
  logger::Union{}
  _module::Union{}
  group::Union{}
  std_level::Union{}
  level::Union{}
  err::Union{}
  iter::Union{}
  ρⱼᵢᴴ::Union{}
  Pⱼᵢᴴ::Union{}
  ρᵢⱼᴴ::Union{}
  Pᵢⱼᴴ::Union{}
  @_48::Union{}
Body::Nothing
1 ── Core.NewvarNode(:(@_9))
│ Core.NewvarNode(:(val@_10))
│ Core.NewvarNode(:(@_11))
│ Core.NewvarNode(:(@_12))
│ Core.NewvarNode(:(@_13))
│ Core.NewvarNode(:(Wvecwidth##))
│ Core.NewvarNode(:(Tloopeltype##))
│ Core.NewvarNode(:(vargsym#233))
│ Core.NewvarNode(:(##grouped#strided#pointer####19###))
│ Core.NewvarNode(:(preserve#buffer#@_19))
│ Core.NewvarNode(:(vptr##_drhopLn))
│ Core.NewvarNode(:(preserve#buffer#@_21))
│ Core.NewvarNode(:(vptr##_drhopLp))
│ Core.NewvarNode(:(#loopconstnumber###12###))
│ Core.NewvarNode(:(preserve#buffer#@_24))
│ Core.NewvarNode(:(vptr##_xᵢⱼᶻ))
│ Core.NewvarNode(:(#temp###5###))
│ Core.NewvarNode(:(#temp###6###))
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ Main.nothing
│ nothing
│ (val@_14 = Main.eachindex(xᵢⱼᶻ))
│ nothing
│ %30 = val@_14::Base.OneTo{Int64}
│ (#looprangeiter###1### = LoopVectorization.canonicalize_range(%30))
│ (#loopleniter###2### = StaticArrayInterface.static_length(#looprangeiter###1###))
│ (#iter_loop_lower_bound###3### = LoopVectorization.maybestaticfirst(#looprangeiter###1###))
│ (#iter_loop_upper_bound###4### = LoopVectorization.maybestaticlast(#looprangeiter###1###))
│ %35 = LoopVectorization.check_args(xᵢⱼᶻ, drhopLp, drhopLn)::Core.Const(true)
└─── goto #15 if not %35
2 ── goto #15 if not true
3 ── %38 = (LoopVectorization.can_turbo)(LoopVectorization.sub_fast, Val{1}())::Core.Const(true)
└─── goto #15 if not %38
4 ── %40 = (LoopVectorization.can_turbo)(LoopVectorization.mul_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %40
5 ── %42 = (LoopVectorization.can_turbo)(LoopVectorization.vfmadd_fast, Val{3}())::Core.Const(true)
└─── goto #15 if not %42
6 ── %44 = (LoopVectorization.can_turbo)(LoopVectorization.pow_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %44
7 ── %46 = (LoopVectorization.can_turbo)(LoopVectorization.sub_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %46
8 ── %48 = (LoopVectorization.can_turbo)(LoopVectorization.mul_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %48
9 ── %50 = (LoopVectorization.can_turbo)(LoopVectorization.sub_fast, Val{1}())::Core.Const(true)
└─── goto #15 if not %50
10 ─ %52 = (LoopVectorization.can_turbo)(LoopVectorization.vfmadd_fast, Val{3}())::Core.Const(true)
└─── goto #15 if not %52
11 ─ %54 = (LoopVectorization.can_turbo)(LoopVectorization.pow_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %54
12 ─ %56 = (LoopVectorization.can_turbo)(LoopVectorization.sub_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %56
13 ─ %58 = (LoopVectorization.can_turbo)(LoopVectorization.mul_fast, Val{2}())::Core.Const(true)
└─── goto #15 if not %58
14 ─ (#temp###6### = LoopVectorization.sub_fast(g))
│ (#temp###5### = LoopVectorization.mul_fast(ρ₀, #temp###6###))
│ %62 = LoopVectorization.stridedpointer_preserve(xᵢⱼᶻ)::Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Vector{Float64}}
│ %63 = Base.indexed_iterate(%62, 1)::Core.PartialStruct(Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Int64}, Any[LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Core.Const(2)])
│ (vptr##_xᵢⱼᶻ = Core.getfield(%63, 1))
│ (@_13 = Core.getfield(%63, 2))
│ %66 = Base.indexed_iterate(%62, 2, @_13::Core.Const(2))::Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(3)])
│ (preserve#buffer#@_24 = Core.getfield(%66, 1))
│ (#loopconstnumber###12### = 1)
│ %69 = LoopVectorization.stridedpointer_preserve(drhopLp)::Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Vector{Float64}}
│ %70 = Base.indexed_iterate(%69, 1)::Core.PartialStruct(Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Int64}, Any[LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Core.Const(2)])
│ (vptr##_drhopLp = Core.getfield(%70, 1))
│ (@_12 = Core.getfield(%70, 2))
│ %73 = Base.indexed_iterate(%69, 2, @_12::Core.Const(2))::Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(3)])
│ (preserve#buffer#@_21 = Core.getfield(%73, 1))
│ %75 = LoopVectorization.stridedpointer_preserve(drhopLn)::Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Vector{Float64}}
│ %76 = Base.indexed_iterate(%75, 1)::Core.PartialStruct(Tuple{LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Int64}, Any[LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}, Core.Const(2)])
│ (vptr##_drhopLn = Core.getfield(%76, 1))
│ (@_11 = Core.getfield(%76, 2))
│ %79 = Base.indexed_iterate(%75, 2, @_11::Core.Const(2))::Core.PartialStruct(Tuple{Vector{Float64}, Int64}, Any[Vector{Float64}, Core.Const(3)])
│ (preserve#buffer#@_19 = Core.getfield(%79, 1))
│ %81 = vptr##_xᵢⱼᶻ::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
│ %82 = Core.tuple(#iter_loop_lower_bound###3###)::Core.Const((static(1),))
│ %83 = LoopVectorization.gespf1(%81, %82)::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}
│ %84 = LoopVectorization.densewrapper(%83, xᵢⱼᶻ)::LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}
│ %85 = vptr##_drhopLp::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
│ %86 = Core.tuple(#iter_loop_lower_bound###3###)::Core.Const((static(1),))
│ %87 = LoopVectorization.gespf1(%85, %86)::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}
│ %88 = LoopVectorization.densewrapper(%87, drhopLp)::LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}
│ %89 = vptr##_drhopLn::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{1}}}
│ %90 = Core.tuple(#iter_loop_lower_bound###3###)::Core.Const((static(1),))
│ %91 = LoopVectorization.gespf1(%89, %90)::LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}
│ %92 = LoopVectorization.densewrapper(%91, drhopLn)::LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}
│ %93 = Core.tuple(%84, %88, %92)::Tuple{LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}, LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}, LayoutPointers.DensePointerWrapper{(true,), Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}, LayoutPointers.StridedPointer{Float64, 1, 1, 0, (1,), Tuple{Static.StaticInt{8}}, Tuple{Static.StaticInt{0}}}}}
│ %94 = Main.Val::Core.Const(Val)
│ %95 = ()::Core.Const(())
│ %96 = Core.apply_type(%94, %95)::Core.Const(Val{()})
│ %97 = (%96)()::Core.Const(Val{()}())
│ %98 = LoopVectorization.grouped_strided_pointer(%93, %97)::Tuple{LayoutPointers.GroupedStridedPointers{Tuple{Ptr{Float64}, Ptr{Float64}, Ptr{Float64}}, (1, 1, 1), (0, 0, 0), ((1,), (1,), (1,)), ((1,), (2,), (3,)), Tuple{Static.StaticInt{8}, Static.StaticInt{8}, Static.StaticInt{8}}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}, Static.StaticInt{0}}}, Tuple{Nothing, Nothing, Nothing}}
│ (##grouped#strided#pointer####19### = (getfield)(%98, 1))
│ %100 = $(Expr(:gc_preserve_begin, :(preserve#buffer#@_24), :(preserve#buffer#@_21), :(preserve#buffer#@_19)))
│ %101 = LoopVectorization.zerorangestart(#looprangeiter###1###)::CloseOpenIntervals.CloseOpen{Static.StaticInt{0}, Int64}
│ %102 = Core.tuple(%101)::Tuple{CloseOpenIntervals.CloseOpen{Static.StaticInt{0}, Int64}}
│ %103 = Core.tuple(##grouped#strided#pointer####19###, ρ₀, #temp###5###, Cb, γ⁻¹)::Tuple{LayoutPointers.GroupedStridedPointers{Tuple{Ptr{Float64}, Ptr{Float64}, Ptr{Float64}}, (1, 1, 1), (0, 0, 0), ((1,), (1,), (1,)), ((1,), (2,), (3,)), Tuple{Static.StaticInt{8}, Static.StaticInt{8}, Static.StaticInt{8}}, Tuple{Static.StaticInt{0}, Static.StaticInt{0}, Static.StaticInt{0}}}, Vararg{Float64, 4}}
│ (vargsym#233 = Core.tuple(%102, %103))
│ %105 = LoopVectorization.eltype(xᵢⱼᶻ)::Core.Const(Float64)
│ %106 = LoopVectorization.eltype(drhopLp)::Core.Const(Float64)
│ %107 = LoopVectorization.eltype(drhopLn)::Core.Const(Float64)
│ (Tloopeltype## = LoopVectorization.promote_type(%105, %106, %107))
│ (Wvecwidth## = LoopVectorization.pick_vector_width(Tloopeltype##::Core.Const(Float64)))
│ %110 = LoopVectorization._turbo_!::Core.Const(LoopVectorization._turbo_!)
│ %111 = Core.apply_type(Main.Val, (false, 0, 0, 0, false, 0xffffffffffffffff, 1, true))::Core.Const(Val{(false, 0, 0, 0, false, 0xffffffffffffffff, 1, true)})
│ %112 = (%111)()::Core.Const(Val{(false, 0, 0, 0, false, 0xffffffffffffffff, 1, true)}())
│ %113 = LoopVectorization.avx_config_val(%112, Wvecwidth##)::Core.Const(Val{(false, 0, 0, 0, false, 4, 32, 15, 64, 0x0000000000000004, 1, true)}())
│ %114 = Main.Val::Core.Const(Val)
│ %115 = Core.tuple(:LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0001, 0x0000), Symbol("##DROPPED#CONSTANT##"), Symbol("##DROPPED#CONSTANT##"), LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0002, 0x0000), Symbol("##DROPPED#CONSTANT##"), Symbol("##DROPPED#CONSTANT##"), LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0003, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0004, 0x0000), :LoopVectorization, :getindex, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.memload, 0x0005, 0x0001), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000005, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0006, 0x0000), :LoopVectorization, :mul_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000040006, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0007, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0008, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0009, 0x0000), :LoopVectorization, :vfmadd_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000700080009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000a, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x000b, 0x0000), :LoopVectorization, :pow_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x000000000000000000000000000a000b, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000c, 0x0000), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x000000000000000000000000000c0009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000d, 0x0000), :LoopVectorization, :mul_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000000001000d, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000e, 0x0000), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000007, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000f, 0x0000), :LoopVectorization, :vfmadd_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000f00080009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0010, 0x0000), :LoopVectorization, :pow_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000000010000b, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0011, 0x0000), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000110009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0012, 0x0000), :LoopVectorization, :mul_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000010012, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0013, 0x0000), :LoopVectorization, :setindex!, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000000000000e, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.memstore, 0x0014, 0x0002), :LoopVectorization, :setindex!, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000013, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.memstore, 0x0015, 0x0003))::Core.Const((:LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0001, 0x0000), Symbol("##DROPPED#CONSTANT##"), Symbol("##DROPPED#CONSTANT##"), LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0002, 0x0000), Symbol("##DROPPED#CONSTANT##"), Symbol("##DROPPED#CONSTANT##"), LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0003, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0004, 0x0000), :LoopVectorization, :getindex, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.memload, 0x0005, 0x0001), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000005, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0006, 0x0000), :LoopVectorization, :mul_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000040006, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x0007, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0008, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x0009, 0x0000), :LoopVectorization, :vfmadd_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000700080009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000a, 0x0000), :LoopVectorization, :LOOPCONSTANTINSTRUCTION, LoopVectorization.OperationStruct(0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.constant, 0x000b, 0x0000), :LoopVectorization, :pow_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x000000000000000000000000000a000b, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000c, 0x0000), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x000000000000000000000000000c0009, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000d, 0x0000), :LoopVectorization, :mul_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x0000000000000000000000000001000d, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000e, 0x0000), :LoopVectorization, :sub_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000007, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, LoopVectorization.compute, 0x000f, 0x0000), :LoopVectorization, :vfmadd_fast, LoopVectorization.OperationStruct(0x00000000000000000000000000000001, ...

```

julia\> @report\_opt target\_modules=(@ **MODULE** ,) loopvectorization\_approach(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)  
No errors detected

Looks good to me?

kind regards

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 8, 2024, 5:00pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/16 "2024-02-08T17:00:29Z")

</div>

Use Cthulhu to descend into `_turbo_!`.  
The problem was in `promote_type` on the example I looked at before, so likely there again.

It is helpful to enable optimizations and warn types within descend.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 6:17pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/18 "2024-02-08T18:17:44Z")

</div>

Okay I found the issue! It is something with the inverse of gamma in the power calculation. If I write it directly no issue anymore on Julia 1.10:

```julia
using BenchmarkTools
using StaticArrays
using LoopVectorization

# Generate the needed constants
N = 6000
Cb = inv(Float64(100000))
xᵢⱼᶻ = rand(Float64,N*4)
xᵢⱼ = rand(SVector{3,Float64},N*4)
ρ₀ = Float64(1000)
g = Float64(9.81)
drhopLp = zeros(Float64, length(xᵢⱼ))
drhopLn = zeros(Float64, length(xᵢⱼ))

function loopvectorization_approach(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
    @tturbo for iter in eachindex(xᵢⱼᶻ)
        # IF I COMMENT '* -xᵢⱼᶻ[iter]' NO ALLOCATIONS, WHY?
        Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼᶻ[iter]
        ρᵢⱼᴴ = ρ₀ * ( ^( 1 + (Pᵢⱼᴴ * Cb), 0.14285714285714285) - 1)  
        Pⱼᵢᴴ = -Pᵢⱼᴴ
        ρⱼᵢᴴ = ρ₀ * ( ^( 1 + (Pⱼᵢᴴ * Cb), 0.14285714285714285) - 1)
        
        drhopLp[iter] = ρᵢⱼᴴ
        drhopLn[iter] = ρⱼᵢᴴ
    end
end

l = @benchmark loopvectorization_approach($Cb,$xᵢⱼᶻ,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn)

display(l)

julia> l = @benchmark loopvectorization_approach($Cb,$xᵢⱼᶻ,$ρ₀,$γ⁻¹,$g,$drhopLp,$drhopLn)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 339.900 μs … 4.994 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 358.600 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 382.900 μs ± 97.110 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

   ▅█▇▅▄▃▂▂▁▁▁▁▁ ▁▁▁▁ ▂
  █████████████████▇▇█▇█████████▇█▇▇▇▇▇▆▄▆▅▆▆▆▅▅▆▄▆▅▁▅▅▅▄▄▆▆▇▇ █
  340 μs Histogram: log(frequency) by time 718 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

```

Now on to see if I can trick it and use the variable.

EDIT: Both ways are Float64 values, I am confused 😅

Kind regards

---

<div class="post-metadata">

**Author:** ![mihalybaci](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mihalybaci/32/13528_2.png) [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)\
**Post date:** [February 8, 2024, 6:51pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/19 "2024-02-08T18:51:31Z")

</div>

If you’re really trying to squeeze every microsecond out of this, then it may be worth it to do

```julia
drhopLp = fill(ρ₀, length(xᵢⱼ))
drhopLn = fill(ρ₀, length(xᵢⱼ))

```

instead of creating them as arrays of zeroes. Then the last two lines in the loop would become

```julia
 drhopLp[iter] *= ρᵢⱼᴴ
 drhopLn[iter] *= ρⱼᵢᴴ

```

It is also possible to do `ρ₀*g` outside the loop as well, and in that same line, there are two multiplies by a negative, so those can probably be eliminated. That’s a couple of operations removed from the loop. Best to check that these things don’t change you’re results though.

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [February 8, 2024, 7:19pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/20 "2024-02-08T19:19:29Z")

</div>

Thanks for the hint, did some of it. Decided against fill, but will reconsider perhaps.

If anyone knows how to solve the issue with `^` in LoopVectorization, see: [Can this be written even faster? (CPU) - #18 by Ahmed\_Salih](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/18)

I would be really happy 🙂

Kind regards

---

<div class="post-metadata">

**Author:** ![Uve](https://avatars.discourse-cdn.com/v4/letter/u/82dd89/32.png) [@Uve](https://discourse.julialang.org/u/Uve)\
**Post date:** [February 8, 2024, 9:07pm UTC](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924/21 "2024-02-08T21:07:46Z")

</div>

I don’t know why but this seems to work

```julia
julia> faux(ρ₀, P, Cb, γ⁻¹) = ρ₀ * ( ^( 1 + (P * Cb), γ⁻¹) - 1)
faux (generic function with 1 method)

julia> function loopvectorization_approach3(Cb,xᵢⱼᶻ,ρ₀,γ⁻¹,g,drhopLp,drhopLn)
           @tturbo for iter in eachindex(xᵢⱼᶻ)
               # IF I COMMENT '* -xᵢⱼᶻ[iter]' NO ALLOCATIONS, WHY?
               Pᵢⱼᴴ = ρ₀ * (-g) * -xᵢⱼᶻ[iter]
               Pⱼᵢᴴ = -Pᵢⱼᴴ
               ρᵢⱼᴴ = faux(ρ₀, Pᵢⱼᴴ, Cb, γ⁻¹)
               ρⱼᵢᴴ = faux(ρ₀, Pⱼᵢᴴ, Cb, γ⁻¹)
               drhopLp[iter] = ρᵢⱼᴴ
               drhopLn[iter] = ρⱼᵢᴴ
           end
       end
loopvectorization_approach3 (generic function with 2 methods)

julia> l3 = @benchmark loopvectorization_approach3($Cb,$xᵢⱼᶻ,$ρ₀, $γ⁻¹, $g,$drhopLp,$drhopLn)
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max): 35.100 μs … 557.100 μs ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 73.200 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 73.923 μs ± 12.997 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▃▂ █▄▂▆▆▃▂▂▁ ▂
  ███▆▅▅▁▁▁▃▁▁▃▃▁▁▁▄▃▁▁▃▁▁▁▁▁▁▁▃▃▁▁▁▁▁▁▄▆██████████▇▆▇▅▆▅▆▅▁▆▆ █
  35.1 μs Histogram: log(frequency) by time 92.5 μs <
 Memory estimate: 0 bytes, allocs estimate: 0.
```

[Next page](https://discourse.julialang.org/t/can-this-be-written-even-faster-cpu/109924.md?page=2)