Content |
Einfache Cache-Optimierung des Matrix-Matrix Produktes
Ziel: Vergleich von 4 algorithmischen Varianten für die GEMM (General Matrix-Matrix Product)
\[C \leftarrow \beta\,C + \alpha\,A B\quad\text{mit}\quad A \in \mathbb{M}^{m \times k},\; B \in \mathbb{M}^{k \times n},\; C \in \mathbb{M}^{m \times n}.\]-
Variante 1 (dgemm_var1): Schul-Methode Zeile-mal-Spalte
-
Variante 2 (dgemm_var2): Alle Matrizen zeilenweise durchlaufen
-
Variante 3 (dgemm_var3): Alle Matrizen spaltenweise durchlaufen
-
Variante 4 (dgemm_var4): Blockweise, gepufferte Multiplikation.
Signatur für GEMM
Exemplarisch für die erste Variante soll die Signatur lauten:
void dgemm_var1(long m, long n, long k, double alpha, const double *A, long incRowA, long incColA, const double *B, long incRowB, long incColB, double beta, double *C, long incRowC, long incColC);
Aufgaben
-
Schreibt für die Varianten 1 bis 3 Algorithmen auf und implementiert diese. Für den Benchmark steht unten wieder eine Vorlage bereit. Diese soll wie folgt ergänzt werden:
-
Die Ergebnisse sollen auf Korrektheit verglichen werden:
-
Variante 1 berechnet \(C_1 \leftarrow \beta\,C_1 + \alpha\,A B\)
-
Variante 2 berechnet \(C_2 \leftarrow \beta\,C_2 + \alpha\,A B\)
-
Variante 3 berechnet \(C_3 \leftarrow \beta\,C_3 + \alpha\,A B\)
Anschließend sollen dann die Differenzen
\[\delta_2 = \| C_2 - C_1 \|_1\quad\text{mit}\quad\delta_3 = \| C_3 - C_1 \|_1\]berechnet und ausgegeben werden
-
-
Ausgabe von Dimensionen, Zeiten, MFLOPS, Differenzen in Tabellenform.
-
-
Erstellt Benchmarks für den Fall von zeilen- oder spaltenweise gespeicherten Matrizen. Dies ist in der Vorlage durch das Macro COLMAJOR steuerbar.
Einfacher Block-Algorithmus für GEMM
Entwickelt einen Algorithmus, der die Multiplikation blockweise durchführt.
Die Grundidee für den Algorithmus ist, dass die Matrizen zunächst quasi-äquidistant partitioniert betrachtet werden:
-
Matrix \(A\) bezüglich \(M_C\) und \(K_C\) in
\[A = \left(\begin{array}{c|c|c} \mathbf{A}_{1,1} & \dots & \mathbf{A}_{1,k_b} \\ \hline \vdots & & \vdots \\ \hline \mathbf{A}_{m_b,1} & \dots & \mathbf{A}_{m_b,k_b} \\ \end{array}\right)\quad\text{mit}\quad m_b = \lceil m/M_c \rceil,\; k_b = \lceil k/K_c \rceil\]und
\[\mathbf{A}_{i,j} \in \mathbb{M}^{M_C \times K_C}\quad 1 \leq i < m_b\; 1 \leq j < k_b.\] -
Matrix \(A\) bezüglich \(M_C\) und \(K_C\) in
\[B = \left(\begin{array}{c|c|c} \mathbf{B}_{1,1} & \dots & \mathbf{B}_{1, n_b} \\ \hline \vdots & & \vdots \\ \hline \mathbf{B}_{k_b,1} & \dots & \mathbf{B}_{k_b,n_b} \\ \end{array}\right)\quad\text{mit}\quad k_b = \lceil k/K_c \rceil,\; n_b = \lceil n/N_c \rceil.\]und
\[\mathbf{B}_{i,j} \in \mathbb{M}^{M_C \times K_C}\quad1 \leq i < k_b\;1 \leq j < n_b\;\] -
Matrix \(C\) bezüglich \(M_C\) und \(M_C\) in
\[B = \left(\begin{array}{c|c|c} \mathbf{C}_{1,1} & \dots & \mathbf{C}_{1, n_b} \\ \hline \vdots & & \vdots \\ \hline \mathbf{C}_{m_b,1} & \dots & \mathbf{C}_{m_b,n_b} \\ \end{array}\right)\quad\text{mit}\quad m_b = \lceil m/M_c \rceil,\; n_b = \lceil n/N_c \rceil.\]und
\[\mathbf{B}_{i,j} \in \mathbb{M}^{M_C \times K_C}\quad1 \leq i < m_b\;1 \leq j < n_b\;\]
Statt einzelner Elemente muliplizieren wir die Blöcke, nachdem diese zuerst zeilenweise oder spaltenweise in einen Puffer kopiert wurden:
-
\(\text{If}\; \beta \neq 1\)
-
\(C \leftarrow \beta\, C\)
-
-
\(\text{For}\; i=1,\dots,m_b\)
-
\(\text{For}\; j=1,\dots,n_b\)
-
\(\text{For}\; l=1,\dots,k_b\)
-
\(\overline{A} \leftarrow \mathbf{A}_{i,l}\)
-
\(\overline{B} \leftarrow \mathbf{B}_{l,j}\)
-
\(\overline{C} \leftarrow \alpha\,\overline{A} \, \overline{B}\)
-
\(\mathbf{C}_{i,j} \leftarrow \mathbf{C}_{i,j} + \overline{C}\)
-
-
-
Aufgabe
-
Ist die Schleifenanordnung im obigen Block-Algorithmus korrekt/ideal?
-
Verbessert den Algorithmus und implementiert diesen:
-
Für die Blockgrößen \(M_C\), \(N_C\), \(K_C\) verwenden wir Macros. Zum Beispiel:
#define M_C 256 #define K_C 256 #define N_C 1024
-
Zum Puffern der Blöcke können wieder globale Arrays benutzt werde. Zum Beispiel:
double A_[M_C*K_C]; double B_[K_C*N_C]; double C_[M_C*N_C];
-
Zum Kopieren und Addieren stehen in der Vorlage die Funktionen dgecopy und dgeaxpy bereit.
-
Vorlage
#include <math.h> #include <stdio.h> #include <stdlib.h> #include <sys/times.h> #include <unistd.h> #ifndef COLMAJOR #define COLMAJOR 0 #else #undef COLMAJOR #define COLMAJOR 1 #endif #ifndef MAXDIM_M #define MAXDIM_M 4000 #endif #ifndef MAXDIM_N #define MAXDIM_N 4000 #endif #ifndef MAXDIM_K #define MAXDIM_K 4000 #endif #ifndef MIN_M #define MIN_M 100 #endif #ifndef MIN_N #define MIN_N 100 #endif #ifndef MIN_K #define MIN_K 100 #endif #ifndef MAX_M #define MAX_M 4000 #endif #ifndef MAX_N #define MAX_N 4000 #endif #ifndef MAX_K #define MAX_K 4000 #endif #ifndef INC_M #define INC_M 100 #endif #ifndef INC_N #define INC_N 100 #endif #ifndef INC_K #define INC_K 100 #endif #ifndef ALPHA #define ALPHA 1.0 #endif #ifndef BETA #define BETA 0.0 #endif double A[MAXDIM_M*MAXDIM_K]; double B[MAXDIM_K*MAXDIM_N]; double C1[MAXDIM_M*MAXDIM_N]; double C2[MAXDIM_M*MAXDIM_N]; double C3[MAXDIM_M*MAXDIM_N]; double C4[MAXDIM_M*MAXDIM_N]; double walltime() { struct tms ts; static double ClockTick=0.0; if (ClockTick==0.0) { ClockTick = 1.0 / ((double) sysconf(_SC_CLK_TCK)); } return ((double) times(&ts)) * ClockTick; } void initMatrix(long m, long n, double *A, long incRowA, long incColA) { int i, j; for (j=0; j<n; ++j) { for (i=0; i<m; ++i) { A[i*incRowA+j*incColA] = ((double)rand() - RAND_MAX/2)*200/RAND_MAX; } } } void dgecopy(long m, long n, const double *A, long incRowA, long incColA, double *B, long incRowB, long incColB) { int i, j; for (j=0; j<n; ++j) { for (i=0; i<m; ++i) { B[i*incRowB+j*incColB] = A[i*incRowA+j*incColA]; } } } void dgeaxpy(long m, long n, double alpha, const double *X, long incRowX, long incColX, double *Y, long incRowY, long incColY) { long i, j; for (i=0; i<m; ++i) { for (j=0; j<n; ++j) { Y[i*incRowY+j*incColY] += alpha*X[i*incRowX+j*incColX]; } } } double asumDiffMatrix(long m, long n, const double *A, long incRowA, long incColA, double *B, long incRowB, long incColB) { int i, j; double asum = 0; for (j=0; j<n; ++j) { for (i=0; i<m; ++i) { asum += fabs(B[i*incRowB+j*incColB] - A[i*incRowA+j*incColA]); } } return asum; } //------------------------------------------------------------------------------ // 1. GEMM Variante: Schulmethode //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ // 2. GEMM Variante: Zeilenweise //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ // 3. GEMM Variante: Spaltenweise //------------------------------------------------------------------------------ //------------------------------------------------------------------------------ // 4. GEMM Variante: Bisschen mit Puffer //------------------------------------------------------------------------------ int main() { initMatrix(MAXDIM_M, MAXDIM_K, A, 1, MAXDIM_M); initMatrix(MAXDIM_K, MAXDIM_N, B, 1, MAXDIM_K); initMatrix(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M); dgecopy(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M, C2, 1, MAXDIM_M); dgecopy(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M, C3, 1, MAXDIM_M); dgecopy(MAXDIM_M, MAXDIM_N, C1, 1, MAXDIM_M, C4, 1, MAXDIM_M); long m, n, k; // Header-Zeile für die Ausgabe for (m = MIN_M, n = MIN_N, k = MIN_K; m <=MAX_M && n <= MAX_N && k <= MAX_K; m += INC_M, n += INC_N, k += INC_K) { double t0, t1, t2, t3, t4; long incRowA = (COLMAJOR) ? 1 : k; long incColA = (COLMAJOR) ? m : 1; long incRowB = (COLMAJOR) ? 1 : n; long incColB = (COLMAJOR) ? k : 1; long incRowC = (COLMAJOR) ? 1 : n; long incColC = (COLMAJOR) ? m : 1; t0 = walltime(); // GEMM Variante 1 hier aufrufen t1 = walltime() - t0; // Weitere Varianten aufrufen .... // Audgabe: Dimensionen, Zeit, MFLOPS, ... } }