tensoroperations.jl

using TensorOperations, BenchmarkTools
n = 20
Γ1 = rand(n, n, n);
Γ2 = rand(n, n, n);

f(Γ1, Γ2) = @tensor Γ[α, β, i, k] := Γ1[α, i, j] * Γ2[β, j, k]

function g(Γ1, Γ2)
    n = size(Γ1, 1)
    Γ = zeros(n, n, n, n)
    for α in 1:n, β in 1:n, i in 1:n, k in 1:n
        Γ[α, β, i, k] = Γ1[α, i, :]'Γ2[β, :, k]
    end
    return Γ
end

function h(Γ1, Γ2)
    n = size(Γ1, 1)
    Γ = zeros(n, n, n, n)
    @inbounds @views for α in 1:n, β in 1:n, i in 1:n, k in 1:n
        Γ[α, β, i, k] = Γ1[α, i, :]'Γ2[β, :, k]
    end
    return Γ
end

@btime f($Γ1, $Γ2); # 403.300 μs (194 allocations: 1.24 MiB)
@btime g($Γ1, $Γ2); # 26.453 ms (320003 allocations: 74.46 MiB)
@btime h($Γ1, $Γ2); # 3.143 ms (3 allocations: 1.22 MiB)