# Julia matrix-multiplication performance

**URL:** https://discourse.julialang.org/t/julia-matrix-multiplication-performance/55175
**Category:** Performance
**Tags:** linearalgebra
**Created:** [February 12, 2021, 5:43pm UTC](https://discourse.julialang.org/t/julia-matrix-multiplication-performance/55175 "2021-02-12T17:43:48Z")
**Posts on this page:** 1
**Showing post:** 12

<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: [February 12, 2021, 10:12pm UTC](https://discourse.julialang.org/t/julia-matrix-multiplication-performance/55175/12 "2021-02-12T22:12:57Z")

</div>

For example, here is a little recursive implementation of a [cache-oblivious matrix multiplication](https://en.wikipedia.org/wiki/Cache-oblivious_algorithm) that stays within a factor of 2 of the single-threaded OpenBLAS performance on my laptop up to 3000×3000 matrices:

```julia
function add_matmul_rec!(m,n,p, i0,j0,k0, C,A,B)
    if m+n+p <= 256 # base case: naive matmult for sufficiently large matrices
        @avx for i = 1:m, k = 1:p
            c = zero(eltype(C))
            for j = 1:n
                @inbounds c += A[i0+i,j0+j] * B[j0+j,k0+k]
            end
            @inbounds C[i0+i,k0+k] += c
        end
    else
        m2 = m ÷ 2; n2 = n ÷ 2; p2 = p ÷ 2
        add_matmul_rec!(m2, n2, p2, i0, j0, k0, C, A, B)
        
        add_matmul_rec!(m-m2, n2, p2, i0+m2, j0, k0, C, A, B)
        add_matmul_rec!(m2, n-n2, p2, i0, j0+n2, k0, C, A, B)
        add_matmul_rec!(m2, n2, p-p2, i0, j0, k0+p2, C, A, B)
        
        add_matmul_rec!(m-m2, n-n2, p2, i0+m2, j0+n2, k0, C, A, B)
        add_matmul_rec!(m2, n-n2, p-p2, i0, j0+n2, k0+p2, C, A, B)
        add_matmul_rec!(m-m2, n2, p-p2, i0+m2, j0, k0+p2, C, A, B)
        
        add_matmul_rec!(m-m2, n-n2, p-p2, i0+m2, j0+n2, k0+p2, C, A, B)
    end
    return C
end

function matmul_rec!(C, A, B)
    m,n = size(A)
    n,p = size(B)
    size(C) == (m,p) || error("incorrect dimensions ", size(C), " ≠ $m × $p")
    fill!(C, zero(eltype(C)))
    return add_matmul_rec!(m,n,p, 0,0,0, C,A,B)
end

matmul_rec(A, B) = matmul_rec!(Array{promote_type(eltype(A), eltype(B))}(undef,
                                     size(A,1), size(B,2)),
                               A, B)

```

In comparison, Octavian.jl is ~100× more code, and OpenBLAS is \> 1000× more code (and only supports 4 arithmetic types, whereas this code is type-generic).

Here larger numbers = faster (vertical axis is rate of floating-point arithmetic operations), timings are for double precision, “naive” is a 3-loop Julia version with `@avx`, and BLAS is `mul!` with OpenBLAS:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/0/90c7668a840fd494fb297469f576df0b74cb09ad.jpeg)

---

_[View the full topic](https://discourse.julialang.org/t/julia-matrix-multiplication-performance/55175)._