python - Resultados decepcionantes en pyCUDA benchmark para computación a distancia entre N puntos

install pycuda

Question

La siguiente secuencia de comandos se configuró para fines de referencia. Calcula la distancia entre N puntos usando una norma euclidiana L2. Se implementan tres rutinas diferentes:

1. Solución de alto nivel que utiliza la función scipy.spatial.distance.pdist .
2. Solución de bajo nivel OpenMP powered scipy.weave.inline .
3. la solución GPGPU impulsada por pyCUDA.

Estos son los resultados de referencia en un i5-3470 (16 GB de RAM) con una GTX660 (2 GB de RAM):

 ------------
Scipy Pdist
Execution time: 3.01975 s
Frist five elements: [ 0.74968684 0.71457213 0.833188 0.48084545 0.86407363]
Last five elements: [ 0.65717077 0.76850474 0.29652017 0.856179 0.56074625]
------------
Weave Inline
Execution time: 2.48705 s
Frist five elements: [ 0.74968684 0.71457213 0.83318806 0.48084542 0.86407363]
Last five elements: [ 0.65717083 0.76850474 0.29652017 0.856179 0.56074625]
------------
pyCUDA
CUDA clock timing: 0.713028930664
Execution time: 2.04364 s
Frist five elements: [ 0.74968684 0.71457213 0.83318806 0.48084542 0.86407363]
Last five elements: [ 0.65717083 0.76850468 0.29652017 0.856179 0.56074625]
------------

Estoy un poco decepcionado con el rendimiento pyCUDA. Como soy nuevo en CUDA, probablemente haya algo que me falta aquí. Entonces, ¿dónde está el quid de la cuestión? ¿Estoy llegando a los límites del ancho de banda de la memoria global? Mala elección de bloques y gridsizes?

import numpy,time,math
import pycuda.autoinit
import pycuda.driver as drv
from pycuda.compiler import SourceModule
from scipy.spatial.distance import pdist
from scipy import weave
def weave_solution(x):
"""
OpenMP powered weave inline.
"""
N,DIM = numpy.shape(x)
L = ((N-1)**2+(N-1))/2
solution = numpy.zeros(L).astype(numpy.float32)
ncpu = 4
weave_omp = {''headers'' : [''<omp.h>''],
''extra_compile_args'': [''-fopenmp''],
''extra_link_args'' : [''-lgomp'']}
code = /
r''''''
omp_set_num_threads(ncpu);
#pragma omp parallel
{
int j,d,pos;
float r=0.0;
#pragma omp for
for (int i=0; i<(N-1); i++){
for (j=(i+1); j<N; j++){
r = 0.0;
for (d=0; d<DIM; d++){
r += (x[i*DIM+d]-x[j*DIM+d])*(x[i*DIM+d]-x[j*DIM+d]);
}
pos = (i*N+j)-(i*(i+1)/2)-i-1;
solution[pos] = sqrt(r);
}
}
}
''''''
weave.inline(code,[''x'',''N'',''DIM'',''solution'',''ncpu''],**weave_omp)
return numpy.array(solution)
def scipy_solution(x):
"""
SciPy High-level function
"""
return pdist(x).astype(numpy.float32)
def cuda_solution(x):
"""
pyCUDA
"""
N,DIM = numpy.shape(x)
N = numpy.int32(N)
DIM = numpy.int32(DIM)
L = ((N-1)**2+(N-1))/2
solution = numpy.zeros(L).astype(numpy.float32)
start = drv.Event()
end = drv.Event()
mod = SourceModule("""
__global__ void distance(float *x,int N,int DIM,float *solution){
const int i = blockDim.x * blockIdx.x + threadIdx.x;
int j,d,pos;
float r=0.0;
if ( i < (N-1) ){
for (j=(i+1); j<N; j++){
r = 0.0;
for (d=0; d<DIM; d++){
r += (x[i*DIM+d]-x[j*DIM+d])*(x[i*DIM+d]-x[j*DIM+d]);
}
pos = (i*N+j)-(i*(i+1)/2)-i-1;
solution[pos] = sqrt(r);
}
}
}
""")
func = mod.get_function("distance")
start.record()
func(drv.In(x),N,DIM,drv.Out(solution),block=(192,1,1),grid=(192,1))
end.record()
end.synchronize()
secs = start.time_till(end)*1e-3
print "CUDA clock timing: ",secs
return solution
if __name__ == ''__main__'':
# Set up data points
N = 25000
DIM = 3
x = numpy.random.rand(N,DIM).astype(numpy.float32)
print "-"*12
# Scipy solution
print "Scipy Pdist"
stime = time.time()
spsolution = scipy_solution(x)
stime = time.time()-stime
print "Execution time: {0:.5f} s".format(stime)
print "Frist five elements:", spsolution[:5]
print "Last five elements:", spsolution[-5:]
print "-"*12
# Weave solution
print "Weave Inline"
wtime = time.time()
wsolution = weave_solution(x)
wtime = time.time()-wtime
print "Execution time: {0:.5f} s".format(wtime)
print "Frist five elements:", wsolution[:5]
print "Last five elements:", wsolution[-5:]
print "-"*12
# pyCUDA solution
print "pyCUDA"
ctime = time.time()
csolution = cuda_solution(x)
ctime = time.time()-ctime
print "Execution time: {0:.5f} s".format(ctime)
print "Frist five elements:", csolution[:5]
print "Last five elements:", csolution[-5:]
print "-"*12

Editar:

He añadido la línea hash bang

#!/usr/bin/env python

en la parte superior del archivo y lo hizo ejecutable. Después de comentar el cálculo con weave.inline y scipy.spatial.distance.pdist , NVIDIA Visual Profiler solicita los siguientes resultados:

Answer 1

En este momento tiene 192 hilos cada uno actualizando N-1 posiciones, podría lanzar fácilmente más bloques / hilos.

Lo que quiere hacer es reemplazar este bucle for (j=(i+1); j<N; j++){ , reemplazarlo por N-1 hilos haciendo solo el bucle interno.

Si desea llevarlo más lejos, podría tener hilos N-1 * DIM haciendo cada uno la declaración en el ciclo interno, almacenar el resultado en la memoria compartida y finalmente hacer una reducción en eso. Ver Optimizar la Reducción Paralela en CUDA

Mirando esta línea:

r += (x[i*DIM+d]-x[j*DIM+d])*(x[i*DIM+d]-x[j*DIM+d]);

Las transacciones y el patrón de memoria no son uniformes y se fusionan. Tampoco sé si nvcc optimizará su expresión a solo dos transacciones de memoria en lugar de cuatro que se muestran aquí, ya que no sé si pycuda pasa -O3 a nvcc. Pon (x[i*DIM+d]-x[j*DIM+d]) en una variable de registro para asegurarte y simplemente cuadrarlo tú mismo.

De lo contrario, también puedes intentar #pragma unroll antes de cada ciclo for para, de ser posible, desenrollarlos.