perf.jl

function test(f)
    tmin = 1e10
    for i in 1:5
        t = time()
        f()
        t = time() - t
        if t < tmin
            tmin = t
        end
    end
    return tmin
end


function get_test_func(N, reps, op; A_complex=false, B_complex=false)
    A = float(rand(N, N))
    if A_complex
        A = complex(A)
    end

    B = float(rand(N, N))
    if B_complex
        B = complex(B)
    end

    @assert eltype(A) == (A_complex ? Complex{Float64} : Float64)
    @assert eltype(B) == (B_complex ? Complex{Float64} : Float64)

    if op == :dot
        @assert eltype(A * B) == (A_complex ? Complex{Float64} : Float64)
        function f()
            for i in 1:reps
                C = A * B
            end
        end
    elseif op == :add
        @assert eltype(A .+ B) == (A_complex ? Complex{Float64} : Float64)
        function f()
            for i in 1:reps
                C = A .+ B
            end
        end
    end

    return f
end


for op in (:dot, :add)
    for (A_complex, B_complex) in ((false, false), (true, false), (true, true))
        ac = A_complex ? "complex" : "float"
        bc = B_complex ? "complex" : "float"
        for N in (10, 100, 1000)
            if op == :add
                reps = int(1e7 / N^2)
            elseif op == :dot
                reps = int(2e6 / N^2)
            end
            f = get_test_func(N, reps, op, A_complex=A_complex, B_complex=B_complex)
            t = test(f)
            println("$op($ac, $bc), $N x $N: $t s")
        end
    end
end

perf.py

import time
import numpy as np

def test(f):
    tmin = 1e10
    for i in range(5):
        t = time.time()
        f()
        t = time.time() - t
        if t < tmin:
            tmin = t
    return tmin


def get_test_func(N, reps, op, A_complex=False, B_complex=False):
    A = np.random.rand(N, N)
    if A_complex:
        A = A.astype(np.complex128)

    B = np.random.rand(N, N)
    if B_complex:
        B = B.astype(np.complex128)

    assert A.dtype == np.complex128 if A_complex else np.float64
    assert B.dtype == np.complex128 if B_complex else np.float64

    if op == 'dot':
        assert (A.dot(B)).dtype == np.complex128 if A_complex else np.float64
        def f():
            for i in range(reps):
                C = A.dot(B)
    elif op == 'add':
        assert (A + B).dtype == np.complex128 if A_complex else np.float64
        def f():
            for i in range(reps):
                C = A + B

    return f


for op in ('dot', 'add'):
    for A_complex, B_complex in ((False, False), (True, False), (True, True)):
        ac = "complex" if A_complex else "float"
        bc = "complex" if B_complex else "float"
        for N in (10, 100, 1000):
            if op == 'add':
                reps = int(1e7 / N**2)
            elif op == 'dot':
                reps = int(2e6 / N**2)
            f = get_test_func(N, reps, op, A_complex=A_complex, B_complex=B_complex)
            t = test(f)
            print("{op}({ac}, {bc}), {N} x {N}: {t} s".format(op=op, ac=ac, bc=bc, N=N, t=t))