viernes, 7 de febrero de 2014

Funciones __device__


En ejercicios anteriores hemos aprendido como usar la palabra reservada __global__ para marcar una función como código que el host puede llamar y que origina una invocación de un kernel paralelo en la GPGPU.

Con una función __global__, cada thread CUDA sigue su propio patrón de ejecución en forma serial. En CUDA, los kernels consisten en mayormente en código C/C++ que pueden ser muy rápidos.

Como programadores en paralelo, tendremos la necesidad de abstraer y encapsular el kernel en funciones. La palabra reservada __device__ nos permite marcar una función que es llamada desde threads ejecutándose en la GPGPU. Por ejemplo:

__device__ float mi_funcion_device(float x)
{
   return x + 1;
}

Aunque las funciones __device__ son similares a las funciones __global__ en que pueden ser ejecutadas por threads CUDA, de hecho se comportan más como funciones normales en C. A diferencia de funciones __global__, las funciones __device__ no pueden ser configuradas con <<B, T>>, y no están sujetas a ninguna restricción especial en los tipos de sus parámetros o sus resultados. El código del host no puede llamar a las funciones __device__ directamente; si se desea acceder a una función __device__, necesitamos escribir una función __global__ que la llame.

Como se podría esperar, las funciones __device__ pueden llamar a otras funciones __device__:

__device__ float mi_segunda_funcion_device(float y)
{
   return mi_funcion_device(y)/2;
}

Siempre y cuando no se llamen a si mismas:

__device__ int float mi_funcion_device_ilegal_recursiva(int x)
{
   if (x == 0) return 1;
   return x * mi_funcion_device_ilegal_recursiva(x-1);
}

El siguiente código muestra como se podría usar la funciones __device__ para empaquetar varios bits de código cuando se desarrolla un kernel CUDA.

#include <stdlib.h>
#include <stdio.h>

__device__ int get_global_index(void)
{
  return blockIdx.x * blockDim.x + threadIdx.x;
}

__device__ int get_constant(void)
{
  return 7;
}

__global__ void kernel1(int *array)
{
  int index = get_global_index();
  array[index] = get_constant();
}

__global__ void kernel2(int *array)
{
  int index = get_global_index();
  array[index] = get_global_index();
}

int main(void)
{
  int num_elements = 256;
  int num_bytes = num_elements * sizeof(int);

  int *device_array = 0;
  int *host_array = 0;

  // reserva memoria
  host_array = (int*)malloc(num_bytes);
  cudaMalloc((void**)&device_array, num_bytes);

  int block_size = 128;
  int grid_size = num_elements / block_size;

  // lanza kernel1 e inspecciona sus resultados
  kernel1<<<grid_size,block_size>>>(device_array);
  cudaMemcpy(host_array, device_array, num_bytes, cudaMemcpyDeviceToHost);

  printf("resultados de kernel1:\n");
  for(int i = 0; i < num_elements; ++i)
  {
    printf("%d ", host_array[i]);
  }
  printf("\n\n");

  // lanza kernel2 e inspecciona sus resultados
  kernel2<<<grid_size,block_size>>>(device_array);
  cudaMemcpy(host_array, device_array, num_bytes, cudaMemcpyDeviceToHost);

  printf("resultados de kernel2:\n");
  for(int i = 0; i < num_elements; ++i)
  {
    printf("%d ", host_array[i]);
  }
  printf("\n\n");

  // liberar memoria
  free(host_array);
  cudaFree(device_array);
  return 0;
}

jueves, 6 de febrero de 2014

Jugando con dimensiones en grids y bloques

En CUDA es posible tener estructura de 1D, 2D y 3D para bloques y grids de threads, lo cual se explicó en Identificación de un thread. Para aclarar la forma de utilizar estas posibles configuraciones se resume en lo siguiente, usando configuraciones 1D y 2D para bloques y 1D, 2D y 3D para threads:


  1. Array 1D de bloques donde cada bloque tiene un array 1D de threads.

  2. UniqueBlockIndex = blockIdx.x;
    UniqueThreadIndex = blockIdx.x * blockDim.x + threadIdx.x;
       
  3. Array 1D de bloques donde cada bloque tiene un array 2D de threads.

    UniqueBlockIndex = blockIdx.x;
    UniqueThreadIndex = blockIdx.x * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
                          
  4. Array 1D de bloques donde cada bloque tiene un array 3D de threads.

    UniqueBlockIndex = blockIdx.x;
    UniqueThreadIndex = blockIdx.x * blockDim.x * blockDim.y + threadIdx.y * blockDim.x + threadIdx.x;
                          
  5. Array 2D de bloques donde cada bloque tiene un array 1D de threads.

    UniqueBlockIndex = blockIdx.y * gridDim.x + blockIdx.x;
    UniqueThreadIndex = UniqueBlockIndex * blockDim.x + threadIdx.x;
                                                             
  6. Array 2D de bloques donde cada bloque tiene un array 2D de threads.

    UniqueBlockIndex = blockIdx.y * gridDim.x + blockIdx.x;
    UniqueThreadIndex = UniqueBlockIndex * blockDim.y * blockDim.x + threadIdx.y * blockDim.x + threadIdx.x;
                          
  7. Array 2D de bloques donde cada bloque tiene un array 3D de threads.

    UniqueBlockIndex = blockIdx.y * gridDim.z * + blockIdx.x;
    UniqueThreadIndex = UniqueBlockIndex * blockDim.z * blockDim.y * blockDim.x + threadIdx.z * blockDim.y * blockDim.z + threadIdz.y * blockDim.x + threadIdx.x;

viernes, 31 de enero de 2014

Grids Bidimensionales

Esta entrada mostrará un código donde el kernel será lanzado con bloques de 2D, y un grid de threads 2D para realizar operaciones sobre una matriz de datos. Las operaciones ejemplifican la conversión de los índices del thread y del bloque en 1D, esto debido a que el manejo de los datos en la GPGPU es 1D.


#include <stdlib.h>
#include <stdio.h>

__global__ void kernel(int *array) {

// se obtienen los índices x y y para cada thread.
int index_x = blockIdx.x * blockDim.x + threadIdx.x;
int index_y = blockIdx.y * blockDim.y + threadIdx.y;

// el array viene en forma 1D (a pesar de ser matriz), por lo
// que hay que covertir de 2D a un índice 1D

int grid_width = gridDim.x * blockDim.x;
int index = index_y * grid_width + index_x;

// convierte el índice de bloque de 2D en un índice 1D.
int result = blockIdx.y * gridDim.x + blockIdx.x;

// escribe el resultado
array[index] = result;
}
int main(void) {
int num_elements_x = 16;
int num_elements_y = 16;

int num_bytes = num_elements_x * num_elements_y * sizeof(int);

int *device_array = 0;
int *host_array = 0;

// reserva memoria en host
host_array = (int*) malloc(num_bytes);
// reserva memoria en GPU
cudaMalloc((void**) &device_array, num_bytes);

// crea bloques bidimensionales de threads de 4 x 4
dim3 block_size;
block_size.x = 4;
block_size.y = 4;

// configura un grid bidimensional de 4 x 4
dim3 grid_size;
grid_size.x = num_elements_x / block_size.x;
grid_size.y = num_elements_y / block_size.y;

// se envían el arreglo de bloques y el grid de threads por bloque
kernel<<<grid_size, block_size>>>(device_array);

// copia al host
cudaMemcpy(host_array, device_array, num_bytes, cudaMemcpyDeviceToHost);

// impresión
for (int row = 0; row < num_elements_y; ++row) {
for (int col = 0; col < num_elements_x; ++col) {
printf("%2d ", host_array[row * num_elements_x + col]);
}
printf("\n");
}
printf("\n");

// liberar memoria
free(host_array);
cudaFree(device_array);
}


Salida del programa:


 0  0  0  0  1  1  1  1  2  2  2  2  3  3  3  3 
 0  0  0  0  1  1  1  1  2  2  2  2  3  3  3  3 
 0  0  0  0  1  1  1  1  2  2  2  2  3  3  3  3 
 0  0  0  0  1  1  1  1  2  2  2  2  3  3  3  3 
 4  4  4  4  5  5  5  5  6  6  6  6  7  7  7  7 
 4  4  4  4  5  5  5  5  6  6  6  6  7  7  7  7 
 4  4  4  4  5  5  5  5  6  6  6  6  7  7  7  7 
 4  4  4  4  5  5  5  5  6  6  6  6  7  7  7  7 
 8  8  8  8  9  9  9  9 10 10 10 10 11 11 11 11 
 8  8  8  8  9  9  9  9 10 10 10 10 11 11 11 11 
 8  8  8  8  9  9  9  9 10 10 10 10 11 11 11 11 
 8  8  8  8  9  9  9  9 10 10 10 10 11 11 11 11 
12 12 12 12 13 13 13 13 14 14 14 14 15 15 15 15 
12 12 12 12 13 13 13 13 14 14 14 14 15 15 15 15 
12 12 12 12 13 13 13 13 14 14 14 14 15 15 15 15 
12 12 12 12 13 13 13 13 14 14 14 14 15 15 15 15 

miércoles, 22 de enero de 2014

Suma de vectores II

El principal propósito de trabajar con una GPGPU es el de potenciar el cálculo con cantidades inmensa de datos. Hemos descrito anteriormente un post donde se realizaba la suma de vectores funcionando perfectamente, sin embargo a continuación se muestra un código que funciona para una cantidad de elementos mucho mayor del vector de números:


__global__ void sumaVector(long *v1, long *v2, long *v3, long N) {

int threadId = blockIdx.x * blockDim.x + threadIdx.x;
while (threadId < N) {
v3[threadId] = v1[threadId] + v2[threadId];
threadId += blockDim.x * gridDim.x;
}
}

int main(int argc, char** argv) {
        long N = 9000000; // 9 millones

long *v1, *v2, *v3;
v1 = (long *) malloc(N * sizeof(long));
v2 = (long *) malloc(N * sizeof(long));
v3 = (long *) malloc(N * sizeof(long));

for (long i = 0; i < N; i++) {
// datos de prueba
v1[i] = 10;
v2[i] = 11;

}

printf("%ld", v1[1000001]); //OK

long *dv1, *dv2, *dv3;

cudaMalloc((void**) &dv1, N * sizeof(long));
cudaMalloc((void**) &dv2, N * sizeof(long));
cudaMalloc((void**) &dv3, N * sizeof(long));

// copiando memoria a la GPGPU
cudaMemcpy(dv1, v1, N * sizeof(long), cudaMemcpyHostToDevice);
cudaMemcpy(dv2, v2, N * sizeof(long), cudaMemcpyHostToDevice);

// número de bloques
int B = 1024;
        int T = 1024;

// Llamando a ejecutar el kernel
sumaVector<<<B, T>>>(dv1, dv2, dv3, N);

// copiando el resultado a la memoria Host
cudaMemcpy(v3, dv3, N * sizeof(long), cudaMemcpyDeviceToHost);

cudaFree(dv1);
cudaFree(dv2);
cudaFree(dv3);

printf("%ld", v3[900001]); //OK

        return (EXIT_SUCCESS);
}

Características importantes a tomar en cuenta

  • Este cálculo se realiza en dimensión 1D.
  • Se define la llamada del kernel para 1024 bloques y cada bloque con 1024 threads. Por lo tanto la GPGPU ejecutará 1,048,576 threads. 
  • El vector tiene un tamaño de 9 millones de elementos (se usó valores long para disponer en cualquier caso de números grandes).
  • Debido a que el número de elementos excede el número de threads en el kernel, la codificación de éste incluye un ciclo que indica que se deberá estar realizando hasta que se terminen de calcular las 9 millones de veces.

martes, 1 de octubre de 2013

Generar números aleatorios en CUDA


La generación de números aleatorios (RNG) tiene diversas aplicaciones en simulaciones computacionales, algoritmos evolutivos, método de Monte-Carlo, entre otros y por lo tanto, serán de importancia para el cálculo en GPGPUs.

En estos problemas, podemos distinguir:
  1. Los Números Aleatorios "Verdaderos" (True Random Number): Los más complicados de generar, se basan en métodos no determinísticos, generalmente en fenómenos físicos (por ejemplo, radioactivos, atmósfera) que se espera tengan resultados aleatorios. La generación de éstos números no es periódica, es decir, que no se repite la secuencia de números generados. El sitio Random.org proporciona servicios gratuitos de generación de números aleatorios de éste tipo.
  2. Los Pseudo Números Aleatorios: Usan algoritmos computacionales capaces de producir secuencias largas de números aparentemente aleatorios, los cuales están determinados por un valor inicial al que se le denomina semilla (seed). En estos números la secuencia eventualmente se repite.
Obviamente es imposible generar números aleatorios verdaderos en una computadora determinística. La función de aleatoriedad aplica algún tipo de transformación sobre otro número, determinando una sucesión que "parece" aleatoria. De cualquier forma, si dos generaciones de números empezaran a partir de la misma semilla, el resultado sería el mismo.

En la programación en C/C++, se hace gran uso de la función rand() para generar estos tipos de números pseudoaleatorios. Las semilla generalmente más usada es el tiempo, lo cual genera resultados aceptables.

Sin embargo,al generar números aleatorios en una GPGPU puede volverse algo complicado. La solución más sencilla e ingenua es crear todos los números aleatorios necesarios en el host y colocarlos en la memoria global de la GPGPU (pre-generación). La desventaja está en el bandwith necesario para transferir dichos números a la memoria del dispositivo.

Por lo tanto es más eficiente generar los números aleatorios directamente en la memoria del dispositivo en un kernel exclusivamente dedicado a ello. Para ello se puede aprovechar la paralelización en la generación de dichos números, obviamente tomando en cuenta que es necesario partir de diferentes semillas, si no, el número generado sería el mismo. La biblioteca NVIDIA CURAND hace más fácil la creación de números dentro del kernel del dispositivo. Dichos números estarán almacenados en la memoria local, y disponibles para el cálculo que se requiera.

Los números generados por CURAND son pseudoaleatorios y/o cuasialeatorios. Una secuencia de pseudoaleatorios satisface la mayoría de las propiedades estadísticas de una secuencia de números verdaderamente aleatorios, sin embargo, es generada por un algoritmo determinista. Una secuencia cuasialeatoria de puntos n-dimensionales es determinada por un algoritmo determinista diseñado para llenar el espacio n-dimensional.

A continuación se muestra un código de ejemplo, genera un vector de números flotantes en el dispositivo. Para efectos de muestra, se copian al host e imprime.

#include <stdio.h>
#include <curand_kernel.h>
#include <time.h>

__global__ void setup_kernel(curandState * state, unsigned long seed) {
int id = threadIdx.x;

/* cada thread tiene la misma semilla, y un diferente número
* de secuencia
*/
curand_init(seed, id, 0, &state[id]);

}

__global__ void generate(curandState* globalState, float *result) {
int ind = threadIdx.x;

// copiar estado a la memoria local para mayor eficiencia
curandState localState = globalState[ind];

// generar número pseudoaleatorio
float r = curand_uniform(&localState);

//copiar state de regreso a memoria global
globalState[ind] = localState;

//almacenar resultados
result[ind] = r;
}

int main(int argc, char** argv) {
int N = 30;

curandState* devStates;
float *devResults;
float hostResults[N];

// reservando espacio para los states PRNG en el device
cudaMalloc(&devStates, N * sizeof(curandState));

// reservando espacio para el vector de resultados en device
cudaMalloc((void**) &devResults, N * sizeof(float));

dim3 tpb(N, 1, 1);

// setup semillas
setup_kernel<<<1, tpb>>>(devStates, time(0));

// generar números aleatorios
generate<<<1, tpb>>>(devStates, devResults);

cudaMemcpy(hostResults, devResults, N * sizeof(float),
cudaMemcpyDeviceToHost);

cudaFree(devStates);
cudaFree(devResults);

for (int i = 0; i < N; i++) {
printf("%f\n", hostResults[i]);
}

return 0;
}


Link: CUDA Curand

viernes, 27 de septiembre de 2013

Modelo de Memoria en CUDA


El modelo de programación en CUDA asume que todos los threads se ejecutan en un dispositivo separado del host que ejecuta la aplicación. Por lo tanto se mantiene implícita la suposición que el host y los dispositivos mantienen sus propios espacios de memoria separados, referidos como la memoria del host (RAM) y la del dispositivo, el cual a su vez está conformado por registros, memoria local, memoria compartida, memoria global, o constantes, como se ve en la siguiente figura:


Cada thread puede:
  • Leer/escribir en registros por thread.
  • Leer/escribir en memoria local por thread.
  • Leer/escribir en memoria compartida por bloque.
  • Leer/escribir en memoria global por grid.
  • Sólo lectura en memoria constante por grid.

Reglas en el manejo de memoria:

  • Actualmente sólo se puede transferir datos desde el host a la memoria global (y memoria constante) y no directamente del host a la memoria compartida.
  • La memoria constante se usa para datos que no cambian (por ejemplo, leídas sólo por la GPU).
  • La memoria compartida llega a tener una velocidad 15x de la memoria global.
  • Los registros tienen velocidad similar a la memoria compartida si lee la misma dirección o no hay conflictos.

Tiempo de vida y alcances de la memoria en CUDA



  • __device__ es opcional cuando se usa con __local__, __shared__, o __constant__
  • Las variables sin identificador residen automáticamente en un registro. Excepto los arrays que residen en memoria local.
  • Las variables escalares residen en registros on-chip de alta velocidad.
  • Las variables compartidas residen en memorias on-chip de alta velocidad.
  • Los arrays locales a un thread y las variables globales residen en memoria off-chip sin caché.
  • Las constantes residen en memoria off-chip sin caché.

3 reglas de la programación en GPGPU

1. Proporcionar los datos a la GPGPU y mantenerlos ahí.

Las GPGPUs son dispositivos que están conectados en un bus PCI Express a la computadora host. El bus PCIe (8 GB/s) es muy lento comparado al sistema de memoria de una GPGPU (160-200 GB/s).

2. Dar suficiente trabajo a la GPGPU.

Debido a que las GPGPU pueden tener un rendimiento en niveles de teraflop, son en muchas ocasiones más rápidos para resolver problemas pequeños de forma más rápida que lo que tarda el host en iniciar el kernel.

3. Enfocarse en el reúso de los datos dentro de la GPGPU para evitar las limitaciones de ancho de banda de la memoria.

Hacer uso de los recursos de memoria internos que dispone CUDA, como son los registros, memoria compartida, y entre otros, para evitar los cuellos de botella en el traspaso de memoria.