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

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.

Embarrasingly Parallel Algorithms

Llevan éste nombre los algoritmos más sencillos de adaptar su implementación en una GPGPU. Para conocerlos mejor, nombraremos las siguientes características:

• También se les conoce como algoritmos naturalmente paralelos.
• Son los algoritmos paralelos más simples debido que casi no requieren comunicación entre procesos.
• Cada proceso puede realizar sus cálculos de forma propia sin necesidad de comunicarse con otros.
• Pueden requerir alguna partición inicial de los datos, o juntar los datos resultantes al final, aunque no siempre.

Caso ideal:

• Todos los subproblemas o tareas son definidas antes de que el cómputo inicie.
• Todas las sub-soluciones son almacenadas en localidades de memoria independientes (variables, elementos de arrays).
• Por lo tanto, el cómputo de cada sub-solución es completamente independiente.
• Si el cómputo requiere alguna comunicación inicial o final, lo llamaremos Nearly embarrasingly parallel.

Algunos ejemplos:

• Renderizado de gráficos para computadora.
• Algoritmos genéticos
• Simulaciones de Monte-Carlo
• Conjuntos de Mandelbrot (a.k.a. Fractales)

Link interesante

Cuda Teaching Center for High Performance Computing

de Wake Forrest University, San Diego, California, USA.


http://users.wfu.edu/choss/CUDA/

Multiplicar matrices en CUDA

Programa para calcular la multiplicación de matrices en CUDA.

#include <stdio.h>
#define N 16

void matrixMultCPU(int a[N][N], int b[N][N], int c[N][N]) {
 int n,m;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
int sum = 0;
for (int k = 0; k < N; k++) {
m = a[i][k];
n = b[k][j];
sum += m * n;
}
c[i][j] = sum;
}
}
}

__global__ void matrixMultGPU(int *a, int *b, int *c) {
int k, sum = 0;
int col = threadIdx.x + blockDim.x * blockIdx.x;
int fil = threadIdx.y + blockDim.y * blockIdx.y;

if (col < N && fil < N) {
for (k = 0; k < N; k++) {
sum += a[fil * N + k] * b[k * N + col];
}
  c[fil * N + col] = sum;
}
}

int main() {
int a[N][N], b[N][N], c[N][N];
int *dev_a, *dev_b, *dev_c;
int cont,i,j;

/* inicializando variables con datos foo*/
for (i = 0; i < N; i++) {
cont = 0;
for (j = 0; j < N; j++) {
a[i][j] = cont;
b[i][j] = cont;
cont++;
}
}

int size = N * N * sizeof(int);

cudaMalloc((void **) &dev_a, size);
cudaMalloc((void **) &dev_b, size);
cudaMalloc((void **) &dev_c, size);

cudaMemcpy(dev_a, a, size, cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, size, cudaMemcpyHostToDevice);

dim3 dimGrid(1, 1);
dim3 dimBlock(N, N);

matrixMultGPU<<<dimGrid, dimBlock>>>(dev_a, dev_b, dev_c);

cudaMemcpy(c, dev_c, size, cudaMemcpyDeviceToHost);

cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);

// imprimiendo
for (int y = 0; y < N; y++) {
for (int x = 0; x < N; x++) {
printf("[%d][%d]=%d ", y, x, c[y][x]);
}
printf("\n");
}

return 0;

}

Notas importantes del código:

• Está codificada la función para hacer la multiplicación en una CPU y una GPGPU para comprobar las diferencias.
• El grid de bloques es de 1D, sólo 1 bloque.
• El número de threads en el bloque es cuadrado, equivalente al número de elementos lineales de la matriz (16 x 16).
• Es muy importante no olvidar el offset de índices para los elementos de la matriz.
• Como se ha definido cada bloque de forma bidimensional, se calculan índices de columnas y filas (después se convertirán a un índice lineal para hacer las operaciones)
• Se lanzan 256 threads, la forma más sencilla de verlo es que en cada uno de ellos se generará el cálculo de su correspondiente elemento en la matriz resultante. Para poder generar éste resultado, es necesario hacer un ciclo hasta N (16 en éste caso) para hacer la sumatoria de las multiplicaciones necesarias.
• La multiplicación en la CPU se realiza con ciclos anidados,dichas multiplicaciones y sumas están en O(N³).
• En el caso de la multiplicación en la GPU, debido a que un único thread se utiliza para calcular el valor de c(i,j), ésta solución es de tipo O(N²).

Suma de matrices (método 2)

Código para sumar una matriz en CUDA

#include "stdio.h"
#define columnas 300
#define filas 300
__global__ void add(int *a, int *b, int *c) {

int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int i = (columnas * y) + x;

c[i] = a[i] + b[i];
}

int main() {
int cont = 0;
int i, j;
// matrices en host
int a[filas][columnas], b[filas][columnas], c[filas][columnas];

// matrices en GPGPU
int *dev_a, *dev_b, *dev_c;

cudaMalloc((void **) &dev_a, filas * columnas * sizeof(int));
cudaMalloc((void **) &dev_b, filas * columnas * sizeof(int));
cudaMalloc((void **) &dev_c, filas * columnas * sizeof(int));

/* inicializando variables con datos foo*/
for (i = 0; i < filas; i++) {
cont = 0;
for (j = 0; j < columnas; j++) {
a[i][j] = cont;
b[i][j] = cont;
cont++;
}
}
cudaMemcpy(dev_a, a, filas * columnas * sizeof(int),cudaMemcpyHostToDevice);
cudaMemcpy(dev_b, b, filas * columnas * sizeof(int),cudaMemcpyHostToDevice);

// definiendo grid
dim3 grid(columnas, filas);

// grid del tamaño de la matriz, con un thread por bloque
add<<<grid, 1>>>(dev_a, dev_b, dev_c);

cudaMemcpy(c, dev_c, filas * columnas * sizeof(int), cudaMemcpyDeviceToHost);

// imprimiendo
for (int y = 0; y < filas; y++)
{
for (int x = 0; x < columnas; x++) {
printf("[%d][%d]=%d ", y, x, c[y][x]);
}
printf("\n");
}
return 0;
}

Notas importantes del código:

• El código anterior realiza la suma de una matriz cuadrada de 300 x 300.
• La diferencia principal con el código del método 1, es que el grid de bloques mapea la matriz definida, y el número de threads por bloque es sólo de 1.
• Los demás cálculos son equivalentes.

Suma de Matrices (método 1)

Código para sumar una matriz en CUDA

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <cuda.h>
#include <math.h>

#define T 10 // max threads x bloque
#define N 300


__global__ void sumaMatrices(int *m1, int *m2, int *m3) {

int col = blockIdx.x * blockDim.x + threadIdx.x;
int fil = blockIdx.y * blockDim.y + threadIdx.y;

int indice = fil * N + col;


if (col < N && fil < N) {
// debido a que en los últimos bloques no se realizan todos los threads
m3[indice] = m1[indice] + m2[indice];
}
}

int main(int argc, char** argv) {

int m1[N][N];
int m2[N][N];
int m3[N][N];
int i, j;
int c = 0;

/* inicializando variables con datos foo*/
for (i = 0; i < N; i++) {
c = 0;
for (j = 0; j < N; j++) {
m1[i][j] = c;
m2[i][j] = c;
c++;
}
}

int *dm1, *dm2, *dm3;

cudaMalloc((void**) &dm1, N * N * sizeof(int));
cudaMalloc((void**) &dm2, N * N * sizeof(int));
cudaMalloc((void**) &dm3, N * N * sizeof(int));

// copiando memoria a la GPGPU
cudaMemcpy(dm1, m1, N * N * sizeof(int), cudaMemcpyHostToDevice);
cudaMemcpy(dm2, m2, N * N * sizeof(int), cudaMemcpyHostToDevice);

// cada bloque en dimensión x y y tendrá un tamaño de T Threads
dim3 dimThreadsBloque(T, T);

// Calculando el número de bloques en 1D
float BFloat = (float) N / (float) T;
int B = (int) ceil(BFloat);

// El grid tendrá B número de bloques en x y y
dim3 dimBloques(B, B);

// Llamando a ejecutar el kernel
sumaMatrices<<<dimBloques, dimThreadsBloque>>>(dm1, dm2, dm3);

// copiando el resultado a la memoria Host
cudaMemcpy(m3, dm3, N * N * sizeof(int), cudaMemcpyDeviceToHost);
//cudaMemcpy(m2, dm2, N * N * sizeof(int), cudaMemcpyDeviceToHost);

cudaFree(dm1);
cudaFree(dm2);
cudaFree(dm3);

printf("\n");

for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
printf(" [%d,%d]=%d", i, j, m3[i][j]);

}
printf("\n\n");

}
printf("\nB = %d", B);
printf("\n%d, %d",dimBloques.x, dimBloques.y);
printf("\n%d, %d",dimThreadsBloque.x, dimThreadsBloque.y);
return (EXIT_SUCCESS);
}

Notas importantes:

• Realiza la suma de una matriz cuadrada de enteros de 300 x 300 (90,000 elementos en total).
• Define cada bloque con Threads de 2D con valores de 10, es decir 100 threads por bloque.
• Define un grid de bloques en 2D. Para que alcancen los elementos de la matriz el grid será de 30 x 30 (900 bloques).
• Los elementos de la matriz se toman como un único vector. C se encarga de hacerlo parecer como una matriz. Por lo tanto es necesario calcular el desplazamiento que se genera al aumentar la fila del elemento que se esté calculando:

int col = blockIdx.x * blockDim.x + threadIdx.x;
int fil = blockIdx.y * blockDim.y + threadIdx.y;

int indice = fil * N + col;

SumaVectores.cu

El siguiente código en CUDA-C realiza una suma básica de dos vectores, y almacena el resultado en un tercer vector.

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <cuda.h>
#include <math.h>

#define T 1024 // max threads x bloque
#define N 100000

__global__ void sumaVector(int *v1, int *v2, int *v3) {

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

if (tid < N) {
// debido a que en el último bloque no se realizan todos los threads
v3[tid] = v1[tid] + v2[tid];
}
}

int main(int argc, char** argv) {

int v1[N];
int v2[N];
int v3[N];
int i;

/* inicializando variables con datos foo*/
for (i = 0; i < N; i++) {
v1[i] = i;
v2[i] = i;
v3[i] = 0;
}

for (i = 0; i < N; i++) {
printf(" %d %d \t", v1[i], v2[i]);
if (i % 20 == 0)
printf("\n");

}

int *dv1, *dv2, *dv3;

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

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

/* Calculando el número de bloques*/
float BFloat = (float) N / (float) T;
int B = (int) ceil(BFloat);

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

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

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

printf("\n");

for (i = 0; i < N; i++) {
printf(" %d=%d", i, v3[i]);
if (i % 40 == 0)
printf("\n");

}

return (EXIT_SUCCESS);
}

Características importantes a tomar en cuenta

• Este cálculo se realiza en dimensión 1D.
• Debido a que el número de threads por dimensión en un bloque es 1024, se define ese valor como T. 
• N representa el tamaño del vector, en este caso se definió en 100 mil.
• Debido al tamaño de N, es necesario dividir el trabajo en diferentes bloques, por lo que se realiza el cálculo para la variable B (en este caso seria 100000/1024 =97.65, convertido a 98 bloques).
• Cada bloque es reservado completamente, aunque parte del último bloque no se utilice.

Arquitectura de una GPGPU NVIDIA Tesla C2075

Número máximo de threads por bloque: 1024

Tamaño máximo de las dimensiones x, y, z de un bloque de threads: 1024 x 1024 x 64

Tamaño máximo de cada dimensión del grid de bloques de threads: 65535 x 65535 x 65535

NVIDIA Nsight Eclipse Edition

Nsight Eclipse es un IDE basado en Eclipse que permite editar, construir y debuggear aplicaciones en CUDA-C. Viene incluido por default en la versión 5.5 de Cuda Toolkit.

Sólo es necesario escribir en la consola

$nsight

Y la aplicación se iniciará:

Características de NVIDIA Nsight Eclipse Edition

Instalar Fedora 18 & CUDA 5.5

A continuación se describen los pasos necesarios para instalar CUDA 5.5 en Fedora 18 kernel version 3.6.10-4.fc18.x86_64

Instalar paquetes necesarios

sudo yum install kernel-devel
sudo yum install gcc-c++
yum update audit

Repositorio CUDA

Descargar el repositorio para Fedora 18 desde  CUDA Donwload e instalarlo en una terminal:

sudo rpm -Uhv cuda-repo-fedora18-5.5-0.x86_64.rpm

Instalar el controlador propietario de Video

El controlador instalado por default "nouveau" es incompatible con el Toolkit de CUDA, por lo que es necesario reemplazarlo:

sudo yum remove xorg-x11-drv-nouveau
sudo yum install nvidia-settings nvidia-kmod xorg-x11-drv-nvidia

Para prevenir que el controlador nouveau se active accidentalmente después, se debe editar el archivo /etc/default/grub. En la línea GRUB_CMDLINE_LINUX_DEFAULT añadir:

"rdblacklist=nouvear nouveau.modeset=0"

Reiniciar el sistema

Ahora ya es posible reiniciar el sistema operativo, debiendo trabajar correctamente.

Instalar CUDA Toolkit

En una terminal instalar CUDA Toolkit mediante:

sudo yum install cuda

* Descargará aproximadamente 700 Mb de datos en la instalación

Añadir variables de ambiente

Editar el archivo .bashrc del home, añadiendo:

export CUDA_HOME=/usr/local/cuda
export LD_LIBRARY_PATH=${CUDA_HOME}/lib64

PATH=${CUDA_HOME]/bin:${PATH}
export PATH

Comprobar la instalación

La instalación debe estar completa, para comprobar es posible realizar un ejemplo básico de Hello world y compilarlo.

O también dentro de /usr/local/cuda/samples compilar usando Make, y correr /usr/local/cuda/samples/bin/x86_64/linux/release/.deviceQuery.

Definiendo las estructuras Grid/Block para un kernel

Para llamar a ejecutar un kernel es necesario proveer una configuración de ejecución, esto es, las dimensiones del grid y del bloque con que se ejecutará dicho kernel. Esta información se proporciona mediante dos estructuras claves:

• Número de bloques en cada dimensión
• Threads por bloque en cada dimensión

Sintaxis de llamada a un kernel:

      nombreKernel <<< B, T >>> (arg1, arg2, ... argN);

donde:

• B es una estructura que define el número de bloques en el grid para cada dimensión (1D o 2D)
• T es una estructura que define el número de threads en un bloque para cada dimensión (1D, 2D o 3D).

Si se desea definir una estructura de 1D, se puede usar un entero para B y T en:

       nombreKernel <<< B, T >>> (arg1, arg2, ... argN);

donde:

• B es un entero que define un grid de dimensión 1D
• T es un entero que define un bloque 1D de ese tamaño.

por ejemplo:

       nombreKernel <<< 1, 100 >>> (arg1, arg2, ... argN);

dim3 en CUDA

dim3 es una estructura de datos utilizada comúnmente en CUDA para pasar al kernel la configuración de ejecución (es decir, las dimensiones del grid y del bloque).

Por ejemplo:

dim3 grid(512);                        //512 x 1 x 1

dim3 block(1024, 1024);              //1024 x 1024 x 1

kernelEjemplo <<grid, block >>>();

Puntos a tomar en cuenta de dim3:

• dim3 tiene 3 elementos x, y y z.
• dim3 es una estructura definida en %CUDA_INC_PATH%/vector_types.h
• dim3 en código C puede ser inicializado como dim3 grid = {512, 512, 1};
• dim3 en código C++ puede ser inicializado como dim3 grid(512,512,1);
• No es necesario proporcionar los 3 elementos. Si alguno de ellos falta, será inicializado por defecto con valor 1.

Identificación de un thread

CUDA identifica a cada thread de forma única. Para ello se apoya de ciertas variables incluidas en el driver:

• blockIdx.x, blockIdx.y, blockIdx.z
• Retornan el valor del ID del bloque en el eje x, y o z, respectivamente, del bloque que esté ejecutando el código dado.
• blockDim.x, blockDim.y, blockDim.z
• Retornan el valor de la dimensión del bloque, es decir, el número de threads en un bloque en el eje x, y o z.
• threadIdx.x, threadIdx.y, threadIdx.z
• Retornan el ID del thread en el eje x, y o z que está siendo ejecutado por el GPGPU en un bloque en particular.

De ésta forma, se puede expresar una colección de bloques, y una colección de threads dentro de un bloque, como un array 1D, 2D o 3D, lo cual puede ser de ayuda con información en 2D o 3D.

Ejemplo:

Suponer un grid de 1D (un grid se mapea directamente con una GPGPU), y una arquitectura de bloques de 1D con 4 bloques y cada uno con 8 threads.

Para identificar el ID Global del thread 26:

GlobalThreadID = blockIdx.x * blockDim.x + threadIdx.x
GlobalThreadID =         3       *         8         +      2
GlobalThreadID = 26       

Conceptos básicos en CUDA

• Una GPGPU corre de forma asíncrona al CPU del host
• Una GPGPU corre en un espacio de memoria separado del CPU del host
•  La velocidad de la memoria del host es 8-20 GB/s, mientras que la velocidad de la memoria en una GPGPU es aproximadamente 160-200 GB/s.
• Los programas en CUDA usan kernels
• Un kernel es una subrutina llamada desde el host que se ejecutara en la GPGPU.
• Un kernel no es una función, ya que no puede retornar un valor.
• Un kernel está definido con una especificador declarativo __global__ , el cual le indica al compilador que el kernel es llamable desde el procesador del host.
• Las llamadas a kernels son asíncronas
•  El host encola un kernel para su ejecución en el GPGPU y no lo espera a que finalice.
• En algún momento después el kernel corre en el GPGPU
• Para mayor eficiencia, se pueden crear pipelines para encolar kernels y mantener el GPGPU ocupado tanto sea posible.
• Más adelante, puede ser necesario determinar cuando el kernel o la pipeline ha sido completada, por lo que se utilizan comúnmente dos mecanismos de sincronización:
• cudaThreadSynchronize(), hace que el host se detenga y espere a que todos los kernels encolados se ejecuten
• Ejecutar una transferencia de datos con bloque mediante  cudaMemcpy(), ya que incluye dentro la llamada a cudaThreadSyncrhonize().
• La unidad básica de trabajo de una GPGPU es un thread. 
• Cada thread actúa como si tuviera su propio procesador con identidad y registros separados
• El hardware de la GPGPU (GPU thread scheduler) define cuando un grupo de threads que pueden correr concurrentemente. También tiene la habilidad para intercambiar entre threads de forma rápida y transparente desde el punto de vista del software.
• Existen algunas instrucciones que permiten que los threads se comuniquen a través de espacios de memoria compartida de CUDA.
• Una configuración de ejecución define el número de threads que ejecutará el kernel y el arreglo del grid computacional en 1D, 2D o 3D. Esta información se encierra en triples brackets "<<<" y ">>>" después del nombre del kernel y antes de la lista de parámetros entre paréntesis.
• La región más grande en una GPGPU es llamada memoria global. 
• La memoria global está sujeta a reglas de fusión que combinan múltiples transacciones de memoria en una sola carga en funciones que pueden obtener una transferencia más alta desde y hacia la memoria.
• En general, el mayor rendimiento de la memoria ocurre cuando se accesa de manera en trozos de 128 bytes consecutivos.
• Existen otras formas de memoria accesible para el programador de una GPGPU, como son constant, cache, shared, local, texture y register memory.
• Aunque el ancho de banda de la memoria global es de 160-200 GB/s, es lento comparado con el rendimiento de 1 teraflop que puede alcanzar. Por lo tanto es esencial hacer reuso de datos para alcanzar el máximo rendimiento.

NVIDIA Tesla C2075

./deviceQuery de NVIDIA Tesla C2075


[cuda@scamlab release]$ ./deviceQuery
[deviceQuery] starting...

./deviceQuery Starting...

 CUDA Device Query (Runtime API) version (CUDART static linking)

Found 1 CUDA Capable device(s)

Device 0: "Tesla C2075"
  CUDA Driver Version / Runtime Version          4.2 / 4.2
  CUDA Capability Major/Minor version number:    2.0
  Total amount of global memory:                 5375 MBytes (5636292608 bytes)
  (14) Multiprocessors x ( 32) CUDA Cores/MP:    448 CUDA Cores
  GPU Clock rate:                                1147 MHz (1.15 GHz)
  Memory Clock rate:                             1566 Mhz
  Memory Bus Width:                              384-bit
  L2 Cache Size:                                 786432 bytes
  Max Texture Dimension Size (x,y,z)             1D=(65536), 2D=(65536,65535), 3D=(2048,2048,2048)
  Max Layered Texture Size (dim) x layers        1D=(16384) x 2048, 2D=(16384,16384) x 2048
  Total amount of constant memory:               65536 bytes
  Total amount of shared memory per block:       49152 bytes
  Total number of registers available per block: 32768
  Warp size:                                     32
  Maximum number of threads per multiprocessor:  1536
  Maximum number of threads per block:           1024
  Maximum sizes of each dimension of a block:    1024 x 1024 x 64
  Maximum sizes of each dimension of a grid:     65535 x 65535 x 65535
  Maximum memory pitch:                          2147483647 bytes
  Texture alignment:                             512 bytes
  Concurrent copy and execution:                 Yes with 2 copy engine(s)
  Run time limit on kernels:                     Yes
  Integrated GPU sharing Host Memory:            No
  Support host page-locked memory mapping:       Yes
  Concurrent kernel execution:                   Yes
  Alignment requirement for Surfaces:            Yes
  Device has ECC support enabled:                Yes
  Device is using TCC driver mode:               No
  Device supports Unified Addressing (UVA):      Yes
  Device PCI Bus ID / PCI location ID:           3 / 0
  Compute Mode:
     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >

CUDA Hello World

Programa Hello World.

Imprime "Hello" a partir de un vector, y realiza una suma de vectores para crear el mensaje "World!".



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

const int N = 16;
const int blocksize = 16;


/**
* Suma vector 'b' al 'a'
*/

__global__ void hello(char *a, int *b) {
    a[threadIdx.x] += b[threadIdx.x];
}

/**
*
*/
int main(int argc, char** argv) {

    char a[N] = "Hello \0\0\0\0\0\0";
    int b[N] = {15, 10, 6, 0, -11, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

    char *ad; // d for device variables
    int *bd;
    const int csize = N * sizeof (char);
    const int isize = N * sizeof (int);

    printf("%s", a);

    cudaMalloc((void**) &ad, csize);
    cudaMalloc((void**) &bd, isize);
    cudaMemcpy(ad, a, csize, cudaMemcpyHostToDevice);
    cudaMemcpy(bd, b, isize, cudaMemcpyHostToDevice);

    dim3 dimBlock(blocksize, 1);
    dim3 dimGrid(1, 1);
    hello << <dimGrid, dimBlock >> >(ad, bd);
    cudaMemcpy(a, ad, csize, cudaMemcpyDeviceToHost);
    cudaFree(ad);
    cudaFree(bd);

    printf("%s\n", a);
    return (EXIT_SUCCESS);
}

¿Cómo compilar CUDA desde consola?

Instrucciones para compilar un programa en CUDA desde consola en Linux

1. Crear archivo "file.cu"
2. Compilar
• nvcc -c file.cu
3. Linkear
• g++ -o nameExecutable -L/usr/local/cuda/lib64/ -lcudart -lcuda file.o
4. Ejecutar
• ./nameExecutable