# How can this code be faster than MATLAB?

**URL:** <https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196>\
**Category:** Performance\
**Tags:** question\
**Created:** [July 5, 2023, 3:44am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196 "2023-07-05T03:44:10Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![XuJingye2022](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xujingye2022/32/38972_2.png) [@XuJingye2022](https://discourse.julialang.org/u/XuJingye2022)\
**Post date:** [July 5, 2023, 3:44am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/1 "2023-07-05T03:44:11Z")

</div>

```julia
using BenchmarkTools

zvec, δvec = randn(10000) * 1e-2, randn(10000) * 1e-3

# Physical constant
const restenergy = 0.511e6
const clight = 3e8

# Computation settings
nturns = 10001

# Machine parameter
struct Par
    v1::Float64
    v2::Float64
    ϕs::Float64
    ϕ2s::Float64
    h1::Int64
    h::Int64
    circum::Float64
    centerenergy::Float64
    αc::Float64
    r::Float64
end
function Par(v1, v2, ϕs, ϕ2s, h1, h, circum, centerenergy, αc)
    Par(v1, v2, ϕs, ϕ2s, h1, h, circum, centerenergy, αc, v2 / v1)
end

function _revolution_cache(p::Par)
    c1 = p.v1 / p.centerenergy
    c2 = sin(p.ϕs) + p.r * sin(p.ϕ2s)
    c3 = 1 / (p.centerenergy / restenergy)^2
    k1 = 2 * p.h1 * pi / p.circum
    k2 = 2 * p.h1 * p.h * pi / p.circum
    (; c1, c2, c3, k1, k2)
end

function revolution_with_cache(z::Number, δ::Number, p::Par, c)
    δ = δ + c.c1 * (sin(p.ϕs - c.k1 * z) + p.r * sin(p.ϕ2s - c.k2 * z) - c.c2)
    η = p.αc - c.c3 / (1 + δ)^2
    z = z - p.circum * η * δ
    z, δ
end

function revolution!(zvec::Vector, δvec::Vector, nturns::Int64, p::Par)
    c = _revolution_cache(p)
    nparticles = length(zvec)
    @inbounds for j = 1:nturns
        @inbounds @simd for i = 1:nparticles
            zvec[i], δvec[i] = revolution_with_cache(zvec[i], δvec[i], p, c)
        end
    end
    zvec, δvec
end

p = Par(3.6395053705870048e+06, 655751.6971214055228756, 2.1575815005097385, 6.0620746573056303, 756, 3, 1360.4, 6e9, 1.56e-5)

@btime revolution!($zvec, $δvec, $nturns, $p);

```

Hello everyone. I’m trying to modify this code faster than MATLAB.

```matlab
zvec = randn(10000,1) * 1e-2;
deltavec = randn(10000,1) * 1e-3;

%%

restenergy = 0.511e6;
clight = 3e8;

nturns = 10001;
nparticle = length(zvec);

v1 = 3.6395053705870048e+06;
v2 = 655751.6971214055228756;
r = v2/v1;
phi_s = 2.1575815005097385;
phi_2s = 6.0620746573056303;
h1 = 756;
h = 3;
circum = 1360.4;
centerenergy = 6e9;
alpha_c = 0.0000156;

%%

tic;

for i = 1:nturns
    deltavec = deltavec + v1/centerenergy * (sin(phi_s - 2 * h1 * pi / circum .* zvec) - sin(phi_s) + r * sin(phi_2s - 2 * h1 * h * pi / circum .* zvec) - r * sin(phi_2s));
    etavec = alpha_c - 1./(centerenergy/restenergy .* (1 + deltavec)).^2;
    zvec = zvec - circum .* etavec .* deltavec;
end

toc

```

On my computer, the Julia code runs for approximately 1.73 seconds, while the MATLAB code runs for 0.63 seconds.

How else can I optimize the code to double or triple its speed?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 3:57am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/2 "2023-07-05T03:57:01Z")

</div>

Please read: [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/)

Great that you are already using BenchmarkTools!

Did you check if your inner function `revolution_with_cache` is allocating any memory, e.g. using the @time macro?

---

<div class="post-metadata">

**Author:** ![XuJingye2022](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xujingye2022/32/38972_2.png) [@XuJingye2022](https://discourse.julialang.org/u/XuJingye2022)\
**Post date:** [July 5, 2023, 4:21am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/3 "2023-07-05T04:21:48Z")

</div>

```julia
julia> c = _revolution_cache(p);

julia> @time revolution_with_cache(0.0, 0.0, p, c);
  0.004438 seconds (133 allocations: 7.297 KiB, 99.78% compilation time)

julia> @time revolution_with_cache(0.0, 0.0, p, c);
  0.000005 seconds (1 allocation: 32 bytes)

```

Hello,

`0.000005`s seems okay. Calling 10000 times takes `0.05`s.

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 4:27am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/4 "2023-07-05T04:27:39Z")

</div>

> [@XuJingye2022](#):
>
> `0.000005`s seems okay. Calling 10000 times takes `0.05`s.

But you are calling it 10000\*10000 times, or am I wrong?

You have a nested loop…

---

<div class="post-metadata">

**Author:** ![XuJingye2022](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xujingye2022/32/38972_2.png) [@XuJingye2022](https://discourse.julialang.org/u/XuJingye2022)\
**Post date:** [July 5, 2023, 4:29am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/5 "2023-07-05T04:29:05Z")

</div>

😵Sorry, my fault

So it takes …`500` seconds? I don’t know, I’m a novice 😵

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 4:39am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/6 "2023-07-05T04:39:04Z")

</div>

You could pass `zvec` and `δvec` and i to this function to avoid the allocation and modify the elements in place…

You could also try to use `@inline` in front of the function to avoid the call overhead…

---

<div class="post-metadata">

**Author:** ![XuJingye2022](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xujingye2022/32/38972_2.png) [@XuJingye2022](https://discourse.julialang.org/u/XuJingye2022)\
**Post date:** [July 5, 2023, 5:22am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/7 "2023-07-05T05:22:26Z")

</div>

```julia
function revolution!(zvec::Vector, δvec::Vector, nturns::Int64, p::Par)
    c = _revolution_cache(p)
    nparticles = length(zvec)
    @inbounds for j = 1:nturns
        @inbounds @simd for i = 1:nparticles
            δvec[i] += c.c1 * (sin(p.ϕs - c.k1 * zvec[i]) + p.r * sin(p.ϕ2s - c.k2 * zvec[i]) - c.c2)
            zvec[i] -= p.circum * (p.αc - c.c3 / (1 + δvec[i])^2) * δvec[i]
        end
    end
    zvec, δvec
end

```

I tried this.  
This function is absolutely inlined and replace values in place.  
But,

```julia
julia> @btime revolution!($zvec, $δvec, $nturns, $p);
  1.735 s (0 allocations: 0 bytes)

julia> @btime revolution!($zvec, $δvec, $nturns, $p);
  1.729 s (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 5:55am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/8 "2023-07-05T05:55:19Z")

</div>

I can reproduce you results. The next question is, is using Matlab multiple threads to calculate the result?

This is the code:

```julia
for i = 1:nturns
    deltavec = deltavec + v1/centerenergy * (sin(phi_s - 2 * h1 * pi / circum .* zvec) - sin(phi_s) + r * sin(phi_2s - 2 * h1 * h * pi / circum .* zvec) - r * sin(phi_2s));
    etavec = alpha_c - 1./(centerenergy/restenergy .* (1 + deltavec)).^2;
    zvec = zvec - circum .* etavec .* deltavec;
end

```

Any Matlab expert around who could answer this question?

---

<div class="post-metadata">

**Author:** ![ufechner7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ufechner7/32/51363_2.png) [@ufechner7](https://discourse.julialang.org/u/ufechner7)\
**Post date:** [July 5, 2023, 6:07am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/9 "2023-07-05T06:07:13Z")

</div>

Ok, if I start matlab with one thread:

```julia
matlab --use-single-comp-thread

```

your code takes 1.284115s to execute, the Julia version needs 955.048 ms, so Julia is faster.

The question now is, how to modify your Julia code such that it uses multi-threading…

I think I leave the answer to others, must work now…

But have a look at: [GitHub - JuliaSIMD/LoopVectorization.jl: Macro(s) for vectorizing loops.](https://github.com/JuliaSIMD/LoopVectorization.jl/tree/main)

UPDATE:  
When using LoopVectorization I get:

```julia
# Matlab
# singlethreaded: 1.28s
# multithreaded: 0.39s
# Julia with @tturbo macro:
# singlethreaded: 0.26s
# multithreaded: 0.035s

```

So on Ryzen 7950x Julia can be 11 times faster than Matlab.

---

<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:** [July 5, 2023, 6:08am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/10 "2023-07-05T06:08:45Z")

</div>

LoopVectorization.jl will help those `sin` calls.

```julia
julia> @btime revolution2!($zvec, $δvec, $nturns, $p);
  1.994 s (0 allocations: 0 bytes)

julia> @btime revolution3!($zvec, $δvec, $nturns, $p);
  251.768 ms (0 allocations: 0 bytes)

julia> @btime revolution4!($zvec, $δvec, $nturns, $p);
  30.650 ms (0 allocations: 0 bytes)

```

This was with

```julia
julia> function revolution3!(zvec::Vector, δvec::Vector, nturns::Int64, p::Par)
           c = _revolution_cache(p)
           nparticles = length(zvec)
           @inbounds for j = 1:nturns
               @turbo for i = 1:nparticles
                   δvec[i] += c.c1 * (sin(p.ϕs - c.k1 * zvec[i]) + p.r * sin(p.ϕ2s - c.k2 * zvec[i]) - c.c2)
                   zvec[i] -= p.circum * (p.αc - c.c3 / (1 + δvec[i])^2) * δvec[i]
               end
           end
           zvec, δvec
       end
revolution3! (generic function with 1 method)

```

While `revolution2!` was `@inbounds @simd` instead, and `revolution4!` with `@tturbo`.  
Results will vary, but it should help.

---

<div class="post-metadata">

**Author:** ![longemen3000](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/longemen3000/32/7298_2.png) [@longemen3000](https://discourse.julialang.org/u/longemen3000)\
**Post date:** [July 5, 2023, 6:12am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/11 "2023-07-05T06:12:36Z")

</div>

can this help?

```julia
@inbounds for j = 1:nturns
    Threads.@threads for i = 1:nparticles
        δvec[i] += c.c1 * (sin(p.ϕs - c.k1 * zvec[i]) + p.r * sin(p.ϕ2s - c.k2 * zvec[i]) - c.c2)
        zvec[i] -= p.circum * (p.αc - c.c3 / (1 + δvec[i])^2) * δvec[i]
    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:** [July 5, 2023, 6:14am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/12 "2023-07-05T06:14:38Z")

</div>

I get

```julia
julia> function revolution5!(zvec::Vector, δvec::Vector, nturns::Int64, p::Par)
           c = _revolution_cache(p)
           nparticles = length(zvec)
           @inbounds for j = 1:nturns
               Threads.@threads for i = 1:nparticles
                   @inbounds begin; δvec[i] += c.c1 * (sin(p.ϕs - c.k1 * zvec[i]) + p.r * sin(p.ϕ2s - c.k2 * zvec[i]) - c.c2)
                   zvec[i] -= p.circum * (p.αc - c.c3 / (1 + δvec[i])^2) * δvec[i]; end
               end
           end
           zvec, δvec
       end
revolution5! (generic function with 1 method)

julia> @btime revolution5!($zvec, $δvec, $nturns, $p);
  437.087 ms (1666342 allocations: 185.60 MiB)

julia> Threads.nthreads(), Sys.CPU_THREADS
(28, 28)

```

Which is an improvement, but on 28 threads, it is still substantially slower than `@turbo` on a single thread.

---

<div class="post-metadata">

**Author:** ![XuJingye2022](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xujingye2022/32/38972_2.png) [@XuJingye2022](https://discourse.julialang.org/u/XuJingye2022)\
**Post date:** [July 5, 2023, 6:22am UTC](https://discourse.julialang.org/t/how-can-this-code-be-faster-than-matlab/101196/13 "2023-07-05T06:22:04Z")

</div>

Thank you all.