Content |
GEMM
Algorithmus
Die Matrizen \(A, B\) und \(C\) werden in Blöcke unterteilt:
-
\(A\) besteht aus maximal \(m_c \times k_c\) großen Blöcken,
-
\(B\) besteht aus maximal \(k_c \times n_c\) großen Blöcken und
-
\(C\) besteht aus maximal \(m_c \times n_c\) großen Blöcken.
Man kann zu Berechnung von \(C\) folgenden Algorithmus verwenden:
-
Sei \(M = \lfloor \frac{m+m_c-1}{m_c} \rfloor\)
-
Sei \(K = \lfloor \frac{k+k_c-1}{k_c} \rfloor\)
-
Sei \(N = \lfloor \frac{n+n_c-1}{n_c} \rfloor\)
-
Für \(i=1, \dots, M\)
-
Für \(l=1, \dots, K\)
-
Kopiere den Block \(B_{i,l}\) in den Puffer \(\tilde{B}\)
-
Für \(j=1, \dots, N\)
-
Kopiere den Block \(A_{i,l}\) in den Puffer \(\tilde{A}\)
-
Falls \(l=1\) ist berechne
-
\(C_{i,j} \leftarrow \beta C_{i,j} +\alpha \cdot \tilde{A} \tilde{B}\)
sonst
-
\(C_{i,j} \leftarrow C_{i,j} +\alpha \cdot \tilde{A} \tilde{B}\)
-
-
-
-
Bemerkungen
-
Wie groß sind die Blöcke?
-
Wieso ist die Fallunterscheidung für \(l=1\) nötig?
-
Man könnte die for-Schleifen auch anders ordnen, z.B.
-
Erst über \(i\), dann über \(j\) und dann über \(k\) iterieren.
-
Erst über \(k\), dann über \(i\) und dann über \(j\) iterieren.
Welche Nachteile hätten diese Strategien? Wieviele Varianten gibt es?
-
Beispiel
Betrachten wir anhand eines kleinen Beispiel die Multiplikation von solchen Block-Matrizen:
\[\begin{eqnarray*}\left(\begin{array}{ccc}C_{1,1} & C_{1,2} & C_{1,3} \\C_{2,1} & C_{2,2} & C_{2,3} \\\end{array}\right)&\leftarrow&\beta\left(\begin{array}{ccc}C_{1,1} & C_{1,2} & C_{1,3} \\C_{2,1} & C_{2,2} & C_{2,3} \\\end{array}\right)+ \alpha\left(\begin{array}{cccc}A_{1,1} & A_{1,2} & A_{1,3} & A_{1,4} \\A_{2,1} & A_{2,2} & A_{2,3} & A_{2,4} \\\end{array}\right)\left(\begin{array}{cccc}B_{1,1} & B_{1,2} & B_{1,3}\\B_{2,1} & B_{2,2} & B_{2,3}\\B_{3,1} & B_{3,2} & B_{3,3}\\B_{4,1} & B_{4,2} & B_{4,3}\\\end{array}\right)\end{eqnarray*}\]Dann gilt
\[C_{i,j} \leftarrow \beta\cdot C_{i,j} + \alpha\cdot A_{i,1} B_{1,j} + \alpha\cdot A_{i,2} B_{2,j} + \alpha\cdot A_{i,3} B_{3,j} + \alpha\cdot A_{i,4} B_{4,j}\]Also
\[\begin{eqnarray*}C_{i,j} &\leftarrow& \beta \cdot C_{i,j} + \alpha \cdot A_{i,1} B_{1,j} \\C_{i,j} &\leftarrow& \quad\, C_{i,j} + \alpha \cdot A_{i,2} B_{2,j} \\C_{i,j} &\leftarrow& \quad\, C_{i,j} + \alpha \cdot A_{i,3} B_{3,j} \\C_{i,j} &\leftarrow& \quad\, C_{i,j} + \alpha \cdot A_{i,4} B_{4,j}\end{eqnarray*}\]Die einzelnen Produkte dieser Matrix Blöcke werden mit dem Macro-Kernel durchgeführt. Hierbei werden statt der Blöcke \(A_{i,l}\) und \(B_{l,j}\) aber jeweils deren Kopien \(\tilde{A}, \tilde{B}\) mit eventuellem Padding benutzt. Das Produkt \(\alpha \cdot \tilde{A} \tilde{B}\) hat also immer die Dimension \(m_c \times n_c\). Falls \(C_{i,j}\) aber ein Rand-Block ist kann dessen Größe aber kleiner sein. In diesem Fall verwendet man auch hier einen Puffer \(\tilde{C}\) und berechnet zunächst \(\tilde{C} \leftarrow \alpha \cdot \tilde{A} \tilde{B}\). Man kann also auch hier den gleichen Macro-Kernel benutzen. Anschliessend kann man \(C_{i,j}\) mit dem relevanten Teil von \(\tilde{C}\) aufdatieren, d.h. im wesentlichen \(C_{i,j} \leftarrow \beta C_{i,j} + \tilde{C}\) berechnen.
Sonderfälle in der BLAS Spezifikation
Für \(A \in \mathbb{R}^{m \times n}\) und \(B \in \mathbb{R}^{k \times n}\) mit \(k=0\) soll das Produkt \(C \leftarrow \beta C + \alpha A B\) als Skalierung \(C \leftarrow \beta C\) interpretiert werden.
Aufgabe
-
Implementiert das Matrix-Matrix Produkt \(C \leftarrow \beta C + \alpha \cdot AB\). Testet es zunächst als stand-alone Einheit.
-
Baut es anschliessend in src/level3/ ein:
-
Alle Funktionen werden wieder mit ULMBLAS(..) ummantelt.
-
In src/level3/dgemm_nn.h wird die Funktion ULMBLAS(dgemm_nn) deklariert. Dieses Header-File wird von src/level3/dgemm.c eingebunden.
-
Die Implementierungen wandern in src/level3/dgemm_nn.c.
-
-
Lasst den Benchmark durchlaufen:
-
Testet was passiert, wenn man in src/ mit -O3 übersetzt (dazu src/Makefile anpassen).
-
#include <stdio.h>
#define M 14
#define K 15
#define N 16
#define MC 8
#define KC 11
#define NC 12
#define MR 4
#define NR 6
double A[M*K];
double B[K*N];
double C[M*N];
double _A[MC*KC];
double _B[KC*NC];
double _C[MC*NC];
void
initMatrix(int m, int n, double *X, int ldX, int counter)
{
// ... your code here ...
}
void
printMatrix(int m, int n, const double *X, int ldX)
{
// ... your code here ...
}
void
pack_MRxk(int k, const double *A, int incRowA, int incColA, double *buffer)
{
// ... your code here ...
}
void
pack_A(int mc, int kc, const double *A, int incRowA, int incColA,
double *buffer)
{
// ... your code here ...
}
void
pack_kxNR(int k, const double *B, int incRowB, int incColB, double *buffer)
{
// ... your code here ...
}
void
pack_B(int kc, int nc, const double *B, int incRowB, int incColB,
double *buffer)
{
// ... your code here ...
}
void
dgemm_micro_kernel(int kc,
double alpha, const double *A, const double *B,
double beta,
double *C, int incRowC, int incColC)
{
// ... your code here ...
}
void
dgemm_macro_kernel(int mc,
int nc,
int kc,
double alpha,
double beta,
double *C,
int incRowC,
int incColC)
{
// ... your code here ...
}
void
dgescal(int m,
int n,
double alpha,
const double *X,
int incRowX,
int incColX)
{
// ... your code here: compute X <- alpha*X
}
void
dgeaxpy(int m,
int n,
double alpha,
const double *X,
int incRowX,
int incColX,
double *Y,
int incRowY,
int incColY)
{
// ... your code here: compute Y <- Y + alpha*X
}
void
dgemm_nn(int m,
int n,
int k,
double alpha,
const double *A,
int incRowA,
int incColA,
const double *B,
int incRowB,
int incColB,
double beta,
double *C,
int incRowC,
int incColC)
{
// ... your code here ...
}
int
main()
{
int i, j, mc, kc, nc;
initMatrix(M, K, A, M, 1);
initMatrix(K, N, B, K, M*K+1);
printf("A = \n");
printMatrix(M, K, A, M);
printf("B = \n");
printMatrix(K, N, B, K);
// ... your code here: Use dgemm_nn to compute C=A*B
// dgemm_nn( ??? );
printf("C = \n");
printMatrix(M, N, C, M);
return 0;
}