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.