# Comparing Numba and Julia for a complex matrix computation

**URL:** <https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703>\
**Category:** Performance\
**Tags:** question\
**Created:** [June 28, 2021, 3:43pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703 "2021-06-28T15:43:24Z")\
**Posts on this page:** 10\
**Page:** 2

<div class="post-metadata">

**Author:** ![brenhinkeller](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brenhinkeller/32/36648_2.png) [@brenhinkeller](https://discourse.julialang.org/u/brenhinkeller)\
**Post date:** [July 18, 2021, 8:45pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/21 "2021-07-18T20:45:21Z")

</div>

This seems to be another case where the best optimization advice for such numeric code is just to write simple loops and arrays, indexing element-wise, and then add LoopVectorization

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.julialang.org/u/jroon)\
**Post date:** [July 18, 2021, 8:51pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/22 "2021-07-18T20:51:44Z")

</div>

> [@brenhinkeller](#):
>
> It looks like @jroon’s simplified and allocation-free approach also makes it easy to add both SIMD vectorization and multithreading with the `@tturbo` macro from @elrod’s always awesome [LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl):
> 
> ```julia
> using LoopVectorization
> function computeSumHours5_lv(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
> height, width = size(height_mat)
> short_elevations = @view enlarged_size[y:y+height-1, x:x+width-1]
> result = 0.0
> 
> @tturbo for x2 in 1:width
> for y2 in 1:height
> overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
> # Note that we have to eliminate the branch to SIMD-vectorize:
> result += overshadowed * weights[angle_mat[y2, x2] ]
> end
> end
> 
> return result
> end
> 
> function computeAllLines5_lv(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
> for x in 1:size(org_sized, 2) - 1
> for y in 1:size(org_sized, 1) - 1
> shadow_time_hrs[x, y] = computeSumHours5_lv(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
> end
> end
> return shadow_time_hrs
> end
> 
> ```
> 
> Which gives us more than a 16x speedup on an 8-core CPU!

Wow I gotta learn about `@tturbo` it seems 😅

---

<div class="post-metadata">

**Author:** ![brenhinkeller](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brenhinkeller/32/36648_2.png) [@brenhinkeller](https://discourse.julialang.org/u/brenhinkeller)\
**Post date:** [July 18, 2021, 8:55pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/23 "2021-07-18T20:55:13Z")

</div>

Gotta go fast 😁

---

<div class="post-metadata">

**Author:** ![brenhinkeller](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/brenhinkeller/32/36648_2.png) [@brenhinkeller](https://discourse.julialang.org/u/brenhinkeller)\
**Post date:** [July 18, 2021, 9:07pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/24 "2021-07-18T21:07:59Z")

</div>

Basically, the magic starts with the `@turbo` macro (originally called `@avx`) which Chris Elrod wrote to take a `for` loop\* or broadcast and emit its own custom llvm that makes good use of your CPU’s vector registers and instructions (e.g. AVX-512 and friends). More recently, they also added `@tturbo`, (with an extra `t` for “threaded”) (originally called `@avxt`) which uses very low overhead threading from [Polyester.jl](https://github.com/JuliaSIMD/Polyester.jl) to also multithread at the same time.

\*Terms and conditions: must be reorderable, must not have branches, etc.

Some recent discussion of the SIMD part on Zulip: [Zulip](https://julialang.zulipchat.com/#narrow/stream/225542-helpdesk/topic/Project.20Euler/near/246350157)

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [July 21, 2021, 7:21pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/25 "2021-07-21T19:21:45Z")

</div>

Thanks @jroon for trying to simplify the code. Your made the code more simple and even faster. As for @brenhinkeller somehow I got the vague impression, that the `@tturbo` marco doesn’t scale so so good on some systems, but I am not completely sure yet.

I try to explain, what this example code should be an abstraction of. Also which mistake I made.

Given there is a matrix of of heights. The question is, low long is every place overshadowed. Seams like a limited submatrix of places is enough to calculate the overshadowing of a given point. I try to prepare/prepocess some matrices to simplify the actual preparation. Which is:

Given a maximum height difference the size of submatrix is static. I have prepared a matrix of height differences. If a point of the sub-matrix is higher than the height of the overshadowed place + the height difference, means it will actually overshadow the place.

Now, when I want to look at the time the place is overshadowed, this is depending on the angle segments. Like the shadowing time in the middle south is lower than the time in the west or east.

And now we come to the part of the code that is missing: Every shadow segment can throw that show only once. Example: A building throws a shadow on me, I don’t care about the shadow behind.

So what is still needed/missing is to split the old codes line:

`shadowed_segments[angle] = weights[angle] * overshadowed`

to have a set of angles and add the angle if overshadowed evaluates to true. Than only in the end for every angle in the set of angles the Time segment of the angle segments can be calculated.

My current implementation: is about six, seven times slower than, @jroon 's . Also `@turbo` does not not work on sets. So the current implementation is SUPER slow.

```nohighlight
function computeSumHours2!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    angle_set = BitSet(0)

    @inbounds for x2 ∈ 1:width
        for y2 ∈ 1:height
            angle = angle_mat[y2, x2]
            @fastmath overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            union!(angle_set, angle*overshadowed)
        end
    end

    result = zero(eltype(weights))
    @inbounds for elem ∈ setdiff!(angle_set, 0)
        @fastmath result += weights[elem]
    end

    return result
end

function computeAllLines5!(org_sized, enlarged_size, angle_mat, height_mat, shadow_time_hrs, weights)
    for x ∈ 1:size(org_sized, 2) - 1
        @threads for y ∈ 1:size(org_sized, 1) - 1
            shadow_time_hrs[x, y] = computeSumHours2!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x)
        end
    end
    return shadow_time_hrs
end

```

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [July 31, 2021, 12:10pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/26 "2021-07-31T12:10:12Z")

</div>

I am still trieing the improve the code for the additional restriction. The current code uses the improvements @jroon and @brenhinkeller added a while ago.  
And adds an additional loop afterwards. I tried different implementations so far, but everyone had a different setback I did not understand.

I wanted to ask you about the problems of the latest two versions.  
The latest version (the one with the fewest setbacks) is 25 % faster with `@turbo` than with `@tturbo` and even this version is 17 % slower than when I use @bounds instead. (Therefore I updated: LoopVectorization, didn’t change anything)

```julia
function computeSumHours3!(org_sized, enlarged_size, angle_mat, height_mat, weights, y, x, angle_set)
    height, width = size(height_mat)
    short_elevations = @view enlarged_size[y:y+height, x:x+width]
    angle_set .= zero(eltype(angle_set))

    for x2 ∈ 1:width
        @inbounds for y2 ∈ 1:height
            angle = angle_mat[y2, x2]
            @fastmath overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
            angle_set[angle] |= overshadowed
        end
    end

    result = zero(eltype(weights));
    @inbounds for elem ∈ 1:length(weights)
        @fastmath result += weights[elem] * angle_set[elem];
    end

    return result
end

```

I added the parameter:  
`angle_set = zeros(Bool, length(weights))`

The version I tried before:

```nohighlight
        0 height, width = size(height_mat)
        0 short_elevations = @view enlarged_size[y:y+height, x:x+width]
        0 angle_set .= zero(eltype(angle_set))
        - 
        0 for x2 ∈ 1:width
        0 @turbo for y2 ∈ 1:height
        - angle = angle_mat[y2, x2]
        - @fastmath overshadowed = (short_elevations[y2, x2] - org_sized[y, x]) > height_mat[y2, x2]
        0 angle_set[angle] |= overshadowed
        - end
        - end
        - 
        - result = zero(eltype(weights))
     7456 angle_list = findall(angle_set .== true)
        0 @inbounds for elem ∈ angle_list
        0 @fastmath result += weights[elem]
        - end
        - 
        0 return result
        - end

```

It fails when the second loop uses `@turbo` because the iteration isn’t of the form: 1:length(…) and it uses tons of memory for `findall(angle\_set .== true). Which I think is very strange, because the array has only 361 entries.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [July 31, 2021, 12:13pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/27 "2021-07-31T12:13:03Z")

</div>

> [@Ordens\_Ritter](#):
>
> `findall(angle\_set .== true). Which I think is very strange, because the array has only 361 entries.

isn’t this just `findall(angle_set)`?.. ` .== true` will allocate a new array for no reason.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [July 31, 2021, 1:39pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/28 "2021-07-31T13:39:27Z")

</div>

> [@Ordens\_Ritter](#):
>
> ```julia
> - result = zero(eltype(weights))
> 7456 angle_list = findall(angle_set .== true)
> 0 @inbounds for elem ∈ angle_list
> 0 @fastmath result += weights[elem]
> - end
> 
> ```

Could also be written like this, which should be easier for the auto vectorizer to understand:

```julia
result = zero(eltype(weights))
for elem in eachindex(angle_set)
    result += weights[elem] * angle_set[elem] # true/false is a number in julia
end

```

Looks like a classic case for `mapreduce` though (or threaded equivalents).

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [July 31, 2021, 1:43pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/29 "2021-07-31T13:43:06Z")

</div>

superbe, reduces the memory usage somewhat from 7,4 kb to 2.95 kb. Still I don’t understand why, there is so much memory usage left. Does’t it mean, that there is something in between allocated saved as Int64?

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [July 31, 2021, 1:50pm UTC](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703/30 "2021-07-31T13:50:07Z")

</div>

> [@Ordens\_Ritter](#):
>
> 361

If these were 64-bit numbers, data alone takes

```julia
julia> 361 * 8
2888 #bytes

```

and there’s some tag data:

```julia
julia> Base.summarysize(rand(361))
2928

```

So I don’t think allocation on this order of magnitude is unacceptable.

[Previous page](https://discourse.julialang.org/t/comparing-numba-and-julia-for-a-complex-matrix-computation/63703.md?page=1)