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.