Content

GEMM

Algorithmus

Die Matrizen \(A, B\) und \(C\) werden in Blöcke unterteilt:

Man kann zu Berechnung von \(C\) folgenden Algorithmus verwenden:

Bemerkungen

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

#include <assert.h>
#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;
}