====
GEMM                                                                    [TOC]
====

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:

---LATEX------------------------------------------------------------------------
\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

---LATEX------------------------------------------------------------------------
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

---LATEX------------------------------------------------------------------------
\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).

---CODE(type=c)-----------------------------------------------------------------
#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;
}
--------------------------------------------------------------------------------


:navigate: __up__    -> doc:index
           __back__  -> doc:day08/page05