# Simple PDE benchmark

**URL:** <https://discourse.julialang.org/t/simple-pde-benchmark/7430>\
**Category:** Numerics\
**Created:** [December 1, 2017, 8:16am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430 "2017-12-01T08:16:49Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 1, 2017, 8:16am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/1 "2017-12-01T08:16:49Z")

</div>

Hi,

I made a simple benchmark to try and convince colleagues to switch to Julia. It might be of interest to people, so here it is : [https://github.com/antoine-levitt/benchmark\_heat](https://github.com/antoine-levitt/benchmark_heat). It’s the most simple PDE code I could cook up (1D, explicit euler, finite differences), and it’s designed to highlight differences between languages for simple kernels, and in particular loops vs vectorized code. I compare C (gcc and clang with various flags), julia (loops and broadcast, different annotations/versions), and numpy/scilab/octave/matlab (loops and vector)

Some highlights of my very rudimentary benchmarking (best of 3 on my laptop):

- There’s still a penalty for writing vector code in julia (doesn’t do SIMD?)
- Julia vectorized is the fastest of any vectorized versions, and in particular has a factor of three speedup between 0.6 and 0.7!
- Julia loops beats everybody, even the C version when compiled with clang (which I don’t understand)
- Loops in matlab have gotten fast in the latest versions, to the point where they even beat the vectorized version. That must have been huge work, kudos to them!

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [December 1, 2017, 8:49am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/2 "2017-12-01T08:49:30Z")

</div>

Thanks!

> [@antoine-levitt](#):
>
> Loops in matlab have gotten fast in the latest versions, to the point where they even beat the vectorized version.

Last time I checked (a few years ago), the restrictions on the loops to allow JITing were pretty severe. E.g. no call to m-file defined functions.

---

<div class="post-metadata">

**Author:** ![lobingera](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lobingera/32/211_2.png) [@lobingera](https://discourse.julialang.org/u/lobingera)\
**Post date:** [December 1, 2017, 9:03am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/3 "2017-12-01T09:03:00Z")

</div>

“fast in the latest versions” - “Last time I checked (a few years ago)” … maybe you should recheck

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [December 1, 2017, 10:01am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/4 "2017-12-01T10:01:01Z")

</div>

It looks like your vectorised version is allocating slices. What if you use views instead?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 1, 2017, 10:12am UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/5 "2017-12-01T10:12:25Z")

</div>

I do:

```julia
@views @. u_next[2:Nx-1] = u[2:Nx-1] + D*dt/(dx*dx)*(u[3:Nx]-2u[2:Nx-1]+u[1:Nx-2]) - v*dt/dx*(u[3:Nx] - u[2:Nx-1])

```

Removing the views and the dot each add to the computing time.

---

<div class="post-metadata">

**Author:** ![cortner](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cortner/32/204_2.png) [@cortner](https://discourse.julialang.org/u/cortner)\
**Post date:** [December 1, 2017, 2:25pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/6 "2017-12-01T14:25:48Z")

</div>

This is the key issue in matlab: cost of function calls

---

<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:** [December 1, 2017, 2:51pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/7 "2017-12-01T14:51:06Z")

</div>

Yup, these numbers are quite similar to the ones I’ve been seeing after writing tons of these things 😄. You even get the Julia beating C part as well. I’ll have to start sharing this.

---

<div class="post-metadata">

**Author:** ![ScottPJones](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/scottpjones/32/146_2.png) [@ScottPJones](https://discourse.julialang.org/u/ScottPJones)\
**Post date:** [December 1, 2017, 3:19pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/8 "2017-12-01T15:19:30Z")

</div>

> [@antoine-levitt](#):
>
> There’s still a penalty for writing vector code in julia (doesn’t do SIMD?)

Have you tried using the `@simd` macro?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 1, 2017, 4:16pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/9 "2017-12-01T16:16:47Z")

</div>

Not sure how to apply simd to a broadcast, isn’t it supposed to be before a loop?

@ChrisRackauckas how do you explain beating C? I don’t see what julia knows that C doesn’t here. Could it be that clang as a front-end does not produce as good a llvm IR as julia? gcc was about the same as julia in my example.

---

<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:** [December 1, 2017, 4:28pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/10 "2017-12-01T16:28:40Z")

</div>

> [@antoine-levitt](#):
>
> @ChrisRackauckas how do you explain beating C?

Here? No idea. Basically, our benchmarks routinely show we beat a lot of the C/Fortran methods by about 8x-10x. In this post:

> <https://stackoverflow.com/questions/47501844/julia-differentialequations-jl-speed/47507933#47507933>

I explain how most of that is explained by late binding (compiling the solver with the ODE function in that. A package which provides a solver can’t do this). However, there’s still this curious ~2x I don’t understand. In your example, you essentially do the “late binding” by writing it in there, but you also re-create the unexplained difference (at least when compiling with Clang). I said the numbers lined up with what I’ve seen (vs MATLAB/Python/etc as well), but there are some pieces of it that I don’t quite understand why they exist yet.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [December 1, 2017, 4:37pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/11 "2017-12-01T16:37:00Z")

</div>

I am a bit confused: there is only one source file for each language. How did you construct the table? Shouldn’t there be multiple versions (vectorized, loops)?

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [December 1, 2017, 4:46pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/12 "2017-12-01T16:46:35Z")

</div>

I think the various versions are derived by fiddling with the comments. Am I right?

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 1, 2017, 5:23pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/13 "2017-12-01T17:23:58Z")

</div>

Yes. I’m a very low-tech person 😉

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [December 1, 2017, 5:27pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/14 "2017-12-01T17:27:42Z")

</div>

If I may make a suggestion: generating multiple versions iin separate files might be preferable if it is to be generally useful to the readers of the discourse. Right now, at least for some of the codes, it is ambiguous how to derive any particular version. What to leave out and what to keep is obviously going to affect the timing (and possibly the correctness of the result as well).

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [December 1, 2017, 5:53pm UTC](https://discourse.julialang.org/t/simple-pde-benchmark/7430/15 "2017-12-01T17:53:19Z")

</div>

I just made this to figure out roughly what the timings were, uploaded this to github to save it in case I want to look at it in a year, and posted it here in case people were interested. Anybody is very welcome to take the code and put it in a benchmarking suite if they wish, but I have no interest in making it systematic.