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.



viernes, 20 de septiembre de 2013

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)


jueves, 19 de septiembre de 2013

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(N3).
  • 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(N2).


miércoles, 4 de septiembre de 2013

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.





martes, 3 de septiembre de 2013

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


lunes, 2 de septiembre de 2013

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.