simplemul.jl

function simplemul!((β₁, Y, β₂), α₁, A, α₂, B, α₃)
    @assert size(Y, 1) == size(A, 1)
    @assert size(Y, 2) == size(B, 2)
    @assert size(A, 2) == size(B, 1)
    @inbounds for i in 1:size(A, 1), j = 1:size(B, 2)
        acc = zero(α₁ * A[i, 1] * α₂ * B[1, j] * α₃ +
                   α₁ * A[i, 1] * α₂ * B[1, j] * α₃)
        for k = 1:size(A, 2)
            acc += A[i, k] * α₂ * B[k, j]
        end
        Y[i, j] = α₁ * acc * α₃ + β₁ * Y[i, j] * β₂
    end
    return Y
end

function simplemul!((Y, β), A, B, α)
    @assert size(Y, 1) == size(A, 1)
    @assert size(Y, 2) == size(B, 2)
    @assert size(A, 2) == size(B, 1)
    @inbounds for i in 1:size(A, 1), j = 1:size(B, 2)
        acc = zero(A[i, 1] * B[1, j] * α +
                   A[i, 1] * B[1, j] * α)
        for k = 1:size(A, 2)
            acc += A[i, k] * B[k, j]
        end
        Y[i, j] = acc * α + Y[i, j] * β
    end
    return Y
end

fullmul!(Y, A, B) = simplemul!((false, Y, false), true, A, true, B, true)
minmul!(Y, A, B) = simplemul!((Y, false), A, B, true)

using LinearAlgebra
k = 50
n = 50
A = [randn(k, k) for _ in 1:n, _ in 1:n]
B = [randn(k, k) for _ in 1:n]
Y = [zeros(k, k) for _ in 1:n]
@assert mul!(copy(Y), A, B) == fullmul!(copy(Y), A, B) == minmul!(copy(Y), A, B)

using BenchmarkTools
@btime mul!($Y, $A, $B)     # 63.845 ms (10400 allocations: 99.74 MiB)
@btime fullmul!($Y, $A, $B) # 80.963 ms (16501 allocations: 158.24 MiB)
@btime minmul!($Y, $A, $B)  # 64.017 ms (10901 allocations: 104.53 MiB)
nothing