# Performance optimisation: Julia vs C++

**URL:** <https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689>\
**Category:** Performance\
**Tags:** unrolling\
**Created:** [April 20, 2021, 6:50pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689 "2021-04-20T18:50:05Z")\
**Posts on this page:** 12\
**Page:** 2

<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:** [April 20, 2021, 7:45pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/21 "2021-04-20T19:45:41Z")

</div>

> [@Kolaru](#):
>
> Simply defining `const L = 64` outside of the `main` function gave me a x10 speed up.

Did you do something else? (compilation options, etc?) I cannot reproduce that.

Ok, with everything constant (L, Nx, Ny), I also get the 6x speedup.

---

<div class="post-metadata">

**Author:** ![Kolaru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kolaru/32/4574_2.png) [@Kolaru](https://discourse.julialang.org/u/Kolaru)\
**Post date:** [April 20, 2021, 7:48pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/22 "2021-04-20T19:48:40Z")

</div>

Yeah sorry I made a stupid mistake (changed the number of iteration when converting it to a constant x\_x). I edited my message. I still get a significant x6 speed up by declaring everything to be constant.

---

<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:** [April 20, 2021, 8:27pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/23 "2021-04-20T20:27:12Z")

</div>

Inspired by your observation, @Kolaru , here is a version that runs in `404ms` (vs. 2s from the original one) without having to make those variables constant. Did that (hopefully correctly) by putting the sizes in the type signature of the spin configuration:

(running with `--check-bounds=no` makes it run in `313ms`).

```julia
using Random
using BenchmarkTools

const rng = MersenneTwister()

struct SpinConf{Nx,Ny,T}
  m::Matrix{T}
end
Base.getindex(M::SpinConf,i,j) = @inbounds M.m[i,j]
Base.setindex!(M::SpinConf,x,i,j) = @inbounds (M.m[i,j] = x) 
Base.sum(M::SpinConf) = sum(M.m)

function get(spin_conf::SpinConf{Nx,Ny,T}, x, y) where {Nx,Ny,T}
    return @inbounds spin_conf[(x-1+Nx)%Nx + 1, (y-1+Ny)%Ny + 1]
end

function deltaE(spin_conf, i, j)
    return 2.0*get(spin_conf,i,j) * (get(spin_conf,i-1,j) + get(spin_conf,i+1,j) +
                                     get(spin_conf,i,j+1) + get(spin_conf,i,j-1))
end

function energy(spin_conf::SpinConf{Nx,Ny,T}) where {Nx,Ny,T}
    res = 0
    for i = 1:Nx
        @simd for j = 1:Ny
            res += -get(spin_conf,i,j)*(get(spin_conf,i+1,j) + get(spin_conf,i,j+1))
        end
    end
    return res
end

magnetism(spin_conf) = sum(spin_conf)

function sweep!(rec, spin_conf::SpinConf{Nx,Ny,T}, beta, M) where {Nx,Ny,T}
    @inbounds @simd for i = 1:M
        x = rand(rng, 1:Nx)
        y = rand(rng, 1:Ny)
        r = rand(rng)
        dE = deltaE(spin_conf, x, y)
        if exp(-beta*dE) > r
            spin_conf[x,y] *= -1
        end
        rec[i] = energy(spin_conf)
    end
    return rec
end

function main()
    beta = 0.2
    L = 64
    total_iteration = 100_000

    spin_conf = SpinConf{L,L,Int64}(rand(rng, (-1,1), L, L))
    rec = Vector{Int64}(undef,total_iteration)

    sweep!(rec, spin_conf, beta, total_iteration)

    s = sum(last(rec, 1000))/1000
    s /= L^2
    return nothing
end

@btime main()

```

obs: with `@inbounds` in these lines, no need to run without bounds check (~310ms):

```julia
Base.getindex(M::SpinConf,i,j) = @inbounds M.m[i,j]
Base.setindex!(M::SpinConf,x,i,j) = @inbounds (M.m[i,j] = x) 

```

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [April 20, 2021, 8:48pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/24 "2021-04-20T20:48:04Z")

</div>

Interesting that the type parameters for Nx Ny help. As you say there are still some boundschecks, `Base.@propagate_inbounds Base.getindex(M::SpinConf,i,j) = M.m[i,j]` helps, then 285.629 ms on my computer. Down to 243.265 ms if you switch the loop order in `energy`.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [April 20, 2021, 8:52pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/25 "2021-04-20T20:52:09Z")

</div>

Some timing updates on my machine (i.e. comparable to the OP)

**`--check-bounds=no` + making L/Nx/Ny `const` (or using a `SpinConf` struct)**

```julia
  295.933 ms (5 allocations: 821.34 KiB)

```

**`--check-bounds=no` + making L/Nx/Ny `const` (or using a `SpinConf` struct) + `@unroll`**

```julia
  217.009 ms (5 allocations: 821.34 KiB)

```

**`--check-bounds=no` + making L/Nx/Ny `const` (or using a `SpinConf` struct) + breaking up the square (@mcabbott )**

```julia
  158.804 ms (5 allocations: 821.34 KiB)

```

**`--check-bounds=no` + making L/Nx/Ny `const` (or using a `SpinConf` struct) + breaking up the square (@mcabbott ) + `@unroll`**

```julia
  154.410 ms (5 allocations: 821.34 KiB)

```

---

<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:** [April 20, 2021, 9:42pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/26 "2021-04-20T21:42:48Z")

</div>

All that shows how hard is to perform “fair” benchmarking. In this case the C++ code has a constraint on the parameters that the Julia code didn’t and the made a huge difference. The cool thing is that the Julia code now can be run efficiently for different sizes of the array, as the functions will get specialized versions for each size, even if run from the same “main”.

My first attempt was to use a static matrix, but compilation times became way too high.

---

<div class="post-metadata">

**Author:** ![sairus7](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sairus7/32/10816_2.png) [@sairus7](https://discourse.julialang.org/u/sairus7)\
**Post date:** [April 20, 2021, 9:51pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/27 "2021-04-20T21:51:09Z")

</div>

Change `for i = 1:Nx @simd for j = 1:Ny` order from `i, j` into `j, i` ???

Also, you can try making `spin_conf` an array of static arrays, so the loop should unroll.

UPD: Also, since `Nx` constant is power of two, you can change  
from `(x-1+Nx)%Nx + 1`  
into `(x-1+Nx)&(Nx-1) + 1`  
where `(Nx-1)` is also a constant.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [April 20, 2021, 11:45pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/28 "2021-04-20T23:45:14Z")

</div>

I’m not sure whether this is relevant to your apples-to-apples C++ vs julia comparison, but…

You don’t need to recompute the energy at every step. You already compute the change in energy for each ising flip, so just keep a running tally:

```julia
julia> function sweep!(rec, spin_conf, beta, M)
           (Nx, Ny) = size(spin_conf)
           E = Float64(energy(spin_conf))
           @inbounds @simd for i = 1:M
               x = rand(rng, 1:Nx)
               y = rand(rng, 1:Ny)
               r = rand(rng)
               dE = deltaE(spin_conf, x, y)
               if exp(-beta*dE) > r
                   spin_conf[x,y] *= -1
                   E += dE
               end
               rec[i] = E
           end
           energy(spin_conf)!= E && error("mismatch -- expected energy $(E), found $(energy(spin_conf))")
           return rec
       end
julia> @btime main()
  3.458 ms (5 allocations: 821.34 KiB)

```

Typically you should first think about what you compute, and how you do so. Then you should lay out your data (and the Ising config cries for a BitMatrix instead of an integer matrix!). Then at last you should care about `@inbounds` and loop vectorization.

(also, the `@simd` in front of the above loop is magical thinking – how is the loop possibly supposed to vectorize? `@simd` only makes sense if subsequent iterations are independent, it is not magical performance pixy dust)

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 21, 2021, 2:40am UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/29 "2021-04-21T02:40:44Z")

</div>

> [@foobar\_lv2](#):
>
> Ising config cries for a BitMatrix instead of an integer matrix!)

It’s difficult to parallelize writes to a BitMatrix, even if the problem is embarrassingly parallel across indices, since many indices share the same memory location, and I guess there is often no atomic operation to set/clear a bit.

Relatedly, there may be situations where the `BitArray` is slower, again because of the bit-level operations that need to be performed. That’s all to say, it isn’t so obvious to me that going with a BitArray is beneficial.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [April 21, 2021, 8:04am UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/30 "2021-04-21T08:04:01Z")

</div>

> [@goretkin](#):
>
> That’s all to say, it isn’t so obvious to me that going with a BitArray is beneficial.

The good part of BitMatrix for this kind of ising model is not the denser layout. The point is: Given two +/- 1 vectors `u, v` of length 64, you can compute `dot(u, v) = 64 - 2 * count_ones(xor(u_chunk, v_chunk))`, where `u_chunk::UInt64` is the bit-representation of the vector (true for +1, false for -1). In other words, you get 64-wide SIMD for free, without even using fancy vector instructions of your CPU.

Out of this, you can build a really fast `energy` implementation.

Depending on the problem, you may need to pad your matrix, possibly store both the matrix and its transpose, etc.

In other words, you don’t use BitMatrix as a black-box replacement for Matrix; instead you directly operate on the `bitmatrix.chunks`, and use the “BitArray” interface only for convenient low-effort slow operations that are either hard to write using the raw memory representation, or are outside of the hot path.

BitArray is a toolkit, not a black-box drop-in replacement for `Array{Bool}`.

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [April 21, 2021, 8:45am UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/31 "2021-04-21T08:45:59Z")

</div>

Thanks for your comments! Indeed, I’m aware of faster 2D Ising model implementation strategies. AFAIK, eventually, the high performance variants are mostly limited by the random number generation. The fact that this is not the case for the implementation here speaks for itself 🙂 But, for this challenge / discussion, moving to BitArrays or changing the algorithm (summing up the `dE`s) would have to be done on the C++ side as well, which would give us a new timing to chase.

(Just to be clear, in a “real application”, I do agree that algorithmic improvements should come first, of course!)

---

<div class="post-metadata">

**Author:** ![goretkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goretkin/32/167_2.png) [@goretkin](https://discourse.julialang.org/u/goretkin)\
**Post date:** [April 21, 2021, 2:55pm UTC](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689/32 "2021-04-21T14:55:24Z")

</div>

Gotcha. Can the updates also be performed in terms of the underlying representation?

```julia
dE = deltaE(spin_conf, x, y)
if exp(-beta*dE) > r
    spin_conf[x,y] *= -1
    E += dE
end

```

I know we’re off topic, but I sure would like to see more of this 64-wide SIMD bit stuff. I played around with some different layouts here:

> **[GitHub - goretkin/BitBlockArrays.jl](https://github.com/goretkin/BitBlockArrays.jl)**
>
> Contribute to goretkin/BitBlockArrays.jl development by creating an account on GitHub.

includes e.g. transposing an 8-by-8 array:

> <https://github.com/goretkin/BitBlockArrays.jl/blob/dcff1d571d8ad09df64401fe3d59b93cb72c8b0c/src/bitblockarray.jl#L43-L63>

This talk has many more examples:

[![](https://global.discourse-cdn.com/julialang/original/3X/b/8/b8406e32ef46083743c2d031358fe0460e46e06b.jpeg "Dyalog16: A Compendium of SIMD Boolean Array Algorithms in APL") ](https://www.youtube.com/watch?v=G2g13YKjWes)

[Previous page](https://discourse.julialang.org/t/performance-optimisation-julia-vs-c/59689.md?page=1)