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²).