matmul.jl

using Random
using BenchmarkTools

function matmul_python(C, A, B)
    for m in 1:C.rows
        for n in 1:C.cols
            for k in 1:A.cols
                C[m, n] += A[m, k] * B[k, n]
            end
        end
    end
end

struct Matrix{T}
    value::Array{T, 2}
    rows::Int
    cols::Int
end

function Matrix(value::Array{T, 2}, rows::Int, cols::Int) where T
    Matrix{T}(value, rows, cols)
end

function Base.getindex(m::Matrix{T}, i::Int, j::Int) where T
    m.value[i, j]
end

function Base.setindex!(m::Matrix{T}, value, i::Int, j::Int) where T
    m.value[i, j] = value
end

function benchmark_matmul_python(M, N, K)
    global A = Matrix(rand(Float64, M, K), M, K)
    global B = Matrix(rand(Float64, K, N), K, N)
    global C = Matrix(zeros(Float64, M, N), M, N)
    matmul_python(C, A, B)  # Perform matrix multiplication before benchmarking
    secs = @belapsed matmul_python(C, A, B)
    println(secs)
    gflops = (2 * M * N * K) / (secs * 1e9)
    println(gflops, " GFLOP/s")
    return gflops
end

python_gflops = benchmark_matmul_python(128, 128, 128)