Cache-Optimiertes Matrix-Matrix Produkt

Wir werden folgende Block-Faktoren verwenden:

Mathe-Notation	Makro-Name
\(M_c\)	DGEMM_MC
\(N_c\)	DGEMM_NC
\(K_c\)	DGEMM_KC
\(M_r\)	DGEMM_MR
\(N_r\)	DGEMM_NR

In level3.c werden wir drei Varianten für die GEMM-Operation und zwei Pack-Funktionen hinzufügen. Wir beschreiben zunächst nur die GEMM-Operation etwas genauer:

Die Funktion dgemm berechnte die Operation \(C \leftarrow \beta C + \alpha A B\) für eine \(m \times n\) Matrix \(C\), eine \(m \times k\) Matrix \(A\) und eine \(k \times n\) Matrix \(B\). Dabei werden die Matrizen \(A\), \(B\) und \(C\) jeweils quasi-äquidistant partitioniert:
- \(C\) bezüglich \(M_c\) und \(N_c\) in Blöcke \(C_{i,j}\),
- \(A\) bezüglich \(M_c\) und \(K_c\) in Blöcke \(A_{i,\ell}\) und
- \(B\) bezüglich \(K_c\) und \(N_c\) in Blöcke \(B_{\ell,j}\).
Mit den Hilfsfunktion dpack_A und dpack_B werden die dabei entstehenden Blöcke \(A_{i,\ell}\) und \(B_{\ell, j}\) für die weitere Verwendung jeweils in Puffer \(\overline{A}\) und \(\overline{B}\) kopiert. Für die gepackten Matrix-Blöcke gilt:
- \(\overline{A}\) ist ein \(m_c \times k_c\) Block mit \(m_c \leq M_C\) und \(k_c \leq K_C\).
- \(\overline{B}\) ist ein \(k_c \times m_c\) Block mit \(k_c \leq K_C\) und \(m_c \leq M_C\).
Nach der \(j\ell{}i\)-Block-Variante wird die Operation \(C \leftarrow \beta C + \alpha A B\) durchgeführt gemäß:
- for \(j=0,\dots\)
  - for \(\ell=0, \dots \)
    - \(\overline{B} \leftarrow B_{j,\ell}\)
    - for \(i=0, \dots\)
      - \(\overline{A} \leftarrow A_{\ell,i}\)
      - \(C_{i,j} \leftarrow \beta\,C_{i,j} + \alpha\, \overline{A} \overline{B}\) (dgemm_macro)
Als Signatur vorgegeben ist:
```
void
dgemm(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)
```
Die Funktion gemm_macro berechnet wiederum \(C \leftarrow \beta C + \alpha A B\) für eine \(m \times n\) Matrix \(C\), eine \(m \times k\) Matrix \(A\) und eine \(k \times n\) Matrix \(B\). Zu beachten ist:
- Die Matrizen \(A\) und \(B\) sind in einem speziell gepackten Format gespeichert (dazu unten mehr).
- Die Matrix \(C\) ist nach wie vor im Full-Storage-Format gespeichert.
- Es gilt \(m \leq M_c\), \(n \leq N_c\) und insbesondere \(k \leq K_C\).
Die Matrizen \(A\), \(B\) und \(C\) werden jeweils quasi-äquidistant partitioniert:
- \(C\) bezüglich \(M_r\) und \(N_r\) in Blöcke \(C_{i,j}\),
- \(A\) bezüglich \(M_r\) (und \(k\)) in horizontale Paneele \(A_{i}\) und
- \(B\) bezüglich (\(k\) und) \(N_r\) in vertikale Paneele \(B_{j}\).
Die Operation in Block-Form liest sich also als:
\[\left(\begin{array}{c|c|c} C_{0,0} & C_{0,N_r} & \dots \\ \hline C_{M_r,0} & C_{M_r,N_r} & \dots \\ \hline \vdots & \vdots & \ddots \\ \end{array}\right)\leftarrow\left(\begin{array}{c|c|c} C_{0,0} & C_{0,N_r} & \dots \\ \hline C_{M_r,0} & C_{M_r,N_r} & \dots \\ \hline \vdots & \vdots & \ddots \\ \end{array}\right)+\alpha\left(\begin{array}{c} A_{0} \\ \hline A_{M_r} \\ \hline \vdots \end{array}\right)\left(\begin{array}{c|c|c} B_{0} & B_{N_r} & \dots \end{array}\right)\]
Das Produkt hat somit die Form eines dyadischen Produktes, das wir berechnen mit
- for \(j=0,\dots\)
  - for \(i=0, \dots\)
    - \(C_{i,j} \leftarrow \beta\,C_{i,j} + \alpha\, A_i B_j\) (dgemm_micro)
Als Signatur vorgegeben ist:
```
void
dgemm_macro(int m, int n, int k,
            double alpha,
            const double *A,
            const double *B,
            double beta,
            double *C, int incRowC, int incColC)
```
Auch die Funktion gemm_micro berechnet \(C \leftarrow \beta C + \alpha A B\). Hier gilt für die Dimensionen:
- \(C\) ist eine \(M_r \times N_r\) Matrix im Full-Storage-Format. Die Dimension ist also fest vorgegeben.
- \(A\) ist eine \(M_r \times k\) Matrix und \(B\) eine \(k \times N_r\) Matrix.
- Für die Dimension \(k\) gilt: \(k \leq K_c\).
Für die Operation wird \(A\) spaltenweise und \(B\) zeilenweise betrachtet, so dass sich die Operation liest als:
\[C \leftarrow \beta C +\alpha\left(\begin{array}{c|c|c} \vec{a}_0 & \vec{a}_1 & \dots \end{array}\right)\left(\begin{array}{c} \vec{b}_0^T \\ \hline \vec{b}_1^T \\ \hline \vdots \end{array}\right)=\beta C + \alpha\left(\vec{a}_0 \vec{b}_0^T + \vec{a}_1 \vec{b}_1^T + \dots\right)\]
Als Signatur vorgegeben ist:
```
void
dgemm_micro(int k,
            double alpha,
            const double *A,
            const double *B,
            double beta,
            double *C, int incRowC, int incColC)
```

Packen von Matrix-Blöcken aus \(A\)

Sei \(A\) ein Matrix-Block der Dimension \(m_c \times k_c\) mit \(m_c \leq M_c\) und \(k_c \leq K_c\). Weiter sei \(\overline{A}\) ein Puffer mit \(M_c \times K_c\) Elementen, die wir mit \(p_{0}, p_{1}, \dots, p_{M_c \cdot K_c-1}\) bezeichnen.

Den Block \(A\) betrachten wir auf unterschiedliche Weise:

Bei der komponentenweisen Betrachtung verwenden wir die Notation
\[A =\left(\begin{array}{cccc} a_{0,0} & a_{0,1} & \dots & a_{0,k_c-1} \\ \vdots & \vdots & & \vdots \\ a_{m_c-1,0} & a_{m_c-1,1} & \dots & a_{m_c-1,k_c-1} \\ \end{array}\right)\]
Betrachten wir die vertikalen Paneelen, dann verwenden wir die Notation
\[A =\left(\begin{array}{c} A_{0} \\ \vdots \\ A_{m_p-1,0} \\ \end{array}\right)\quad\text{mit}\quad m_p = \left\lceil \frac{m_c}{M_r} \right\rceil\]

Der Puffer soll nach folgendem Algorithmus mit \(M_r \cdot k_c \cdot m_p\) Elementen gefüllt werden:

for \(j=0, \dots, k_c-1\)
- for \(\ell = 0, \dots, m_p-1\)
  - for \(i_0=0, \dots, M_r-1\)
    - \(i = \ell \cdot M_r + i_0\)
    - \(\nu = \ell \cdot M_r \cdot k + j\cdot M_r + i_0\)
    - if i<m_c
      - \(p_\nu \leftarrow a_{i,j}\)
    - else
      - \(p_\nu \leftarrow 0\)

In der nachfolgenden Testvorlage wird eine \(10\times 12\) Matrix \(A\) bezüglich \(M_c = 8\) und \(K_c = 10\) in vier Blöcke partitioniert. Diese sollen bezüglich \(M_r = 4\) wiederum in horizonatale Paneele unterteilt werden:

Implementiert die Pack-Funktion
In der Vorlage wird nur das Packen des Blockes links oben getestet. Ergänzt die weiteren Test.

Das Programm ist stand-alone benötigt also nicht das Projekt der letzten Session.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <sys/times.h>
#include <unistd.h>

//-- setup and print matrices --------------------------------------------------

void
initGeMatrix(int m, int n, double *A, int incRowA, int incColA)
{
    int i, j;

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            A[i*incRowA+j*incColA] = i*n + j + 1;
        }
    }
}

void
printGeMatrix(int m, int n, const double *A, int incRowA, int incColA)
{
    int i, j;

    for (i=0; i<m; ++i) {
        for (j=0; j<n; ++j) {
            printf("%10.4lf ", A[i*incRowA+j*incColA]);
        }
        printf("\n");
    }
    printf("\n\n");
}

//------------------------------------------------------------------------------

#ifndef DGEMM_MC
#define DGEMM_MC 8
#endif

#ifndef DGEMM_NC
#define DGEMM_NC 9
#endif

#ifndef DGEMM_KC
#define DGEMM_KC 10
#endif

#ifndef DGEMM_MR
#define DGEMM_MR 4
#endif

#ifndef DGEMM_NR
#define DGEMM_NR 3
#endif



//------------------------------------------------------------------------------

void
dpack_A(int mc, int kc, const double *A, int incRowA, int incColA, double *p)
{
    int i, i0, j, l, nu;
    int mp = (m+DGEMM_MR-1) / DGEMM_MR;

    // ... your code here ...
}

//------------------------------------------------------------------------------

#ifndef DIM_M
#define DIM_M 10
#endif

#ifndef DIM_K
#define DIM_K 12
#endif

#ifndef ROWMAJOR
#define INCROW_A 1
#define INCCOL_A DIM_M
#else
#define INCROW_A DIM_K
#define INCCOL_A 1
#endif

#define MIN(X,Y) (X) < (Y) ? (X) : (Y)

double A[DIM_M*DIM_K];
double A_buffer[DGEMM_MC*DGEMM_KC];

int
main()
{
    initGeMatrix(DIM_M, DIM_K, A, INCROW_A, INCCOL_A);

    printf("A=\n");
    printGeMatrix(DIM_M, DIM_K, A, INCROW_A, INCCOL_A);

    printf("Nach packen von Block links oben\n");
    dpack_A(MIN(DIM_M, DGEMM_MC), MIN(DIM_K, DGEMM_KC),
            A, INCROW_A, INCCOL_A, A_buffer);

    printf("A_buffer=\n");
    printGeMatrix(1, DGEMM_MC*DGEMM_KC, A_buffer, 0, 1);

    printf("Nach packen von Block rechts oben\n");
    // ... your code here ...

    printf("A_buffer=\n");
    printGeMatrix(1, DGEMM_MC*DGEMM_KC, A_buffer, 0, 1);

    printf("Nach packen von Block links unten\n");
    // ... your code here ...

    printf("A_buffer=\n");
    printGeMatrix(1, DGEMM_MC*DGEMM_KC, A_buffer, 0, 1);

    printf("Nach packen von Block rechts unten\n");
    // ... your code here ...

    printf("A_buffer=\n");
    printGeMatrix(1, DGEMM_MC*DGEMM_KC, A_buffer, 0, 1);

    return 0;
}

Packen von Matrix-Blöcken aus \(B\)

Analog soll ein Block der Dimension \(k_c \times n_c\) mit \(k_c \leq K_c\) und \(n_c \leq N_c\) bezüglich \(N_r\) in vertikale Paneele partioniert und gepackt werden (inklusive Zero-Padding):

Leitet dazu wieder einen Algorithmus zum Packen her.
Fügt einen entsprechenden Test in obige Vorlage hinzu.