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