# Code optimization using SIMD or LoopVectorization

**URL:** <https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782>\
**Category:** Performance\
**Tags:** question\
**Created:** [November 4, 2023, 2:45am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782 "2023-11-04T02:45:32Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [November 4, 2023, 2:45am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/1 "2023-11-04T02:45:32Z")

</div>

Dear Julia experts,

I have implemented the probability density function for the Pearson-IV distribution. This implementation is based on the R-package “PearsonDS”. Here is the code:

```julia
import SpecialFunctions: lgamma, gamma
using SIMD

const M_LN_SQRT_PI = 0.5log(pi)

@inline C_logPearsonIVnorm0(m, nu) = 
    (-M_LN_SQRT_PI - lgamma(m) - lgamma(m-0.5) + 2real(lgamma(m+0.5nu*im)))

function dpearsonIV_simd(
    x::Vector{Float64}; 
    params=[],
    log_p=false
    )::Vector{Float64}
    n = length(x)
    (m, nu, location, scale) = params
    @assert ((scale > 0) && (m > 0.5))
    inv_scale = 1/scale
    ret = similar(x)
    t = Vec{4,Float64}((0,0,0,0))
    s = Vec{4,Float64}((0,0,0,0))
    k = C_logPearsonIVnorm0(m, nu) - log(scale)
    @inbounds for i = 1:4:n
            t = vload(Vec{4,Float64}, x, i)
            @fastmath t = inv_scale * (t-location)
            @fastmath s = Vec{4,Float64}((atan(t[1]),atan(t[2]),atan(t[3]),atan(t[4])))
            # @fastmath s = (-nu) * s
            # @fastmath t = t*t+1
            # @fastmath t = log(t)
            @fastmath t = (-nu) * s + (-m) * log(t*t+1) + k
            if !log_p
                @fastmath t = exp(t)
            end
            vstore(t, ret, i)
    end
    return ret
end

```

I have benchmarked this code and found that the function costs 31 microseconds on a 1000-element input array x. But the essential computations, namely the multipilcations and `atan`, cost approximately 8 microseconds.

I want to learn some suggestions on my code to make further improvements, since it will be used frequently in very complicated fitting algorithms. Thanks!

---

<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:** [November 4, 2023, 3:08am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/2 "2023-11-04T03:08:06Z")

</div>

The main suggestion I have is that you are likely this at too low a level. This looks like it would be a good place for `LoopVectorization`. With that, this would simplify to

```julia
using LoopVectorization
function dpearsonIV_simd(
    x::Vector{Float64}; 
    params,
    log_p=false)
    n = length(x)
    (m, nu, location, scale) = params
    @assert ((scale > 0) && (m > 0.5))
    inv_scale = 1/scale
    k = C_logPearsonIVnorm0(m, nu) - log(scale)
    ret = @turbo @. atan(inv_scale * (x-location))
    @turbo @. ret = (-nu) * s + (-m) * log(ret*ret+1) + k
    if !log_p
        @turbo @. ret = exp(ret)
    end
    return ret
end

```

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [November 4, 2023, 3:33am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/3 "2023-11-04T03:33:30Z")

</div>

> [@Oscar\_Smith](#):
>
> ```julia
> using LoopVectorization
> function dpearsonIV_simd(
> x::Vector{Float64}; 
> params,
> log_p=false)
> n = length(x)
> (m, nu, location, scale) = params
> @assert ((scale > 0) && (m > 0.5))
> inv_scale = 1/scale
> k = C_logPearsonIVnorm0(m, nu) - log(scale)
> ret = @turbo @. atan(inv_scale * (x-location))
> @turbo @. ret = (-nu) * s + (-m) * log(ret*ret+1) + k
> if !log_p
> @turbo @. ret = exp(ret)
> end
> return ret
> end
> 
> ```

Oh this is great! `dpearsonIV_lv` is your version with `LoopVecterization.jl`. You can see that it is indeed helpful by cutting the execution time more than half. I have monitored the CPU usage and confirmed that the main job is done within one core.

Thank you @Oscar_Smith !

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/0/5044de9a0d005c9214fddeff9d36969a65b5c2e1.png)

---

<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:** [November 4, 2023, 3:34am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/4 "2023-11-04T03:34:36Z")

</div>

Oh yeah, this is probably slow enough that you might want to multithread it by using `@tturbo` instead of `@turbo`

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [November 4, 2023, 4:06am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/6 "2023-11-04T04:06:18Z")

</div>

That would be even nicer. I have checked the `LoopVectorization,jl` source code, but I don’t have a sharp eye to pin-point the place where `@turbo` or rather `turbo_macro()` does the real job on my case to improve performance. But I wish to learn that. Can you give me some hints?

---

<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:** [November 4, 2023, 4:21am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/7 "2023-11-04T04:21:20Z")

</div>

It’s just

```julia
function dpearsonIV_simd(
    x::Vector{Float64}; 
    params,
    log_p=false)
    n = length(x)
    (m, nu, location, scale) = params
    @assert ((scale > 0) && (m > 0.5))
    inv_scale = 1/scale
    k = C_logPearsonIVnorm0(m, nu) - log(scale)
    ret = @tturbo @. atan(inv_scale * (x-location))
    @tturbo @. ret = (-nu) * s + (-m) * log(ret*ret+1) + k
    if !log_p
        @tturbo @. ret = exp(ret)
    end
    return ret
end

```

---

<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:** [November 4, 2023, 5:21am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/8 "2023-11-04T05:21:20Z")

</div>

> [@Lian\_Yunlong](#):
>
> pin-point the place where `@turbo` or rather `turbo_macro()` does the real job on my case to improve performance. But I wish to learn that. Can you give me some hints?

There’s not really a specific “the place”.  
But, compared to SIMD.jl, it uses SLEEFPirates.jl for faster `log` and `astandard.

Using VectorizationBase.jl instead of SIMD.jl might give you similar performance, but LoopVectorization.jl will try and do more clever things (that may or may not pay off).

@Oscar, might be worth fusing loops if using `@tturbo`? I haven’t benchmarks or checked assembly, but the tradeoff is probably register spills vs threading overhead.  
Register spill cost will obviously be different for 16 vs 32 architectural/ named registers.

---

<div class="post-metadata">

**Author:** ![Lian\_Yunlong](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lian_yunlong/32/14437_2.png) [@Lian\_Yunlong](https://discourse.julialang.org/u/Lian_Yunlong)\
**Post date:** [November 4, 2023, 7:13am UTC](https://discourse.julialang.org/t/code-optimization-using-simd-or-loopvectorization/105782/9 "2023-11-04T07:13:21Z")

</div>

> [@Elrod](#):
>
> But, compared to SIMD.jl, it uses SLEEFPirates.jl for faster `log` and `astandard.

A simple replacement of `atan` by `SLEEFPirates.atan_fast` give me an improvement of ~1.5x:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/2/4/24aaf06247031b58c67bcb9700a8f66a38cc660c.png)  
(To get the above benchmark result, I have just replaced `atan` by `atan_fast` in my original implementation of `dpearsonIV_simd()` and all other settings remain the same.)