Multiplicando matrices

Introducción al Cálculo Numérico
en procesadores gráficos

Multiplicando matrices

GPGPU @ CABIB

1

GPGPU @ CABIB

Icnpg (2012): Multiplicando matrices

El Problema
Vamos a pensar en el producto de dos matrices A y B de forma tal que C = A × B. Entonces un elemento Cij de
C es
N−1
∑
Cij =
Aik Bkj
k=0

donde N es el número de columnas de A, o el número de filas de B.Asumimos también que los índices de las
matrices comienzan en cero. El algoritmo para calcular C es sencillo: debemos recorrer toda la matriz C con los
índices i y j, y para cada elemento de C tenemos que hacer la suma indicada.

El Código
El código tiene dos archivos fuentes, matrixMul.cu y matrixMul.h. También van a encontrar el Makefile correspondiente. Los tamaños de las matrices estándefinidos en matrixMul.h con \verb#define…=, al igual que
el tamaño de bloque a utilizar BLOCK_SIZE y un parámetro SCALE que permite agrandar el tamaño de las matrices
sencillamente. Los kernels están definidos en matrixMul.cu siguiendo el orden de esta guía. Tratamos de mantener el código lo más sencillo posible, aunque quizás para probar performance haya que pasar por la compilación
frecuentemente.1. Indexado
La primer dificultad que uno encuentra al tratar con este problema es que la indexación de las matrices en la documentación de CUDA (Tutorials, clases, etc. etc.) es al revés que en matemática. Veamos.
Supongamos que A tiene hA filas y wA columnas. En matemática definimos Aij al elemento de la fila i y la columna
j. Sin embargo, aquí la notación usual es A(x, y) donde x es el eje'horizontal' de la matriz (cuenta las columnas!),
mientras que y es el eje vertical (cuenta filas!). De esta forma, por ejemplo, A21 = A(1, 2) es el elemento de la tercera
fila (contamos desde cero a hA-1) y de la segunda columna (contamos desde cero a wA-1).
Por otra parte, vamos a manejar las matrices como arrays de una dimensión. Entonces, al elemento A(x, y) se
llega recorriendo lasprimeras y filas, y luego avanzando x columnas:
A(x, y) = A[y * wA + x]
o también recorriendo las primeras x columnas, y avanzando y filas:
A(x, y) = A[y + x * hA]
Cualquiera de las dos alternativas es correcta, pero recorrer la matriz por filas o columnas puede cambiar drásticamente la eficiencia del código.

2. Multiplicando matrices: Host
Consideremos ahora también a B tiene hB filas y wBcolumnas. Para poder multiplicar ambas matrices, el número de
columnas wA de A tiene que ser igual al número de filas hB de B. Por lo tanto, el resultado de A × B es una matriz
C de hA filas y wB columnas. Estas condiciones están definidas en matrixMul.h:
matrixMul.h
26
27
28

/ / Basic Matrix dimensions
/ / ( chosen as m u l t i p l e s of the thread block s i z e f o r s i m p l i c i t y )# d e f i n e WA ( 4 ∗ BLOCK SIZE ∗ SCALE ) / / M a t r i x A w i d t h

2. Multiplicando matrices: Host continua en la página siguiente. . .

2/7

GPGPU @ CABIB

29
30
31
32
33

#define
#define
#define
#define
#define

HA ( 6 ∗ BLOCK SIZE
WB ( 4 ∗ BLOCK SIZE
HB WA / / M a t r i x B
WC WB / / M a t r i x C
HC HA / / M a t r i x C

Icnpg (2012): Multiplicando matrices∗ SCALE ) / / M a t r i x A h e i g h t
∗ SCALE ) / / M a t r i x B w i d t h
height
width
height

de manera tal que podemos usar hC y wC como el número de filas y columnas de C, respectivamente.
El código siguiente tiene la implementación del producto de matrices en el Host. A la función computeGold le
pasamos las tres matrices, el número de filas y columnas de A, y el número de columnas deB. El índice i recorre las
filas de A, el índice j recorre las columnas de B.
matrixMul.h
83
84

85
86
87
88
89
90
91
92
93
94
95
96

void
computeGold ( f l o a t ∗ C , c o n s t f l o a t ∗ A , c o n s t f l o a t ∗ B , u n s i g n e d i n t hA , u n s i g n e d i n t wA, u n s i g n e d
i n t wB )
{
f o r ( u n s i g n e d i n t i = 0 ; i < hA ; ++ i )
f o r ( u n s...