Una ordenación matricial en GPU

Hola a todos,

hoy quiero ilustraros un método de ordenación que implementé para la ordenación de vectores en GPU, la verdad es que primero lo ideé y como suele pasar en estas cosas luego busqué en internet y efectivamente el método ya existía y para mi sorpresa desde hace muchos años, en este caso aplicado a lás máquinas matriciales.

La solución de la que hablo está basada en el método de ordenación paralelo Parallel Neighbor-Sort propuesto por Nico Habermann(1), tambien conocido como odd-even sort, se diseño inicialmente para computación en arrays de procesadores a principios de los 70. Se trata realizar comparaciones de los elementos del vector dos a dos, siendo estos vecinos, realizando paralelamente Elementos/2 comparaciones en cada ejecución del kernel, las llamadas se intercalarán variando la paridad del elemento comparado, repitiendo este proceso hasta que no se produzca ningún cambio, esto determinará el final de la ordenación.

 

Es un método sencillo pero muy eficaz en sistemas multiprocesador con capacidades de paralelismo.

A continuación os ilustro el kernel para la ordenación en CUDA, está ordenando 3 vectoresunas reglas especificas demás complejo del que ya os hablaré, sirva ahora de ejemplo:

__global__ void ordenar_especies(int par,int *cambios,int *contar, int M, int D,  float *Puntos, int *Nivel, int *Orden, float *Valores) {

     __shared__ int cuenta;

     int tid = threadIdx.x + blockIdx.x * blockDim.x;

     int elem = ((3*M*3)+M)/2;

     float puntoA, puntoB, aux = 0.0f;

     int ordenA, ordenB = 0;

     int nivelA, nivelB = 0;

     int cambiar_puntos = 0;

     if (tid == 0){

          cambios = 0;

          if (contar[0] == 0){

              cuenta = 0;

          }

     }

     tid = (tid*2)+par;

     if (tid < elem) {

          puntoA =  tex1Dfetch(valoresTex,tid);

          puntoB =  tex1Dfetch(valoresTex,tid+1);

          ordenA =  tex1Dfetch(ordenTex,tid);

          ordenB =  tex1Dfetch(ordenTex,tid+1);

          nivelA =  tex1Dfetch(nivelTex,tid);

          nivelB =  tex1Dfetch(nivelTex,tid+1);

           if (ordenB > 0) {

               if (ordenA == 0) {

                  cambiar_puntos = 1;

               } else {

                  if (nivelB > nivelA){

                       cambiar_puntos = 1;

                  } else if ((nivelB == nivelA) && (puntoB > puntoA)){

                       cambiar_puntos = 1;

                  }

               }

           }

           if (cambiar_puntos == 1) {

                 for (int x=0;x<D;++x){

                     aux = Puntos[tid+(x*elem)];

                     Puntos[tid+(x*elem)] = Puntos[tid+1+(x*elem)];

                     Puntos[tid+1+(x*elem)] = aux;

                  }

                  Valores[tid] = puntoB;

                  Valores[tid+1] = puntoA;

                  Orden[tid] = ordenB;

                  Orden[tid+1] = ordenA;

                  Nivel[tid] = nivelB;

                  Nivel[tid+1] = nivelA;

                  cambios[0] = 1;

           }

           if (contar[0] == 0){

               if ((ordenA == 0) && (ordenB == 0)) {

                 __syncthreads();

                 cuenta = cuenta + 2;

               }

              else if ((ordenA == 0) || (ordenB == 0)) {

                 __syncthreads();

                 cuenta = cuenta + 1;

               }

           }

     }

    if (tid == (elem-1)){   //el ultimo que ponga la cuenta de su bloque

         if (contar[0] == 0){

            __syncthreads();

            contar[blockIdx.x] = elem-cuenta;

         }

     }

}

Hasta la próxima!!
____________________
(1)Parallel Neighbor-Sort (or the Glory of the Induction Principle). 5285 Port Royal Rd Sprigfield VA : National Technical Information Service, US Department of Commerce, 1972. AD0759248.

Bienvenido Mr. Kepler

En esta entrega compartiré con vosotros la nueva noticia: por fin tengo disponible en la plataforma del departamento un dispositivo con arquitectuta Kepler.

Según informa nVidia:

Kepler es la arquitectura de GPU más avanzada de nuestra historia. Optimizada para DirectX 11, Kepler no es solo rápida, también es supereficiente. 

cSegún parece se Consigue 3 veces más velocidad de procesamiento con NVIDIA® Kepler, aseguran que es la arquitectura de alta computación (HPC) más rápida y eficiente del mundo. Está dotada de tecnología y funciones de computación avanzadas, puede utilizarse en una extensa variedad de aplicaciones de cálculo científico y pone los sistemas de computación híbridos al alcance de mayor número de programadores e investigadores.

Sus principales novedades son:

SMXProporciona mayor velocidad de procesamiento y eficiencia gracias a su innovador multiprocesador de streaming, que permite dedicar más espacio a los núcleos de procesamiento que a la lógica de control.

GPU dinámica con Kepler
Simplifica la programación en la GPU ya que facilita la aceleración de bucles anidados paralelos, lo que significa que una GPU puede iniciar nuevos subprocesos de forma dinámica por sí misma, sin necesidad de volver a la CPU.

Hyper-Q de KeplerReduce el tiempo de inactividad de la CPU al permitir que múltiples núcleos de ésta utilicen una misma GPU Kepler, lo que mejora drásticamente la programabilidad y la eficiencia.

Estas tres mejoras en la arqutectura prometen ser un salto cualitativo y cuantotativo en la aplicación de GPGPU.

Yo por el momento cuento con una modesta GeForce GTX 660, las pruebas que he realizado en principio no me han dejado muy impresionado. Supongo que CUDA 5 permitirá exprimir a la perfección esta arquitectura.

Ya os contaré como se comporta esta nueva arquitectura con CUAD 5


Hasta entoces, feclices desarrollos!!!


 

Generación de números aleatorios en CUDA

Uno de los caballos de batalla más empleado en sistemas de optimización, minería de datos, sistemas expertos probabilísticos y simulación, es la utilización de números aleatorios, los números aleatorios nos permiten generar sobre la marcha poblaciones dentro de un espacio muestral controlado, que nos van a permitir, aplicando ciertas heurísticas, realizar inferencias y/o llegar a resultados. Normalmente estos métodos están clasificados dentro de lo que se conoce como "métodos Montecarlo" o basados en Montecarlo.

Por la intensa aplicación de estos métodos, se hace necesario disponer de mecanismos para la generación de números aleatorios que sea eficiente y con poco costo de procesamiento.

Existen varias soluciones para abordar el problema de la generación de números pseudoaleatorios con CUDA. Así, la librería CURAND provee al desarrollador de un importante surtido de funciones de alto nivel para la obtención de números aleatorios.

La librería CURAND contiene dos APIs: una para uso en código host, y la otra para uso en código de dispositivo.

Host API: para utilizar este API es necesario incluir curand.h en la lista de librerías en nuestro programa. Esta librería oculta al programador toda la problemática relacionada con el estado del generador de números aleatorios y su progresión de estado dentro del kernel de dispositivo. El funcionamiento básico es el siguiente:

1º) Se realiza la inicialización del estado del generador, para ello se utiliza la función curandCreateGenerator().

2º) Se especifican las opciones del generador, indicando además la semilla a utilizar. CURAND ofrece 2 tipos de generadores: CURAND_RNG_XORWOW, es un generador de números pseudoaleatorios que ha sido implementado usando el algoritmo XORWOW (32) introducido por George Marsaglia es una implementación de un generador del tipo xorshif también propuesto por él. Es un algoritmo muy rápido, calcula el siguiente elemento de la secuencia aplicando repetidamente un OR exclusivo sobre el resultado de aplicar un desplazamiento de bit al propio número. El código en lenguaje C publicado por Marsaglia es el siguiente:

unsigned long xor128(){

static unsigned long x=123456789,y=362436069,z=521288629,w=88675123;

unsigned long t;

t=(xˆ(x<<11));x=y;y=z;z=w;

return( w=(wˆ(w>>19))ˆ(tˆ(t>>8)) ) }

El Segundo tipo de generador que ofrece CURAND es un tipo  de generador números quasialeatorios: CURAND_RNG_SOBOL32, que permite generar secuencias de 32-bit y hasta 20.000 dimensiones. Es un algoritmo algo más complejo ideado por I.M. Sobol en 1967, que permite crear secuencias uniformes de baja discrepancia utilizando muchas dimensiones.

Una vez seleccionado el tipo de RNG[1] realizaremos la llamada a la función de la librería: curandSet-PseudoRandomGeneratorSeed() indicando la semilla y tipo de generador deseado.

3º) Reservar el espacio de memoria en el dispositivo.

4º) Invocar una función de generación de las disponibles pasando como parámetro el puntero a la memoria del dispositivo donde se generará la secuencia. Las distintas funciones de generación son:

- curandSetPseudoRandomGeneratorSeed(), permite generar números pseudoaleatorios y quasialeatorios cada elemento es de 32 bits, sin signo, donde todos sus bits son aleatorios.

-   curandGenerateUniform(), genera una secuencia uniforme-mente distribuida de números en coma flotante comprendidos entre (0.0, 1.0].

-  curandGenerateNormal(), permite obtener una secuencia de distribución normal, dados los valores de la media y la desviación estándar.

Para obtener secuencias con valores utilizando doble precisión se cuenta con las funciones homologas a las dos últimas: curandGenerateUniformDouble(), para la distribución uniforme  y curandGenerateNormaDouble() para la normal.

5º) Una vez terminados los procesos deseados, liberamos los recursos solicitados mediante el uso de la instrucción curandDestroyGenerator().

Device API: para utilizar esta modalidad, deberemos incluir la librería curand_kernel.h en nuestro programa. Como es evidente, las funciones aquí definidas se invocarán en los kernels de dispositivo o en funciones de dispositivo, diseñadas para ser ejecutadas por los kernel. En general, el uso de esta API permite una mejor optimización del código y ofrece mejores tiempos de respuesta. La librería permite, al igual que su versión para host, generar secuencias de números pseudoaleatorios y quasialeatorios.

Las funciones principales son las siguientes: curand_init(), para iniciar la secuencia, la función curand(), que permite obtener un número pseudoaleatorio entero sin signo. curand_uniform(), obtiene un número en coma flotante de una distribución uniforme en el intervalo (0.0, 1.0]. curand_normal(), obtiene un número en coma flotante de una distribución normal en el intervalo (0.0, 1.0] utilizando un valor de 0 para la media y de 1 para la desviación estándar.

También se dispone de versiones de estas dos últimas funciones para obtener secuencias con números de doble precisión: curand_uniform_double() y curand_normal_double(). Adicionalmente contamos con dos funciones que permiten obtener pares de elementos utilizando la estructura float2 y double2, que componen de un campo x y otro y, para obtener secuencias de distribución normal: curand_normal2() y curand_normal2_double() respectiva-mente.

El valor que se pasa a estas funciones es el estado del generador, que puede ser de dos tipos, dependiendo del generador que deseemos utilizar: curandState si deseamos utilizar un RNG pseudoaleatorio, o curandStateSobol32 si deseamos obtener valores de la secuencia quasialeatorios. Para el primer caso, en el que se trabaja con números pseudoaleatorios hay que pasar una semilla a la función, mientras que en el caso de quasialeatorios hay que pasar como argumento los vectores de dirección  y el desplazamiento. Finalmente cabe destacar una función para pasar por alto elementos de la secuencia, que es lo mismo que llamar varias veces a la función curand() pero mucho más rápido: skipahead(), en esta función indicamos el número de elementos a descartar y el estado.

Para la inicialización del generador de números aleatorios en CUDA, utilizamos la función curand_init(). Como particularidad hay que utilizar una estructura de datos que mantenga el estado del generador y lo actualice durante el proceso de generación de números.  Necesitamos mantener y pasar esa estructura en cada kernel que quiera generar números aleatorios:

curandState* devStates;

cudaMalloc ( &devStates, Elementos*sizeof( curandState ) );

/* inicializacion de numeros aletorios  */

randomize<<< ( Elementos + ( hilos – 1 )) / hilos, hilos >>> ( devStates, time(NULL) );

En el código anterior, en primer lugar se declara la estructura de estados de tipo curandState ya que se van a utilizar números pseudoaleatorios.

En la siguiente sentencia se reserva la memoria que necesita el dispositivo para actualizar los estados a cada hilo del kernel. Finalmente la última instrucción consiste en una llamada al kernel reandomize (la hemos llamado igual que la función existente en muchos lenguajes de alto nivel para la inicialización de la generación de números aleatorios). Como semilla se emplea el valor del tiempo actual del sistema, preparando así los estados para la obtención de números aleatorios.

__global__ void randomize( curandState *state, unsigned long semilla )

{

    int id = threadIdx.x + blockIdx.x * blockDim.x;

    curand_init ( semilla, id, 0, &state[id] );

}

El kernel randomize se encarga de iniciar el estado para cada hilo al llamar a curand_init. Separamos esta tarea en un kernel independiente con el objetivo de aumentar el rendimiento global del algoritmo. Por un lado para reducir el uso de la pila de memoria por hilo, ya que la generación de los estados de números aleatorios, en el peor de los casos puede llegar a consumir hasta 16KB. Por otro lado al separar esta tarea de la obtención de números aleatorios curand() permitimos al kernel, que está manejando la obtención, disponer de más velocidad y aumentar sus posibilidades de concurrencia al requerir menos recursos del procesador. Esto permite manejar más bloques por procesador.

Para obtener números aleatorios dentro del kernel utilizaremos la función curand_uniform( &localState ), que como hemos indicando anteriormente obtienen un elemento de la secuencia (en nuestro caso pseudoaleatoria) de elementos en coma flotante dentro del intervalo (0.0, 1.0].

---

[1] Las siglas en inglés de: Random Number Generator (Generador de Números Aleatorios)

CUDA 5.0 beta

 

Buenas Noticias, parece que el lanzamiento de las primeras GPU con arquitectura Kepler, estan acelerando la salida del nuevo SDK CUDA y sus drivers.

Ya está disponible la versión 5.0 Preview. Al parecer, vienen con nuevas herramientas y mejoras que harán más facil la implementación de aceleración de programas y procesos mediante GPU, y como no podía ser de otra manera, se incluyen los aceleradores basados en la arquitectura de cálculo NVIDIA Kepler™, la arquitectura de procesamiento más rápida, más eficiente y de mayor rendimiento que se haya inventado jamás.

Entre las principales características de CUDA 5 se incluyen (citando fuentes de nVidia):

Paralelismo dinámico: lleva la aceleración en la GPU a nuevos algoritmos
Los procesos de la GPU pueden generar dinámicamente nuevos procesos, lo que permite a la GPU adaptarse a los datos. Al reducir las interacciones con la CPU, el paralelismo dinámico simplifica enormemente la programación paralela. Además, posibilita la aceleración en la GPU de un número más amplio de algoritmos conocidos, tales como los utilizados en muchas aplicaciones de dinámica de fluidos computacional y métodos de refinamiento adaptable de mallas.

Llamadas a otras librerías desde la GPU: el origen de un ecosistema de desarrolladores
Una nueva librería BLAS de CUDA permite a los desarrolladores utilizar el paralelismo dinámico para hacer llamadas a sus propias librerías desde la GPU. El desarrollador puede diseñar API que permitan a otros desarrolladores ampliar la funcionalidad de sus kernels e implementar retrollamadas (callbacks) en la GPU para adaptar la funcionalidad de otras librerías externas que admitan llamadas desde la GPU. La posibilidad de “vincular objetos” proporciona una forma eficiente y conocida de desarrollar aplicaciones de GPU grandes ya que otorga a los desarrolladores la capacidad de compilar múltiples archivos CUDA de origen en diferentes archivos de objetos y vincularlos a librerías y aplicaciones de mayor tamaño.

Soporte de GPUDirect para RDMA: minimiza los cuellos de botella en la memoria del sistema
La tecnología GPUDirect proporciona comunicación directa entre las GPU y otros dispositivos PCI-E, y acceso directo a la memoria entre las tarjetas de red y la GPU. También reduce considerablemente la latencia de MPISendRecv entre los nodos de GPU de un cluster y mejora el rendimiento general de las aplicaciones.

NVIDIA Nsight Eclipse Edition: una forma rápida y fácil de generar código CUDA
Con NVIDIA Nsight Eclipse Edition, los programadores pueden programar, depurar y analizar el rendimiento de las aplicaciones de GPU dentro del conocido entorno de desarrollo integrado basado en Eclipse para las plataformas Linux y Mac OS X.Un editor de CUDA integrado y los ejemplos de código CUDA ayudan a agilizar la programación, y la refactorización automática del código facilita el traslado de bucles basados en la CPU a kernels CUDA.Un sistema de análisis integrado proporciona análisis automatizado del rendimiento y una guía de instrucciones detalladas para eliminar los cuellos de botella del código. Por último, la diferenciación visual de la sintaxis ayuda a distinguir mejor entre el código de la GPU y el de la CPU.

Estoy deseando probarlo!!  Un saludo a todos.

Paralelizar en CPU vs GPU

Hoy vamos a tratar dos cuestiones clave a la hora de abordar una implementación en GPU:

los que estén acostumbrados a implementar algoritmos paralelos en CPU tendrán un acceso algo mas llano a la programación GPGPU, pero hay ciertas diferencias a tener en cuenta: la paralelización en CPU difiere considerablemente de la paralelización sobre GPU. La diferencia fundamental se encuentra en el número de procesadores con el que contamos en una y otra arquitectura, si bien en CPU contamos con pocos procesadores pero muy rápidos, en GPU disponemos de muchos procesadores algo más lentos.

La consecuencia directa de esta diferencia implica adoptar un enfoque diferente, orientado a la unidad mínima de procesamiento independiente, evitando en la medida de lo posible, la dependencia entre estas unidades, esto lo tenemos que reflejar tanto en el diseño de las estructuras de datos, como en las tareas que realicemos sobre estas.

Además en GPU podemos adoptar 2 enfoques diferentes:

A.     Enfoque orientado al bloque

Cuando implementamos un kernel en CUDA, éste puede diseñarse de tal manera que deba existir cierta comunicación entre los hilos que ejecutan las tareas. De alguna manera estos hilos colaboran entre sí y pueden pasarse información entre ellos, estableciendo mecanismos de sincronización para acceder a la memoria local donde residen las variables que son comunes al proceso. En este enfoque existe una limitación, la comunicación o colaboración sólo es posible entre hilos que pertenecen al mismo bloque de ejecución. La ventaja principal es que es posible realizar una comunicación muy rápida mediante la memoria local, ya que su tiempo de latencia es mínimo (en comparación con el acceso a memoria global).

La implementación de este enfoque requiere entonces de dos fases, una primera en la que todos los hilos de cada bloque procesan en colaboración obteniendo un resultado por bloque, y una segunda fase en la que son procesados los resultados de cada bloque. Esta tarea se puede ejecutar en otro kernel de GPU o bien directamente en CPU.

B.     Enfoque orientado al hilo

Este enfoque implica una estructura más plana, en la que para una determinada ejecución existe independencia total entre los elementos de cómputo y la tarea que se aplica. Por ejemplo, una tarea que asigne un número aleatorio a cada elemento de un vector puede ser un caso típico.

Aunque inicialmente puede parecer que sí tienen dependencia dada la naturaleza pseudo-aleatoria de estos números, es posible hacer una descomposición de tareas de forma que se pueda dar la independencia de procesos.

Nuestro algoritmo de reducción de vectores propuesto en post anteriores, está implementado utilizando este enfoque. Así, para desligar la dependencia se han dividido las ejecuciones del kernel en niveles, que son ejecutados secuencialmente pero para cada ejecución todas las operaciones son independientes.

El enfoque orientado al hilo permite de una manera más sencilla alinear las lecturas de memoria con cada uno de los hilos que hacen uso de esa memoria, reduciendo el “lag” que produce la latencia de la memoria sobre todo cuando esa relación es de uno a uno.

En resumen: dependiendo del tipo de problema una estrategia será más beneficiosa que la otra, aunque en muchos casos lo idóneo será una mezcla de todas.

Un ejemplo básico, parte III: Análisis del rendimiento

Para realizar las pruebas se ha utilizado la plataforma "Davinci" que es una computadora gestionada por el Grupo de Investigación TIC-146 Supercomputación: Algoritmos de la Universidad de Almería. Este servidor multicore dedicado a computación GPGPU cuenta con:

 2 Intel Xeon Quad Core E5640 2.67 GHz 64 bits

 GeForce GTX 295 (2x240 cores)

 2 Tesla C2050 (Fermi) (448 cores)

 GeForce GTX 480 (480 cores)

Todas las GPU son Compatibles con CUDA, si bien sólo las Tesla C2050 y la GeForce GTX 480 ofrecen las prestaciones de cómputo 2.0 necesarias para aprovechar al máximo las prestaciones de CUDA 4.

El sistema operativo de este servidor es un Linux Version  de kernel 2.6.26-2-amd64, con una distribución Debian GNU/Linux 5.0.8 (lenny). La versión de Cuda instalada es la 4.0.

Los programas se han desarrollado en lenguaje C y CUDA, utilizando los compiladores gcc (The GNU Compiler Collection) y nvgc (el compilador CUDA de nVidia).

Las pruebas se han realizado sobre diferentes tamaños de vectores partiendo de 10.000 elementos, y con incrementos de ese mismo número hasta llegar a los 1.000.000 de elementos del vector.

A.     Reducción de matrices

En las pruebas no se ha tenido en cuenta el tiempo de introducir valores aleatorios a los elementos del vector, ya que es común para los tres escenarios. Los tres escenarios comparados son los siguientes: por un lado la versión secuencial de la suma ejecutada en la CPU, nuestro algoritmo de reducción de vectores propuesto, y la reducción de vectores que la nueva librería Thrust® ofrece y que simplifica la programación paralela para CUDA.

Cada uno de los test se ha realizado 10 veces consecutivas y se ha calculado la media aritmética de los resultados obtenidos.

Para medir los tiempos de respuesta de cada uno de los algoritmos hemos utilizado la gestión de eventos de CUDA. Dentro del tiempo de cómputo para las versiones paralelas, se ha incluido el tiempo de copia de memoria principal a memoria de dispositivo y viceversa. La siguiente tabla muestra los resultados obtenidos en el test:

N	CPU	GPU1	Thrust
10000	0.0000382	0.00007563	0.00031225
20000	0.0000758	0.00010429	0.00031417
30000	0.0001138	0.00013749	0.00033243
40000	0.0001512	0.00016749	0.00035026
50000	0.0001884	0.00019928	0.00036693
60000	0.0002266	0.00022946	0.0003838
70000	0.0002640	0.00025606	0.00040174
80000	0.0003020	0.00028915	0.00041772
90000	0.0003392	0.00031836	0.00043793
100000	0.0003770	0.00035073	0.00046362
110000	0.0004146	0.00038296	0.00049638
120000	0.0004536	0.00041684	0.00051473
130000	0.0004896	0.00044498	0.00056173
140000	0.0005274	0.00046391	0.00058575
150000	0.0005678	0.00048342	0.0005885
160000	0.0006090	0.00050197	0.00060389
170000	0.0006442	0.00051199	0.00061768
180000	0.0006782	0.00052612	0.00063185
190000	0.0007172	0.00054233	0.00064476
200000	0.0007548	0.00055549	0.00066341

Tabla 1. Tiempo (en segundos) de respuesta para los algoritmos

La tabla anterior muestra en la columna N el número de elementos del vector tratado, la columna CPU muestra los tiempos de respuesta obtenidos para el algoritmo en serie sobre la CPU, la columna GPU1 se refiere al algoritmo de reducción propuesto y finalmente Thrust muestra los tiempos de respuesta obtenidos utilizando las librerías Thrust®.

Como se puede observar el rendimiento obtenido con nuestra propuesta de algoritmo es bastante superior al rendimiento obtenido por el algoritmo de reducción de la librería Thrust de nVidia.

Figura 1. Comparativa rendimientos obtenidos

Se realizó una segunda comparativa entre CPU y GPU1 con un crecimiento exponencial de elementos de la forma 1024ⁱ para i = 1 hasta 12 con los siguientes resultados:

N	CPU	GPU1	SpeedUp
1024	0.0000042	0.00004873	0.08618921
2048	0.0000082	0.00004486	0.18279091
4096	0.000016	0.00005397	0.296461
8192	0.0000318	0.00006936	0.45847751
16384	0.0000624	0.00009197	0.67848211
32768	0.000124	0.00014606	0.84896618
65536	0.0002474	0.00024655	1.00344758
131072	0.000494	0.00045282	1.09094121
262144	0.000992	0.00070173	1.41364912
524288	0.0019882	0.00116573	1.70554073
1048576	0.0039892	0.00213745	1.86633605
2097152	0.0109062	0.00441923	2.46789599

Tabla 2. SpeedUp respecto a GPU1

En la siguiente figura se puede observar cómo la pendiente de la tendencia creciente exponencial de la GPU1 es minorada por el factor de SpeedUp obtenido en cada tramo:

Figura 2. Evolución del SpeedUp al aumentar N

B.     Determinación del Punto Muerto Computacional[1]

Trataremos de determinar el punto donde es indiferente realizar los cálculos en CPU respecto de la GPU ó

GPU respecto de la CPU porque el rendimiento sea el mismo. En las pruebas anteriores, hemos podido observar cómo en los tramos donde el número de elementos que componen el vector a reducir es bajo, se daba la situación que el algoritmo en serie ejecutado en la CPU obtenía mejores resultados que sus homólogos paralelizados sobre GPU. Evidentemente debemos de tener esto en cuenta a la hora de decidir en donde debemos ejecutar cada trabajo.

Dicho punto coincidirá con el punto de corte de las funciones que representan el tiempo de ejecución que emplea cada algoritmo en procesar exactamente el mismo trabajo (incluidos en este caso para los algoritmos de GPU, los tiempos de proceso consumidos en la transferencia de memoria principal a memoria del dispositivo y una vez realizado el cálculo, vuelta a recibir los datos a memoria principal) y obtener exactamente el mismo resultado.

Si realizamos un zoom sobre los datos tratados apreciamos

Figura 3. Punto de corte función CPU y GPU1

un punto de corte entre las funciones de tiempo de respuesta sobre el tamaño N del vector. Así, en la Figura 18 podemos observar cómo ese punto de corte está cerca de los 70000 elementos. Podemos decir que en este punto es indiferente la utilización de un algoritmo u otro en términos de tiempos de ejecución. Por debajo de este valor sería mejor computar sobre la CPU, mientras que por encima de ese valor sería mejor hacerlo en GPU.

Evidentemente, esta visión es algo ingenua, y se deberían tener otros factores en cuenta a la hora de tomar esta decisión.

- Coste de utilización: es el coste que implica la utilización de cada escenario. En muchos casos (como en grandes centros de computación) se calcula como la suma del coste de adquisición  más el coste en potencia (watios) por euro de cada recurso.

- Coste de oportunidad[2]: es un concepto más sutil pero muy interesante a la hora de tomar decisiones económicas. Aplicado a nuestro caso, se puede definir como el coste que supone la decisión de computar un algoritmo en un recurso (que es escaso) y no computar o descartar la computación de otro algoritmo en ese mismo recurso.

Teniendo en cuenta los costes anteriores se podrán tomar decisiones por ejemplo: de utilizar una GPU para realizar un determinado proceso aunque el tiempo de respuesta sea similar o peor al obtenido en CPU, debido a que el coste por watio de la GPU es muy inferior, y a que el coste del trabajo que puede realizarse en la CPU sea mucho mayor.

C.      Determinación del Factor de Desenrollado

Ya se habló en post anteriores de la importancia de este factor, y de cómo ha sido un caballo de batalla en los compiladores desde los inicios de la informática.

El objetivo fundamental que perseguimos con esta técnica es compensar los efectos negativos en el rendimiento que producen los tiempos de latencia que se producen en los accesos a memoria.

De forma intuitiva, podemos expresar este efecto como el que se produce en las líneas de autobuses urbanos: un autobús puede transportar 50 pasajeros a la vez, tarda prácticamente lo mismo si trasporta 1 o si trasporta a 50. El problema está en que la compañía de autobuses espera optimizar su beneficio, y establece unos horarios de salida y llegada, por eso, aunque un pasajero suba al autobús, deberá esperar cierto tiempo hasta ser trasportado a su destino. En nuestro caso los pasajeros son los datos que necesitamos tratar, la parada de origen es la memoria global del dispositivo, y la parada final, los registros locales del hilo que ejecuta un kernel. La arquitectura CUDA, está diseñada para esperar en la parada de origen entre 400 y 600 ciclos de reloj, además hay que tener en cuenta el número de pasajeros que se permite trasportar en cada viaje.

Como podemos deducir, para optimizar el proceso, cada vez que realicemos un acceso a memoria global del dispositivo, lo haremos de tal manera que llenemos el autobús, y todos esos datos trasportados, los procesaremos en el mismo hilo.

Para nuestro algoritmo, hemos realizado una serie de pruebas partiendo de un factor de desenrollado de 2 hasta 10, 25, 50 y 100, la siguiente tabla muestra los tiempos de respuesta obtenidos en cada caso.

UF	GPU1
2	0.00010693
3	0.00008918
4	0.00008308
5	0.00007948
6	0.00008024
7	0.00007638
8	0.00007826
9	0.00008013
10	0.00007762
25	0.0000807
50	0.00008548
100	0.00011343

Tabla 3. Tiempos de respuesta para distintos valores de UF

Hemos determinado de forma empírica que el mejor factor es el de 7. Para otros algoritmos puede variar ya que todo esto depende del tipo de dato que se esté transportando, del tamaño de los registros necesarios en el kernel, etc.

Figura 4. Selección del mejor factor de Unroll

---

[1] Este término ha sido acuñado para el presente trabajo como analogía al término económico “Punto Muerto” o “Umbral de rentabilidad” donde los ingresos son iguales a los gastos 

[2] Es el término económico utilizado para referirnos al coste de la mejor opción no realizada. Este término fue inventado por Friedrich von Wieser en su Theorie der gesellschaftlichen Wirtschaft (Teoría de la economía social). El coste de la oportunidad es aquello a lo que renunciamos cuando tomamos una decisión económica.