Un ejemplo básico, parte II: Implementación del Kernel

Como os prometí en el anterior post, hoy veremos la implementación del kernel para la reducción de vectores:

Nuestro kernel básicamente consiste en realizar la suma de una serie de elementos y guardarla en el primero de ellos. Cada ejecución se encargará de una posición de memoria diferente, con lo que no es necesario preocuparnos de implementar ningún mecanismo de sincronización. El número de elementos a sumar viene dado por el factor de desenrollado (UF). 

Los parámetros de entrada del Algoritmo son: Datos: puntero a la estructura que contiene los datos a sumar, Elementos: número de elementos que se suman, Impar: indica si el número de elementos es impar (para evitar la realización del cálculo del módulo, dentro del kernel, ya que es computacionalmente costoso), y finalmente UF, que se utilizará para indicar el número de elementos a sumar dentro de cada hilo.

Algoritmo 3: kernel de reducción de vectores
1	tid ← threadIdx.x + blockIdx.x * blockDim.x
2	t ← (ene/UF)
3	si tid < t hacer
4	para i←1 hasta UF hacer
5	a[tid] ← a[tid]  + a[(tid+(t*i))]
6	Si tid = 0 y impar > 0 hacer
7	para i←1 hasta impar+1 hacer
8		a[tid] ← a[tid]  + a[(ene-i)]

Algoritmo 1 Kernel de reducción de vectores

1.-La primera línea del algoritmo, “tid ← threadIdx.x + blockIdx.x * blockDim.x”, es la típica asignación para el cálculo de la posición global de ejecución del kernel. Así se identifica exactamente al hilo que se está referenciando en cada uno de los kernels y es la forma que tenemos para organizar el proceso en paralelo. Como en nuestro caso estamos trabajando únicamente con una dimensión, sólo se referencia al valor del campo “x” de cada variable (threadIdx, blockIdx y blockDim), pero podrían ser hasta tridimensionales, para ello se utilizan los campos “y” y “z” para referirnos a la segunda y tercera dimensión respectivamente.

2.- Dado que el número real de cálculos que tenemos que ejecutar no siempre coincide con el total de hilos ejecutados en los Kernels, debemos controlar esta situación para que el dispositivo no haga trabajo inútil (o redundante). Por este motivo la segunda línea del algoritmo calcula el número total de threads a emplear y en la tercera sólo se permite la ejecución del proceso de reducción a los hilos que están dentro del rango necesario.

4.-El siguiente paso define al bucle que empleamos para la suma de elementos del vector, que viene definido por el factor de desenrollado que empleamos en la llamada.

5.-La siguiente línea realiza la suma de elementos y la acumula en el primer elemento.

6.-Finalmente, tenemos una condición especial: sólo si se trata del primer hilo y el número de elementos a procesar es impar, adicionamos esos elementos impares también al primero.

Las operaciones realizadas en el siguiente ejemplo son: (3+6=9), (8+5=13), (4+2=6), y en la siguiente iteración (9+13+6=28).

En la siguiente Figura se puede observar gráficamente el proceso, para un factor UF de 2:

Figura 1. Reducción con UF = 2

Las propiedades fundamentales de este método de reducción que acabamos de proponer son:

·         No hay relación entre bloques ni hilos del mismo stream, con lo que no es necesario establecer ningún tipo de mecanismo para sincronizar el algoritmo en ningún punto. Es un algoritmo orientado al hilo y no al bloque.

·         Adecúa la petición de recursos del sistema en función al estado de su ejecución.

·         Sólo se realiza un desplazamiento de datos de memoria central a memoria de dispositivo por el total de elementos a sumar, y además solo se necesita desplazar el primer elemento a memoria principal para obtener el resultado de la reducción.

·         Con el control de UF se reduce el efecto de la latencia de memoria global del dispositivo, no siendo por tanto necesario realizar las lecturas alineadas con el orden de los hilos.

Ya hemos presentado un algoritmo para la reducción de vectores, que como sabréis, es una de las operaciones más típicas en multitud de escenarios, tanto en análisis de datos, como en optimización. El próximo post presentará los datos obtenidos en las pruebas. 

Un ejemplo básico, parte I: la llamada al kernel

Vamos a diseñar un algoritmo que permita realizar llamadas a nuestro kernel de tal forma que conociendo N número de elementos a procesar o tamaño del bucle en versión serie, utilice los mínimos recursos (de procesador gráfico) minimizando también las transferencias de memoria. Nos centraremos en un típico caso que en el que se emplea la programación paralela: Reducción de vectores. En la literatura hay multitud de ejemplos y diversos enfoques, nVidia también realiza varias implementaciones hasta 7 diferentes, van ilustrando los diferentes factores para mejorar el rendimiento. Yo implementé mi propio método, que ha resultado ser más rápido que las implementaciones de nVidia y de la librería Trust en varios escenarios. En el presente post y el siguiente lo presento:

Algoritmo para la llamada al kernel

Primero necesitamos obtener las características del dispositivo donde se van a ejecutar los programas, eso lo podemos obtener en tiempo de ejecución interrogando las capacidades de cada dispositivo. Las tres variables cuyo valor necesitamos saber son: N  número de elementos a procesar, UF factor de desenrollado (Unroll factor), que determinará la velocidad de segmentación de los datos y el número de instrucciones serie que ejecuta el kernel en cada hilo, y datos que es el puntero al bloque de datos a procesar.

Una vez determinados los valores de estas variables podremos calcular el factor de ocupación. El consumo máximo de tiempo se emplea en la migración de los datos desde memoria principal a la memoria del dispositivo, por lo que realizaremos sólo un proceso de carga y otro de descarga. El tamaño de bloques necesarios para la ejecución se irá reduciendo de forma logarítmica en cada iteración, ya que este algoritmo está orientado a las operaciones de reducción de vectores.

Para acumular los valores obtenidos se sobrescribe el vector de entrada con los resultados que se van obteniendo en cada cálculo, evitando así crear dos estructuras de datos, una de entrada y otra de salida. Además para cada llamada se aprovecha la información residente en la memoria del dispositivo, reduciendo la cantidad de datos transferida entre la memoria central y la del dispositivo.

Algoritmo 2: Algoritmo de llamada a un kernel
1	Obtener_caracteristicasGPU (max_bloques, max_hilos)
2	Determinar_factor_ocupacion(N,max_bloques, max_hilos, registros_kernel, nhilos)
3	UF ← 7
4	auxN  ← N
5	huerfanos  ←  (auxN modulo UF)
6	Datos_a_dispositivo(datos_dispositivo, datos)
7	Para i ← 1 hasta (log(N) / Log(UF))
8		Ejecutar en paralelo Reduce<<<(auxN+nhilos-1)/nhilos,nhilos> >>(datos_dispositivo, auxN, huerfanos, UF)
9		auxN ← auxN / UF
10		huerfanos  ←  (auxN modulo UF)
11	Datos_a_memoria_principal(datos, datos_dispositivo)

Algoritmo 1.Llamada al Kernel

Pasaremos a comentar cada una de las líneas de l algoritmo:

1.- Obtener_caracteristicas GPU (max_bloques, max_hilos): Indica los límites de la arquitectura seleccionada, que se utilizarán para aplicar las reglas que obtendrán el mejor factor de ocupación. Adicionalmente se obtienen otras cualidades de los dispositivos como la capacidad computacional (Compute capability) que indicará cuales son los dispositivos que permiten aprovechar las nuevas capacidades de CUDA 4, caracterizados por tener una capacidad igual o superior a 2.0.

2.- Si se puede estimar el tamaño N de nuestro proceso, y esto se puede traducir en el número de hilos necesarios para invocar al kernel, se podrá determinar el mejor factor de ocupación. Este cálculo óptimo de la relación hilos/bloques se puede realizar fácilmente porque al compilar el kernel se puede conocer el tamaño de de los registros necesarios para cada hilo,   y consecuentemente, se puede determinar el número máximo de hilos por bloque y de bloques por procesador.

3.- En la línea 3 hay que proporcionar un valor de Unroll Factor (UF). Para ello se deberán de realizar varias pruebas para obtener el mejor valor para este parámetro, ya que dependerá mucho de la estructura de datos en memoria que estemos utilizando, el tipo de memoria que utilicemos para su gestión y la coalescencia de los accesos que realicemos. En nuestras pruebas hemos utilizado la memoria global del dispositivo, y los accesos no han sido alineados con la ejecución de los kernels, Como veremos más adelante para nuestro algoritmo el mejor valor ha sido 7.

4.- En la línea 4 utilizamos una variable para guardar el tamaño de procesamiento remanente. El valor de esta variable se va reduciendo logarítmicamente conforme avanza el proceso.

5.- La variable "huérfanos" se utiliza para determinar cuándo se produce un elemento huérfano, el deberá de tratarse de forma especial.

6.- La línea 6 copia los datos desde memoria principal a la memoria global del dispositivo para ser procesados por el kernel.

7.- El bucle realiza una serie de llamadas al kernel en las que se ajustan los bloques que demandamos al dispositivo en función del remanente de ejecuciones. El número de iteraciones se corresponde con el cociente entre el logaritmo  de N y el logaritmo de UF-.

11.- Una vez concluidas las ejecuciones se descargará el dato deseado a memoria principal para obtener el resultado final.  La utilización del las llamadas descendentes, la vamos a entender a continuación en cuanto se exponga el Kernel.

Os dejo analizando el algoritmo, en el próximo post, veremos la implementación del Kernel. Hasta pronto

En busca del Kernel perdido

Con todo lo expuesto en los post anteriores debemos confeccionar nuestra lista de prioridades a la hora de diseñar los kernel que ejecuten las GPU. Resumiendo, debemos tener en cuenta los siguientes factores:

·         Factor de “unroll”

·         ILP: Nivel de Paralelización por instrucción

·         Factor de ocupación de la GPU, que se compone de:

o    Nº de bloques que se necesitan

o    Nº de hilos por cada bloque

o    Registros utilizados en el kernel

·         Transferencias de memoria CPU <-> GPU

·         Tipos de memoria a utilizar, Global, compartida,…

Teniendo todo esto en cuenta, nuestro kernel debe tener un número de instrucciones  mínimo que ocupe la latencia de la arquitectura, y se intentará invocar al kernel con el tamaño justo de bloques e hilos por bloques que maximicen la ocupación del tamaño  requerido. Se minimizará el número de transferencias de memoria entre CPU y GPU. Se intentará evitar las latencias de acceso a memoria global mediante accesos ordenados en memoria compartida.

En el próximo post, ilustraré con un ejemplo todo esto que hemos comentado

Los diversos usos de la memoria

Es muy importante saber que ventajas e inconvenientes tiene el uso de cada tipo de memoria de dispositivo, antes de realizar cualquier diseño, deberías tener claro cómo CUDA gestiona cada tipo de memoria, y cual es el coste o penalización que supone su utilización.

La siguiente Figura muestra los tipos de memoria en las GPUs y el ámbito de acceso de cada elemento

Figura 1. Los diferentes tipos de memoria en la GPU

Cada una de ellas tiene unos propósitos concretos que pasamos a detallar: Los registros y la memoria local son espacios de memoria accesibles únicamente por el hilo al que pertenecen. Por otro lado, a la memoria compartida  pueden acceder todos los hilos del mismo bloque, utilizando para ello mecanismos de sincronización que mantienen el orden en las lecturas y escrituras. Al resto de memorias  se puede acceder desde toda la malla que ejecuta un kernel. Se distinguen en que la memoria global no dispone de memoria caché, la memoria de constantes es de sólo lectura, y la memoria de texturas dispone de una caché espacial. 

El cuello de botella principal lo encontramos en la transferencia de memoria entre la CPU y la GPU, de modo que la penalización nos la encontramos en la latencia del acceso a la memoria global (entre 400 y 600 ciclos). Como solución posible se plantea el alineamiento de la memoria global para poder hacer lecturas ordenadas según la posición de cada hilo dentro de la tanda, se puede utilizar la memoria compartida para hacer la ordenación dentro de cada bloque. No siempre es posible hacer esta ordenación y el problema se complica cuando aumentan las dimensiones (matrices por ejemplo).

Entorno respetuoso de ejecución

A simple vista podríamos pensar que si un Grid solo puede ejecutar el mismo kernel y este es controlado solo por un dispositivo el cual solo es accesible desde el mismo thread del host, no es necesario tener en consideración otras interacciones entre distintos procesos. Hasta la versión 3.2 de CUDA esto era así,

Figura 1. CUDA 4.0 revoluciona la computación paralela en GPGPU

pero con la aparición de CUDA 4.0, se han eliminado ésta y muchas otras barreras: en cuanto a la plataforma GPUDirect™ se pueden realizar transferencias de memoria directamente de una GPU a otra sin pasar por memoria host; arquitectura 64-bit nativa; gestión de clúster; soporte multi-GPU en la que se pueden compartir múltiples GPUs a través de múltiples hilos (de host), un solo hilo (de host) puede acceder a todas las GPUs. En cuanto al modelo de programación: direccionamiento de memoria unificado, ahora se podrá contemplar la memoria de sistema y la memoria de cada una de las GPUs como una sola unidad de memoria simplificando así la migración de programas existentes. 

Además el nuevo toolkit viene acompañado de librerías OpenSource y plantillas de primitivas que ayudan y facilitan el manejo de vectores, matrices, etc. Un nuevo sistema de análisis de rendimiento permitirá obtener todas las métricas necesarias para afinar nuestros desarrollos.

Por eso es importante diseñar no sólo kernels, sino también las rutinas de host que los llaman, para que hagan un uso razonable de los recursos, por ejemplo modificando el tamaño de los bloques e hilos en las llamadas dinámicamente, en lugar de dejar un valor fijo para todas las llamadas.

No es “Speed Up” todo lo que reluce

Cuando realizamos optimizaciones en nuestro código, debemos ser conscientes de las limitaciones que impone la arquitectura, en muchos casos lo que puede parecer una mejora acaba siendo un cuello de botella.

 Aunque la técnica de “unroll” o desenrollado de bucles descrita en el post anterior parezca la “Panacea[1]”, hay que ser cautelosos, porque no siempre será así. Es necesario hacer un estudio detallado de cada caso para aplicar factores de desenrollado apropiados a cada uno de los algoritmos que se pretenden computar, ya que un desenrollado de bucles agresivo, puede aumentar demasiado las tareas de planificación (scheduling) y aumentar demasiado el uso de registro en el cuerpo del bucle, que en las CPU puede implicar una migración de los valores de los registros a memoria y esto puede desembocar en una degradación del rendimiento del programa. Por otro lado, si se utilizan factores de desenrollado muy altos, el tamaño del cuerpo del bucle resultante podría desbordar la caché de instrucciones, seguido de pérdidas de caché, y en consecuencia reduciéndose el rendimiento.

En programación GPGPU nos encontramos además con una serie de restricciones que limitan el rendimiento de los programas. Antes de proceder a la implementación deberemos estudiar el impacto que supone el desenrollado de bucles para los programas GPGPU en el contexto de las restricciones de esos recursos y hallar el factor de “unroll” óptimo para ese programa y dispositivo.

ILP y Factor de Ocupación (OF): como ya se ha visto, el nivel de paralelismo por instrucción juega un papel importante a la hora de aumentar el rendimiento, ya que el tamaño de la caché de instrucciones “I-cache” es limitado. Es importante no sobrecargar el cuerpo de la función con demasiadas instrucciones que provoquen un fallo de caché y la consiguiente penalización en el rendimiento.  Por otro lado en las GPU el factor de ocupación puede ser igual o más importante que el anterior a la hora de ocultar las latencias de acceso a memoria y de flujo de instrucciones. Debido a la interrelación que hay entre ILP y OF, se hace muy complicado estimar los valores óptimos. A continuación se muestra un caso concreto que sirve de ejemplo.

En el apartado B de esta sección se expuso brevemente el modelo de programación de CUDA compuesto por Mallas (Grids), Bloques (Blocks) e Hilos (Threads), ahora vamos a compararlo con la arquitectura de las GPUs para comprender mejor el problema del factor de ocupación. Una GPU está compuesta por x SM “stream multiprocessor”, cada uno dispone de c SP “stream processor”, en cada uno hay w Warps[2], que a su vez pueden ejecutar t Threads. La clave para comprender el problema del factor de ocupación reside en la diferencia real que existe entre la capacidad computacional de la arquitectura que utilicemos y la capacidad que solicitamos en la llamada al kernel. La Figura 1 muestra las capacidades para distintos dispositivos con arquitectura CUDA.

Device Name	GeForce GTX 480	Tesla C2050	GeForce GTX 295
Compute capability	2.0	2.0	1.3
Memory Informat.
Total global mem	1.536 MB	2.687 MB	896 MB
Total constant Mem	64 Kb	64 Kb	64 Kb
Max mem pitch	2048 MB	2048 MB	2048 MB
MP Information
Multiprocessor count	15	14	30
Shared mem per mp	48	48	16
Registers per mp	32768	32768	16384
Threads in warp	32	32	32
Warps/MP	48	48	32
Blocks/MP	8	8	8
Max threads/block	1024	1024	512
Max thread dim	(1024,1024, 64)	(1024,1024, 64)	(512, 512, 64)
Max grid dim	(65535, 65535, 65535)	(65535,  65535, 65535)	(65535, 65535, 1)
Warps/block in SP	6	6	4
threads in block	192	192	128
total pararell threads	2880	2688	3840

Figura 1. Capacidades de algunos dispositivos

De todos los valores de la tabla ahora mismo nos interesan: “Multiprocessor count”, que es el número de multiprocesadores que tiene el dispositivo, “Blocks/MP” bloques por multiprocesador, “”Warps/MP” que son agrupaciones de hilos dentro de cada bloque y “Threads in warp” que son el número de hilos que ejecutan en cada warp. Sin tener nada más en cuenta, podemos calcular el número total de hilos que cada dispositivo es capaz de ejecutar de forma concurrente: por ejemplo para el dispositivo GeForce GTX 480 tenemos 15 procesadores x 8 bloques x 6 warps x 48 hilos, en total 34.560 hilos concurrentes ejecutando el mismo kernel. Evidentemente la arquitectura nos está imponiendo además ciertas restricciones: “Max threads per block” está indicando un número (1024) mayor al teórico. Si en un bloque hay 6 warps y 32 hilos (192 hilos en total) ¿cómo es posible alcanzar esa cifra? CUDA puede utilizar hasta 8 bloques por procesador, es decir, que si utilizamos 192 hilos por bloque, dispondremos de 8 bloques por procesador, pero si necesitamos 256 hilos, dispondremos de 6 bloques, para 512 hilos solo 3 bloques, para 1024 solo un bloque por procesador. Cuando se solicitan los recursos en las llamadas a los kernels, CUDA realiza una segmentación y planificación de ejecución de todos los bloques solicitados, ahí es donde se pueden producir holguras si no hemos diseñado bien nuestra llamada, reduciéndose el factor de ocupación.

Otra cuestión a considerar, también muy importante es el uso de los registros que necesita nuestro Kernel, a mayor uso de registros, menor disponibilidad de bloques por procesador. Estas dos variables son las que van a modificar el factor de ocupación que, como hemos podido inferir, es el aprovechamiento físico de la arquitectura para cada tanda de instrucciones. El primero se maneja afinando los parámetros de llamada al núcleo y el segundo, mediante la creación de kernels lo más sencillos que se pueda (sin perder de vista el tiempo de latencia). Podemos saber el uso de registros que hace nuestro kernel añadiendo el parámetro (31) – Xptxas –v en la línea de compilación. 

nvcc -Xptxas –v ejemplo.cu ptxas info   : Compiling entry function 'acos_main' ptxas info   : Used 4 registers, 60+56 bytes lmem, 44+40 bytes smem,  20 bytes cmem[1],  12 bytes cmem[14]

Figura 2. Obteniendo información de nuestro kernel

---

[1] De panacea universal: Remedio que buscaban los antiguos alquimistas para curar todas las enfermedades

[2] El término “Warp” puede asimilarse al de “tanda” ya que hace referencia al conjunto de hilos de un procesador que se ejecutan en el mismo instante por pertenecer a la misma programación del scheduler. 

Unroll y Vencerás

Este, es el primero de una serie de 5 post en los cuales quiero presentar las claves básicas para desarrollar kernels CUDA de calidad.

Desde hace décadas, uno de los caballos de batalla en el mundo de los compiladores para CPU es el análisis de bucles y la reducción de iteraciones con el objetivo de optimizar el uso de recursos y mejorar los tiempos de ejecución. El resultado suele ser un bucle de menor número de repeticiones con instrucciones paralelas que son independientes unas de otras, en muchos casos esta simple acción puede reducir drásticamente el tiempo de procesamiento. El número de instrucciones serializadas, dentro del nuevo bucle coincide con el factor de desenrollado “unroll factor”. Este valor es clave a la hora de afinar y sintonizar los mejores valores para la paralelización como veremos más adelante. La Figura 1 muestra un ejemplo sencillo de esta técnica para UF = 4.

/* Antes del desenrollado */ for (i = 0; i < N; ++i) { c[i] = a[i] + b[i]; }

/* Despues del desenrollado */ for (i = 0; i < N - (4 - 1); i += 4) { c[i] = a[i] + b[i]; c[i+1] = a[i+1] + b[i+1]; c[i+2] = a[i+2] + b[i+2]; c[i+3] = a[i+3] + b[i+3]; } /* Bucle para el resto */ for (; i < N; ++i) { c[i] = a[i] + b[i]; }

                                                    Figura 1. Ejemplo de "unroll"

Los beneficios inmediatos de este proceso son: 

1. Se reduce el contador de instrucciones dinámicas, se reducen el número de comparaciones y operaciones ramificadas, para la misma carga de trabajo. 

2. Aumentan las posibles combinaciones para la planificación, debido a que surgen más instrucciones independientes, aumentando así el nivel de paralelización por instrucción (ILP). 

3. Aumenta la oportunidad para explotar la herencia de los registros y memoria local cuando los bucles externos se desenrollan y los internos se fusionan en el mismo juego de instrucciones iguales para todos.

En las arquitecturas GPGPU esta técnica también es habitualmente usada, aunque con otro enfoque, la diferencia de arquitecturas y modelos de programación necesita de una nueva forma de contemplar el impacto del desenrollado de bucles en el rendimiento en los programas para GPGPU.

A diferencia de los compiladores para CPU, no es tan fácil poder realizar optimizaciones automáticas de código usando esta técnica, por lo que el programador es el que deberá diseñar un buen código paralelizado en GPU. El modelo de programación CUDA y las restricciones de los dispositivos nos dan la guía para realizar la paralelización de bucles, como no es el objetivo de este post exponer los pormenores de la arquitectura y su forma de uso, sólo haremos una breve reseña con objeto de aportar claridad a esta exposición, para ampliar conceptos puede referirse el lector a las guías que acompañan el toolkit de CUDA.

 Básicamente podemos decir que un Kernel[1] equivale a una iteración del cuerpo de un bucle (en programación secuencial). La GPU tiene la habilidad de ejecutar de forma paralela el mismo kernel tantas veces como se defina en su llamada. Cada kernel corre en un hilo, y existen unas variables locales a cada uno de los kernels, pudiendo así cada kernel disponer de la información del índice o contador de iteración equivalente en caso de los bucles secuenciales. Los hilos (Thread), están organizados en bloques (Block), y los bloques en mallas (Grid). Los hilos pueden llegar a tener hasta un máximo de 3 dimensiones, y los bloques y las mallas hasta 2 o 3 dependiendo de la arquitectura. La siguiente línea de código muestra como un kernel puede obtener su identificador único (utilizando una dimensión en cada nivel):

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

Cómo pueden compartir y acceder a la memoria cada uno de los hilos y su ámbito se describe en el post "Los diversos usos de la memoria".


 El tamaño de la malla a utilizar y del bloque de hilos, es una elección que se deja al programador en función del problema a resolver. Esta decisión tendrá impacto en el rendimiento final del programa, que como veremos a continuación, está ligada al problema y sujeta a las restricciones que impone la arquitectura GPU y los recursos propios del dispositivo que se utilice.

---

[1] El Kernel es una porción de código que es invocada desde un proceso en CPU (host) y se ejecuta en un dispositivo GPU.