Бенчмаркинг (python vs. c++ using BLAS) и (numpy)

Я хотел бы написать программу, которая широко использует функции линейной алгебры BLAS и LAPACK. Поскольку производительность является проблемой, я сделал несколько бенчмарков и хотел бы знать, является ли подход, который я использовал, законным.

У меня есть, так сказать, три участника, и я хочу проверить их производительность с помощью простого умножения матрицы на матрицу. Участники конкурса:

1. Numpy, использующий только функциональность dot.
2. Python, вызов функций BLAS через общий объект.
3. C++, вызов функций BLAS через общий объект.

Сценарий

Я реализовал умножение матрицы на матрицу для разных размерностей i i имеет размерность от 5 до 500 с шагом 5, а матрицы m1 и m2 заданы следующим образом:

m1 = numpy.random.rand(i,i).astype(numpy.float32)
m2 = numpy.random.rand(i,i).astype(numpy.float32)

1. Numpy

Используемый код выглядит следующим образом:

tNumpy = timeit.Timer("numpy.dot(m1, m2)", "import numpy; from __main__ import m1, m2")
rNumpy.append((i, tNumpy.repeat(20, 1)))

2. Python, вызов BLAS через общий объект

С помощью функции

_blaslib = ctypes.cdll.LoadLibrary("libblas.so")
def Mul(m1, m2, i, r):

    no_trans = c_char("n")
    n = c_int(i)
    one = c_float(1.0)
    zero = c_float(0.0)

    _blaslib.sgemm_(byref(no_trans), byref(no_trans), byref(n), byref(n), byref(n), 
            byref(one), m1.ctypes.data_as(ctypes.c_void_p), byref(n), 
            m2.ctypes.data_as(ctypes.c_void_p), byref(n), byref(zero), 
            r.ctypes.data_as(ctypes.c_void_p), byref(n))

тестовый код выглядит так:

r = numpy.zeros((i,i), numpy.float32)
tBlas = timeit.Timer("Mul(m1, m2, i, r)", "import numpy; from __main__ import i, m1, m2, r, Mul")
rBlas.append((i, tBlas.repeat(20, 1)))

3. c++, вызов BLAS через общий объект

Теперь код c++, естественно, немного длиннее, поэтому я сокращаю информацию до минимума.
Я загружаю функцию с помощью

void* handle = dlopen("libblas.so", RTLD_LAZY);
void* Func = dlsym(handle, "sgemm_");

Я измеряю время с помощью gettimeofday вот так:

gettimeofday(&start, NULL);
f(&no_trans, &no_trans, &dim, &dim, &dim, &one, A, &dim, B, &dim, &zero, Return, &dim);
gettimeofday(&end, NULL);
dTimes[j] = CalcTime(start, end);

где j - цикл, выполняющийся 20 раз. Я вычисляю прошедшее время с помощью

double CalcTime(timeval start, timeval end)
{
double factor = 1000000;
return (((double)end.tv_sec) * factor + ((double)end.tv_usec) - (((double)start.tv_sec) * factor + ((double)start.tv_usec))) / factor;
}

Results

Результат показан на графике ниже:

Вопросы

1. Считаете ли вы, что мой подход справедлив, или есть какие-то ненужные накладные расходы, которых я могу избежать?
2. Ожидали ли вы, что результат покажет такое огромное расхождение между подходами на c++ и python? Оба используют общие объекты для своих вычислений.
3. Поскольку я предпочел бы использовать python для своей программы, что я мог бы сделать для увеличения производительности при вызове процедур BLAS или LAPACK?

Скачать

Полный бенчмарк можно скачать здесь. (J.F. Sebastian сделал эту ссылку возможной^^)