=================
GEMM Micro Kernel                                                        [TOC]
=================

Der Micro Kernel soll eine $m_r \times k_c$ Paneel $\hat{A}$ mit einer
$k_c \times n_r$ Paneel $\hat{B}$ multiplizieren.  Genauer gesagt soll die
Operation $C \leftarrow \beta C + \alpha \; \hat{A} \cdot \hat{B}$ durchgeführt
werden.  Die Werte $m_r$ und $n_r$ sind dabei zur Compile-Zeit bereits bekannt
der Wert $k_c$ erst zur Laufzeit.

Algorithmus
===========
Ein möglicher Algorithmus für die Berechnung von

---LATEX------------------------------------------------------------------------
\begin{eqnarray*}
    C &\leftarrow& \beta C + \alpha \; \hat{A} \cdot \hat{B}
      & = & \beta C + \alpha
        \left(
            \begin{array}{ccc}
            \vec{a}_1 & \dots & \vec{a}_{k_c}
            \end{array}
        \right)
        \cdot
        \left(
            \begin{array}{c}
            \vec{b}_{ 1}^T  \\
            \vdots          \\
            \vec{b}_{k_c}^T \\
            \end{array}
        \right)
      & = & \beta C + \alpha
            \left(
            \vec{a}_1 \cdot \vec{b}_{ 1}^T +
            \dots +
            \vec{a}_{k_c} \cdot \vec{b}_{k_c}^T
            \right)
\end{eqnarray*}
--------------------------------------------------------------------------------

Könnte sich an folgender Idee orientieren:
 - Lege einen lokalen Puffer $AB$ mit $m_r \cdot n_r$ Elemente an und
   initialisiere ihn mit Null.
 - Überschreibe den Puffer $AB$ mit
   $\vec{a}_1 \cdot \vec{b}_{ 1}^T+\dots + \vec{a}_{k_c} \cdot \vec{b}_{k_c}^T$
   (in diesem Schritt wird also richtig gerechnet!)
 - Falls $\beta=0$ ist initialisiere $C$ mit Nullen und sonst mit $\beta C$.  Im
   letzteren Fall hat man also $m_r \cdot n_r$ zusätzliche Multiplikationen.
 - Falls $\alpha=1$ ist überschreibe $C$ mit $C+AB$ sonst mit
   $C+\alpha \cdot AB$.  Im letzteren Fall hat man also $m_r \cdot n_r$
   zusätzliche Multiplikationen.

Weitere Bemerkung zur Implementierung
=====================================
Die entscheidende Micro-Micro Operation ist also die Berechnung von

---LATEX------------------------------------------------------------------------
\begin{eqnarray*}
\vec{a}_l \cdot \vec{b}_l^T &=&
\left(
      \begin{array}{c}
      a_1^{(l)} \\ \vdots \\ a_{m_r}^{(l)}
      \end{array}
\right)
\left(
      \begin{array}{ccc}
      b_1^{(l)} & \dots & b_{n_r}^{(l)}
      \end{array}
\right)\\[0.5cm]
&=&
\left(
      \begin{array}{c}
      a_1^{(l)} \\ \vdots \\ a_{m_r}^{(l)}
      \end{array}
\right)
b_1^{(l)}
+ \dots +
\left(
      \begin{array}{c}
      a_1^{(l)} \\ \vdots \\ a_{m_r}^{(l)}
      \end{array}
\right)
b_{n_r}^{(l)} \\[0.5cm]
&=&
 \vec{a}_l \cdot b_1^{(l)} + \dots + \vec{a}_l \cdot b_{n_r}^{(l)}
\end{eqnarray*}
--------------------------------------------------------------------------------

Wenn man $m_r$ richtig wählt, dann können die Elemente von $\vec{a}_l$ also
während der gesamten Micro-Micro Operation in den Registern bleiben!

Aufgabe
=======
Implementiert einen Micro-Kernel in C. Unten gibt es wieder ein Grundgerüst mit
einem Mini-Test im Hauptprogramm:
  - Was wird im Hauptprogramm genau berechnet bzw. getestet?
  - Ihr könnt mit Matlab nachrechnen ob euer Ergebnis stimmt.
  - Verwendet auch verschiedene Werte für $\alpha$ und $\beta$.

---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 _A[MC*KC];
double _B[KC*NC];
double _C[MR*NR];

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 ...
}

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);


    pack_A(MC, KC, A, 1, M, _A);
    pack_B(KC, NC, B, 1, K, _B);

    printf("Buffer: _A = \n");
    printMatrix(MR, KC*MC/MR, _A, MR);

    printf("Buffer: _B = \n");
    printMatrix(NR, KC*NC/NR, _B, NR);

    dgemm_micro_kernel(KC, 1.0, _A, _B, 0.0, _C, 1, MR);

    printf("Buffer: _C = \n");
    printMatrix(MR, NR, _C, MR);

    return 0;
}
--------------------------------------------------------------------------------

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