# Performance of a transpose compared to c and cache friendliness

**URL:** <https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883>\
**Category:** Performance\
**Created:** [March 10, 2023, 7:00pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883 "2023-03-10T19:00:20Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![sneakyturtle](https://avatars.discourse-cdn.com/v4/letter/s/d6d6ee/32.png) [@sneakyturtle](https://discourse.julialang.org/u/sneakyturtle)\
**Post date:** [March 10, 2023, 7:00pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/1 "2023-03-10T19:00:20Z")

</div>

In a quest to learn good HPC principles, I am trying to understand how to implement cache friendly loops in Julia.

I found this code [loop\_blocking](https://github.com/ZiCog/loop_blocking) that compares performance of loop blocking techniques between c and Rust.  
I tried implementing the same thing in Julia, but I’m confused by the performance difference I’m getting. I’m writing here to seek help to understand the performance of this code.

```julia
module LoopBlocking
const MAX::Int64 = 8192

function Fill(A::Array{Int64}, B::Array{Int64})
    count::Int64 = 0
    @inbounds for i = 1:MAX
        for j = 1:MAX
            A[i, j] = count
            B[i, j] = -count
            count += 1
        end
    end
end

function Add(A::Array{Int64}, B::Array{Int64})
    for i = 1:MAX, j = 1:MAX
        @inbounds A[i, j] += B[j, i]
    end
    return nothing
end

function Add_2(A::Array{Int64}, B::Array{Int64})
    @inbounds A .+= transpose(B)
    return nothing
end

const BLOCK_SIZE::Int64 = 32
function Add_1(A::Array{Int64}, B::Array{Int64})
    for i = 1:BLOCK_SIZE:MAX, j = 1:BLOCK_SIZE:MAX
        for ii = i:i+BLOCK_SIZE-1, jj = j:j+BLOCK_SIZE-1
            @inbounds A[ii, jj] += B[jj, ii]
        end
    end
    return nothing
end

function Add_3(A::Array{Int64}, B::Array{Int64})
    for i = 1:BLOCK_SIZE:MAX, j = 1:BLOCK_SIZE:MAX
        @inbounds @views A[i:i+BLOCK_SIZE-1, j:j+BLOCK_SIZE-1] .+= transpose(B[j:j+BLOCK_SIZE-1, i:i+BLOCK_SIZE-1])
    end
    return nothing
end

function main()
    A::Array{Int64} = Array{Int64}(undef, MAX, MAX)
    B::Array{Int64} = Array{Int64}(undef, MAX, MAX)

    Fill(A, B)
    Add(A, B)

    println("Simple loop")
    @time Fill(A, B)
    @time Add(A, B)

    Correct = copy(A)

    Fill(A, B)
    Add_1(A, B)
    
    println("Loop blocking")
    @time Fill(A, B)
    @time Add_1(A, B)

    @assert A == Correct

    Fill(A, B)
    Add_2(A, B)

    println("Built-in broadcasting and transpose")
    @time Fill(A, B)
    @time Add_2(A, B)

    @assert A == Correct

    Fill(A, B)
    Add_3(A, B)
    
    println("Loop blocking with views")
    @time Fill(A, B)
    @time Add_3(A, B)

    @assert A == Correct
end

end

```

These are the results I get

```bash
Simple loop
  1.253246 seconds
  0.841699 seconds
Loop blocking
  1.238872 seconds
  0.458186 seconds
Built-in broadcasting and transpose
  1.165040 seconds
  0.461195 seconds
Loop blocking with views
  1.266697 seconds
  0.286722 seconds

```

And here is typical c performance

```bash
MAX: 8192
BLOCK_SIZE: 32
transpose_0: 528ms
transpose_1: 216ms

```

It seems that broadcasting leads to different performance than simple loops and is the closes to c performance.  
Can someone help me understand why simple for loops cannot achieve the same performance as broadcasting?

P.S.: I’m not entirely sure the last function is correct, but the assert seems to work.

---

<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:** [March 10, 2023, 7:15pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/2 "2023-03-10T19:15:29Z")

</div>

My first guess without looking closely is you’re seeing the difference between row and column major arrays.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 10, 2023, 7:24pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/3 "2023-03-10T19:24:45Z")

</div>

Not specific to your question but whenever you benchmark Julia code, be sure to use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) and follow their instructions, this will yield more accurate results for comparison

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 10, 2023, 7:31pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/4 "2023-03-10T19:31:40Z")

</div>

> [@sneakyturtle](#):
>
> ```julia
> @inbounds for i = 1:MAX
> for j = 1:MAX
> 
> ```

Your loops are in the [wrong order](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major).

~~Also, there is not much point in blocking loops like this that make only a single pass over the array, in order, e.g. to fill it.~~ You only want to block in order to increase temporal locality (e.g. in matrix multiplication, see e.g. [this Julia notebook](https://github.com/mitmath/18335/blob/spring20/notes/Memory-and-Matrices.ipynb)) and/or spatial locality (e.g. for [matrix transposition](https://discourse.julialang.org/t/function-on-matrix-transpose-and-performance/20068/4)). _Update_: sorry, I missed that you are computing A+B^T, see below.

---

<div class="post-metadata">

**Author:** ![sneakyturtle](https://avatars.discourse-cdn.com/v4/letter/s/d6d6ee/32.png) [@sneakyturtle](https://discourse.julialang.org/u/sneakyturtle)\
**Post date:** [March 11, 2023, 1:31pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/5 "2023-03-11T13:31:00Z")

</div>

Thanks for the reply, changing the order gives the same performance as broadcasting (although, both are slightly slower than c ☹).

I thought that since one matrix is transpose, there shouldn’t matter which way I do the loops. I guess it’s better to follow the matrix I assign to than the one I read?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 11, 2023, 3:18pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/6 "2023-03-11T15:18:22Z")

</div>

> [@sneakyturtle](#):
>
> I thought that since one matrix is transpose, there shouldn’t matter which way I do the loops. I guess it’s better to follow the matrix I assign to than the one I read?

Oh, I was only looking at your `Fill` function, and missed that you swapped the indices in subsequent functions. If you do A + B^T, then indeed there is a spatial-locality (cache-line) benefit to blocking (or, alternatively, a cache-oblivious recursive strategy), but I would still order the loops for `A` (which you both read and write) rather than `B` (which you only read).

---

<div class="post-metadata">

**Author:** ![mstewart](https://avatars.discourse-cdn.com/v4/letter/m/b5a626/32.png) [@mstewart](https://discourse.julialang.org/u/mstewart)\
**Post date:** [March 11, 2023, 3:33pm UTC](https://discourse.julialang.org/t/performance-of-a-transpose-compared-to-c-and-cache-friendliness/95883/7 "2023-03-11T15:33:01Z")

</div>

> I would still order the loops for A (which you both read and write) rather than B (which you only read).

That is, in fact, what makes the difference between `Add_3` and `Add_1` on my machine. I used to have some C code in which I used a lot of this sort of blocking, but I don’t think I worried about distinguishing reads or writes. It’s not really all that small difference. I probably left some easy performance on the table.