En la primer parte de las notas, dedicada a CUDA tuvimos la oportunidad de ver códigos para llevar a cabo la multiplicación de matrices. Ahora es tiempo de hacer lo propio, pero en PyCUDA. Mostramos dos códigos, uno con el algoritmo usual y el segundo con el algoritmo de multiplicación por teselas. A continuación mostramos los dos códigos con sus respectivos tiempos de ejecución.
In [1]:
%%time
import numpy as np
from pycuda import driver, compiler, gpuarray, tools
# -- inicializamos el dispositivo
import pycuda.autoinit
plantilla_codigo_kernel = """
__global__ void MatrizMulKernel(float *a, float *b, float *c)
{
// 2D Thread ID (suponiendo que solo se ejecuta *un* bloque)
int tx = threadIdx.x;
int ty = threadIdx.y;
// Pvalor se usa para guardar el elemento de la matriz
// que es calculado por el thread
float Pvalor = 0;
// Cada thread carga un renglon de M y una columna de N,
// para producir un elemento de P.
for (int k = 0; k < %(TAMANIO_MATRIZ)s; ++k) {
float Aelemento = a[ty * %(TAMANIO_MATRIZ)s + k];
float Belemento = b[k * %(TAMANIO_MATRIZ)s + tx];
Pvalor += Aelemento * Belemento;
}
// Escribe la matriz a la memoria del device;
// cada thread escribe un elemento.
c[ty * %(TAMANIO_MATRIZ)s + tx] = Pvalor;
}
"""
# definimos en tamaño de la matriz
TAMANIO_MATRIZ = 4
# creamos dos matrices aleatorias
a_cpu = np.random.randn(TAMANIO_MATRIZ, TAMANIO_MATRIZ).astype(np.float32)
b_cpu = np.random.randn(TAMANIO_MATRIZ, TAMANIO_MATRIZ).astype(np.float32)
# calculamos una referencia en el CPU para verificar el cálculo en el GPU
c_cpu = np.dot(a_cpu, b_cpu)
# transferimos del host (CPU) al device (GPU)
a_gpu = gpuarray.to_gpu(a_cpu)
b_gpu = gpuarray.to_gpu(b_cpu)
# creamos un arreglo en gpu para guardar el resultado (C = A*B)
c_gpu = gpuarray.empty((TAMANIO_MATRIZ, TAMANIO_MATRIZ), np.float32)
# jalamos el código del kernel
# especificando la constante TAMANIO_MATRIZ
kernel_codigo = plantilla_codigo_kernel % {
'TAMANIO_MATRIZ': TAMANIO_MATRIZ
}
# compilamos el código del kernel
mod = compiler.SourceModule(kernel_codigo)
# jalamos la función del kernel compilado
matrixmul = mod.get_function("MatrizMulKernel")
# llamamos al kernel en la tarjeta
matrixmul(
# entradas
a_gpu, b_gpu,
# salidas
c_gpu,
# bloque de TAMANIO_MATRIZ X TAMANIO_MATRIZ threads
block = (TAMANIO_MATRIZ, TAMANIO_MATRIZ, 1),
)
# imprimimos los resultados
print ("-" * 80)
print ("Matriz A (GPU):")
print (a_gpu.get())
print ("-" * 80)
print ("Matriz B (GPU):")
print (b_gpu.get())
print ("-" * 80)
print ("Matriz C (GPU):")
print (c_gpu.get())
print ("-" * 80)
print ("Diferencia CPU-GPU:")
print (c_cpu - c_gpu.get())
np.allclose(c_cpu, c_gpu.get())
In [2]:
%%time
from __future__ import division
import numpy as np
from numpy import linalg as la
from pycuda import driver, compiler, gpuarray, tools
# -- inicializamos el dispositivo
import pycuda.autoinit
plantilla_kernel_codigo = """
__global__ void MatrizMulKernel(float *A, float *B, float *C)
{
const uint wA = %(TAMANIO_MATRIZ)s;
const uint wB = %(TAMANIO_MATRIZ)s;
// Indice de bloque
const uint bx = blockIdx.x;
const uint by = blockIdx.y;
// Indice de thread
const uint tx = threadIdx.x;
const uint ty = threadIdx.y;
// Indice de la primer submatriz de A procesada por el bloque
const uint aBegin = wA * %(TAMANIO_BLOQUE)s * by;
// Indice de la ultima submatriz de A procesada por el bloque
const uint aEnd = aBegin + wA - 1;
// Tamanio de paso para iterar sobre las submatrices de A
const uint aStep = %(TAMANIO_BLOQUE)s;
// Indice de la primer submatriz de B procesada por el bloque
const uint bBegin = %(TAMANIO_BLOQUE)s * bx;
// Tamanio de paso para iterar sobre las submatrices de B
const uint bStep = %(TAMANIO_BLOQUE)s * wB;
// El elemento de cada submatriz que el calculado por el thread
float Csub = 0;
// Ciclo sobre las submatrices de A y B
for (int a = aBegin, b = bBegin;
a <= aEnd;
a += aStep, b += bStep)
{
// Memoria compartida para la submatriz de A
__shared__ float As[%(TAMANIO_BLOQUE)s][%(TAMANIO_BLOQUE)s];
// Memoria compartida para la submatriz de B
__shared__ float Bs[%(TAMANIO_BLOQUE)s][%(TAMANIO_BLOQUE)s];
// Pasar las matrices de memoria global a compartida
As[ty][tx] = A[a + wA * ty + tx];
Bs[ty][tx] = B[b + wB * ty + tx];
// Sincronizar para asegurarnos de que se cargan completamente
// las matrices antes de comenzar a calcular
__syncthreads();
// Multiplicar las dos matrices
for (int k = 0; k < %(TAMANIO_BLOQUE)s; ++k)
Csub += As[ty][k] * Bs[k][tx];
// Sincronizar para asegurarse de que el calculo anterior
// ha finalizado
__syncthreads();
}
// Escribir la submatriz a la memoria global
const uint c = wB * %(TAMANIO_BLOQUE)s * by + %(TAMANIO_BLOQUE)s * bx;
C[c + wB * ty + tx] = Csub;
}
"""
# definimos el tamaño de la matriz
TAMANIO_MATRIZ = 4
# definimos el tamaño de los bloques y las teselas
TAMANIO_TILE = 2
TAMANIO_BLOQUE = TAMANIO_TILE
# creamos dos matrices cuadradas aleatorias
a_cpu = np.random.randn(TAMANIO_MATRIZ, TAMANIO_MATRIZ).astype(np.float32)
b_cpu = np.random.randn(TAMANIO_MATRIZ, TAMANIO_MATRIZ).astype(np.float32)
# calculamos la referencia
c_cpu = np.dot(a_cpu, b_cpu)
# transferimos memoria del host al device
a_gpu = gpuarray.to_gpu(a_cpu)
b_gpu = gpuarray.to_gpu(b_cpu)
# creamos un arreglo en blanco para guardar el resultado
c_gpu = gpuarray.empty((TAMANIO_MATRIZ, TAMANIO_MATRIZ), np.float32)
# jalamos el código del kernel especificando las constantes
kernel_codigo = plantilla_kernel_codigo % {
'TAMANIO_MATRIZ': TAMANIO_MATRIZ,
'TAMANIO_BLOQUE': TAMANIO_BLOQUE,
}
# compilamos el kernel
mod = compiler.SourceModule(kernel_codigo)
# jalamos la función del kernel
matrixmul = mod.get_function("MatrizMulKernel")
# llamamos al kernel
matrixmul(
# entradas
a_gpu, b_gpu,
# salidas
c_gpu,
# malla de varios bloques
grid = (TAMANIO_MATRIZ // TAMANIO_TILE, TAMANIO_MATRIZ // TAMANIO_TILE),
# bloque de varios threads
block = (TAMANIO_TILE, TAMANIO_TILE, 1),
)
# Imprimimos los resultados
print "-" * 80
print "Matrix A (GPU):"
print a_gpu.get()
print "-" * 80
print "Matrix B (GPU):"
print b_gpu.get()
print "-" * 80
print "Matrix C (GPU):"
print c_gpu.get()
print "-" * 80
print "Diferencia CPU-GPU:"
print c_cpu - c_gpu.get()
print "Norma L2:", la.norm(c_cpu - c_gpu.get())
np.allclose(c_cpu, c_gpu.get())
Como podemos apreciar, la multiplicación por teselas es más eficiente.
Algo que es importante remarcar es el hecho de que uno estaría tentado a usar double en lugar de float en el kernel, sin embargo aún hay problemas para utilizar variables double en algunas tarjetas (tal es el caso de la tarjeta en que corremos esto).
En las notas de Roberto Zamora podemos ver operaciones más complicadas.