Gostaria de escrever um programa que faça uso extensivo das funcionalidades de álgebra linear BLAS e LAPACK. Como o desempenho é um problema, fiz alguns testes e gostaria de saber se a abordagem adotada é legítim

Tenho, por assim dizer, três concorrentes e quero testar seu desempenho com uma simples multiplicação matriz-matriz. Os concorrentes são:

Numpy, utilizando apenas a funcionalidade dedot. Python, chamando as funcionalidades do BLAS através de um objeto compartilhad C ++, chamando as funcionalidades do BLAS através de um objeto compartilhadCenári

mplementei uma multiplicação matriz-matriz para diferentes dimensõei. i corre de 5 a 500 com um incremento de 5 e as matrizesm1 em2 são configurados assim:

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

1. Numpy

O código usado se parece com o seguinte:

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

2. Python, chamando BLAS através de um objeto compartilhado

Com a função

_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))

o código de teste é assim:

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 ++, chamando BLAS através de um objeto compartilhado

gora, naturalmente, o código c ++ é um pouco mais longo, então reduzo as informações ao mínim
Carrego a função com

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

Meu tempo comgettimeofday como isso

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);

Ondej é um loop executando 20 vezes. Calculo o tempo passado com

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;
}

Resultado

O resultado é mostrado no gráfico abaixo:

Questõe Você acha que minha abordagem é justa ou existem algumas despesas indiretas desnecessárias que posso evita Você esperaria que o resultado mostrasse uma discrepância tão grande entre a abordagem c ++ e python? Ambos estão usando objetos compartilhados para seus cálculo Desde que eu prefiro usar python no meu programa, o que posso fazer para aumentar o desempenho ao chamar as rotinas BLAS ou LAPACBaixa

O benchmark completo pode ser baixadoaqu. (J.F. Sebastian tornou possível esse link ^^)