# Julia can be dramatically slower than Matlab when solving 2D PDE

**URL:** <https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321>\
**Category:** New to Julia\
**Tags:** performance, matlab, pde\
**Created:** [November 6, 2024, 8:46am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321 "2024-11-06T08:46:27Z")\
**Posts on this page:** 20\
**Page:** 2

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 6, 2024, 1:51pm UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/22 "2024-11-06T13:51:24Z")

</div>

Dear Oscar,  
Thanks for your suggestion. For now, I only work with partial differential equations

---

<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 6, 2024, 2:02pm UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/23 "2024-11-06T14:02:29Z")

</div>

You can use an arbitrary ODE solver to timestep your PDEs. I believe something like this should work.

```julia
using OrdinaryDiffEqTsit5
# formulate energy as ODE
@views function energy3(dT, T, p, t)
	Fo, h, N = p
	@. dT[2:N-1,2:N-1] = (Fo / h) * ((T[3:N,2:N-1]) - 2 * T[2:N-1,2:N-1] + T[1:N-2,2:N-1]) + (T[2:N-1,3:N] - 2 * T[2:N-1,2:N-1] + T[2:N-1,1:N-2])
    dT[1:N,1] .= 0;
    dT[1:N,N] .= 0;
    dT[1,1:N] .= dT[2,1:N];
    dT[N,1:N] .= dT[N-1,1:N];
end

# parameters
N = 101;
Fo = 0.01;
Th = 0.5;
Tc = -0.5;
h = 1/(N-1);

# initial T
T0 = zeros(N, N)
T0[1:N,1] .= Th;
T0[1:N,N] .= Tc;
prob = ODEProblem(energy3, T0, (0, 1000), (Fo, h, N))
sol = solve(prob, Tsit5())

```

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [November 6, 2024, 2:14pm UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/24 "2024-11-06T14:14:42Z")

</div>

> [@lmiq](#):
>
> while in Matlab (it seems) it _does_ copy the array

Semantically it does copy the array. Implementation-wise it’s done with copy-on-write, so they share the data until either reference gets mutated.  
The same happens with arguments to functions.

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 8, 2024, 3:35am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/25 "2024-11-08T03:35:03Z")

</div>

I dare to ask another question. Now let’s switch to GPU computing. Firstly, let’s execute the code on CPU

```julia
N = 2001;
Fo = 0.01;
h = 1/(N-1);
T=(zeros(N,N));
#T=CuArray((zeros(N,N)));
Th = 0.5;
Tc = -0.5;
Time = 100;
dt=0.01;
function energy(Time, N, T, Th, Fo, dt, Tc, h)
    for iter = 1:Time
        Tn = T;
@. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
T[1:N,1] .= Th;
T[1:N,N] .= Tc;
T[1,1:N] .= T[2,1:N];
T[N,1:N] .= T[N-1,1:N];
    end
end
@time energy(Time, N, T, Th, Fo, dt, Tc, h)
8.872417 seconds (1.09 M allocations: 20.899 GiB, 4.28% gc time, 3.45% compilation time)

```

Now, let’s swtich to GPU

```julia
N = 2001;
Fo = 0.01;
h = 1/(N-1);
#T=(zeros(N,N));
T=CuArray((zeros(N,N)));
Th = 0.5;
Tc = -0.5;
Time = 100;
dt=0.01;
function energy(Time, N, T, Th, Fo, dt, Tc, h)
    for iter = 1:Time
        Tn = T;
@. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
T[1:N,1] .= Th;
T[1:N,N] .= Tc;
T[1,1:N] .= T[2,1:N];
T[N,1:N] .= T[N-1,1:N];
    end
end
@time energy(Time, N, T, Th, Fo, dt, Tc, h)
0.685531 seconds (1.17 M allocations: 64.837 MiB, 3.00% gc time, 60.79% compilation time)

```

It seems like a good speed up. However, Matlab executed this code in 0.21 seconds on GPU. Can someone suggest any easy-to-understand optimization tips? Or should i straightaway try to write a CUDA kernel?

---

<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 8, 2024, 4:31am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/26 "2024-11-08T04:31:36Z")

</div>

On a GPU, you definitely want to use Float32 for all your numbers. You also definitely should be using an ODE solver as previously mentioned.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [November 8, 2024, 4:38am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/27 "2024-11-08T04:38:15Z")

</div>

Why does Julia not use such a mechanism?

---

<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 8, 2024, 4:58am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/28 "2024-11-08T04:58:27Z")

</div>

Copy on write makes multithreading, and getindex more expensive. It’s very far from a free optimization.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 8, 2024, 6:26am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/29 "2024-11-08T06:26:45Z")

</div>

> [@photor](#):
>
> Why does Julia not use such a mechanism?

Also, if `Arr2 = Arr1`, we _want_ mutations to `Arr2` to be visible in `Arr1`. This sharing is central to Julia semantics, and is sorely missed in Matlab.

Or do you mean that slices could be views until they are mutated?

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [November 8, 2024, 6:55am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/30 "2024-11-08T06:55:05Z")

</div>

Yes, absolutely.

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 8, 2024, 7:38am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/31 "2024-11-08T07:38:19Z")

</div>

In future, of course it will be the float32. I just did the test case. In matlab, the double precision was also used. What concerns the ODE solver, i don’t feel that it is suitable in my case since i study different numerical schemes for computaional fluid dynamics in my research

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 8, 2024, 11:54am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/32 "2024-11-08T11:54:50Z")

</div>

There are two important optimization:

1. Use `@views` as mentioned before. That will make slices (like `Tn[2:N-1,2:N-1]`) not allocate new temporary arrays.

2. Again: the line `Tn = T` does not really makes sense here. And this turns out to be important for performance. With `Tn = T` you are just saying that `Tn` and `T` are **the same** array. This is **wrong** when you use the loop version of the code. In this last version, with the broadcast, it is not wrong, but because, to avoid aliasing (to avoid concurrent writing to the array `T`), the compiler will create temporary copies of the slices of the array where it cannot prove that concurrency does not occur.

In summary, use the function like this:

```julia
@views function energy(Time, N, T, Th, Fo, dt, Tc, h)
      Tn = similar(T)
      for iter = 1:Time
          Tn .= T;
# the rest is identical

```

I get, before these changes:

```julia-repl
julia> @btime energy($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h)
  109.071 ms (66900 allocations: 1.82 MiB)

```

After these two changes:

```julia-repl
julia> @btime energy($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h)
  2.984 ms (43405 allocations: 978.33 KiB)

```

Thus a ~30x speedup.

Edit: The path for obtaining optimal code is to eliminate the allocation of intermediates. Here, for example, I add `Tn` as an additional variable, to not allocate inside the function. Then, I get this:

```julia-repl
julia> @views function energy(Time, N, T, Th, Fo, dt, Tc, h, Tn) # Additional parameter Tn
           for iter = 1:Time
               Tn .= T;
       @. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
       T[1:N,1] .= Th;
       T[1:N,N] .= Tc;
       T[1,1:N] .= T[2,1:N];
       T[N,1:N] .= T[N-1,1:N];
           end
       end
energy (generic function with 2 methods)

julia> T = zeros(N,N); Tn = similar(T); # running on CPU

julia> @btime energy($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h, $Tn)
  1.382 s (0 allocations: 0 bytes)

```

Very well, the function does not allocate any intermediate when running on the CPU. This is probably a reasonably good implementation, and running it on the GPU confirms it. The previous version, even with the use of `@views`, was allocating and, thus was slower:

```julia-repl
julia> @views function energy_old(Time, N, T, Th, Fo, dt, Tc, h) # Additional parameter Tn
           for iter = 1:Time
               Tn = T;
       @. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
       T[1:N,1] .= Th;
       T[1:N,N] .= Tc;
       T[1,1:N] .= T[2,1:N];
       T[N,1:N] .= T[N-1,1:N];
           end
       end
energy_old (generic function with 1 method)

julia> @btime energy_old($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h)
  3.193 s (800 allocations: 11.91 GiB)

```

The allocations coming from the fact that the compiler was unable to guarantee, in the broadcast, that concurrent writing to T was not occurring.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 8, 2024, 12:18pm UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/33 "2024-11-08T12:18:46Z")

</div>

> [@Alex90](#):
>
> What concerns the ODE solver, i don’t feel that it is suitable in my case since i study different numerical schemes for computaional fluid dynamics in my research

You can use an ODE solver for that of course: there’s about 300 already available and tested for correctness.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [November 8, 2024, 11:57pm UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/34 "2024-11-08T23:57:54Z")

</div>

Just to expand on the recommendations about functions, it is idiomatic in Julia to indicate a mutating function by adding ! to its name, and to have the mutated variable as first argument. Even better, follow the function template of `OrdinaryDifferentialEquations.jl` like @Oscar_Smith did, but with !:  
`energy3!(dT, T, p, t)`

You may be interested in PDE packages, many of which are organized under [`JuliaFEM.jl`](https://www.juliafem.org/JuliaFEM.jl/latest/), for example `HeatTransfer.jl`.

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [November 9, 2024, 12:03am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/35 "2024-11-09T00:03:26Z")

</div>

> You may be interested in PDE packages, many of which are organized under [`JuliaFEM.jl`](https://www.juliafem.org/JuliaFEM.jl/latest/), for example `HeatTransfer.jl`.

JuliaFEM.jl is extremely outdated now. HeatTransfer.jl hasn’t even been updated in 4 years. A slightly more up to date reference of the available PDE packages is probably(?) [GitHub - JuliaPDE/SurveyofPDEPackages: Survey of the packages of the Julia ecosystem for solving partial differential equations](https://github.com/JuliaPDE/SurveyofPDEPackages)

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 9, 2024, 9:37am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/36 "2024-11-09T09:37:05Z")

</div>

Many thanks for clear and concise explanation! With these changes, Julia executes the GPU code ~1.6 faster than Matlab. However, when i use Float32 variables, i get almost the same performance as Float64.

```julia
N = 2001;
Nx = N; Ny = N;
Fo = 0.01;
h = 1/(N-1);
#T=(zeros(N,N));
T=CuArray((zeros(N,N)));
Tn = similar(T)
Th = 0.5;
Tc = -0.5;
Time = 10000;
dt=0.01; 
@views function energy_2(Time, N, T, Th, Fo, dt, Tc, h, Tn)
    for iter = 1:Time
        Tn .= T;
@. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
T[1:N,1] .= Th;
T[1:N,N] .= Tc;
T[1,1:N] .= T[2,1:N];
T[N,1:N] .= T[N-1,1:N];
    end
end
##@time energy_2(Time, N, T, Th, Fo, dt, Tc, h, Tn)
@btime energy_2($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h, $Tn)
17.028 s (3079034 allocations: 517.09 MiB)

```

And Float32

```julia
N = 2001;
Fo = 0.01;
h = 1/(N-1);
T=CUDA.zeros(Float32,N,N);
Tn = similar(T)
Th = 0.5;
Tc = -0.5;
Time = 10000;
dt=0.01;
Fo_gpu = CUDA.fill(Float32(Fo),1)
dt_gpu = CUDA.fill(Float32(dt),1)
Th_gpu = CUDA.fill(Float32(Th),1)
Tc_gpu = CUDA.fill(Float32(Tc),1)
h_gpu = CUDA.fill(Float32(h),1);
@views function energy(Time, N, T, Th, Fo, dt, Tc, h, Tn)
    for iter = 1:Time
        Tn .= T;
@. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
T[1:N,1] .= Th;
T[1:N,N] .= Tc;
T[1,1:N] .= T[2,1:N];
T[N,1:N] .= T[N-1,1:N];
    end
end
##@time energy(Time, N, T, Th, Fo, dt, Tc, h, Tn)
@btime energy($Time, $N, $T, $Th, $Fo, $dt, $Tc, $h, $Tn)
15.284 s (3628806 allocations: 577.51 MiB)

```

Is it an expectable result in this test case?

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 9, 2024, 9:41am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/37 "2024-11-09T09:41:43Z")

</div>

Dear Chris,  
Thank you for this note. I develop hybrid lattice Boltzmann schemes which i believe still unavailable in commercial softwares and free libraries or packages

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [November 9, 2024, 10:24am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/38 "2024-11-09T10:24:56Z")

</div>

Unexpectedly low performance might come from `Float64` in the intermediate computation. Maybe a better pattern would be like this:

```julia
function energy(
    Time::Int32,
    T::AbstractMatrix{F},
    Th::AbstractMatrix{F},
    Fo::F,
    dt::F,
    Tc::F,
    h::F,
    Tn::AbstractMatrix{F}
) where {F<:AbstractFloat}
    N = size(T, 1)
    @views for iter = 1:Time
        Tn .= T;
        @. T[2:N-1,2:N-1] = Tn[2:N-1,2:N-1] + Fo * dt * ((Tn[3:N,2:N-1] - 2 * Tn[2:N-1,2:N-1] + Tn[1:N-2,2:N-1]) / h + (Tn[2:N-1,3:N] - 2 * Tn[2:N-1,2:N-1] + Tn[2:N-1,1:N-2]) / h)
        T[1:N,1] .= Th;
        T[1:N,N] .= Tc;
        T[1,1:N] .= T[2,1:N];
        T[N,1:N] .= T[N-1,1:N];
    end
end

function energy(
    Time::Integer,
    T::AbstractMatrix{F},
    Th::AbstractMatrix{F},
    Fo::Real,
    dt::Real,
    Tc::Real,
    h::Real,
    Tn::AbstractMatrix{F}
) where {F<:AbstractFloat}
    Fo1, dt1, Tc1, h1 = convert.(F, (Fo, dt, Tc, h))
    return energy(Int32(Time), T, Th, Fo1, dt1, Tc1, h1, Tn)
end

```

This way, the scalar arguments will be first converted into `eltype(T)`, to avoid mixed-precision operations.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [November 9, 2024, 10:27am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/39 "2024-11-09T10:27:17Z")

</div>

> [@Alex90](#):
>
> `Fo_gpu = CUDA.fill(Float32(Fo),1)`

I think you need only `Fo = 1.0f0` (a scalar of 32bits).

And you forgot to change the function call to use the new values.

That said, I tested with 32bits in my machine and didn’t get a significant improvement either. Maybe some more experienced CUDA user can explain.

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 9, 2024, 10:38am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/40 "2024-11-09T10:38:18Z")

</div>

I noticed that too. But after changing the values inside the function, i got the same performance result. Anyway, many thanks for this useful discussion

---

<div class="post-metadata">

**Author:** ![Alex90](https://avatars.discourse-cdn.com/v4/letter/a/6de8d8/32.png) [@Alex90](https://discourse.julialang.org/u/Alex90)\
**Post date:** [November 9, 2024, 10:48am UTC](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321/41 "2024-11-09T10:48:38Z")

</div>

It is very unclear for me what intermediate computations can be in Float64 since all variables are Float32 🤔 After computation, i just in case checked every variable and they all were Float32.  
Many thanks for your code modification. Unfortunately, i did not understand how to run it properly.

[Previous page](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321.md?page=1)

[Next page](https://discourse.julialang.org/t/julia-can-be-dramatically-slower-than-matlab-when-solving-2d-pde/122321.md?page=3)