=======================================
Cache-Optimiertes Matrix-Matrix Produkt                                                             [TOC]
=======================================

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:

  ---- CODE (type=c) ---------------------------------------------------------
  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:

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

  ---- CODE (type=c) ---------------------------------------------------------
  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:

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

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

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

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

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


:navigate: up    -> doc:index
           next  -> doc:session07/page06
           back  -> doc:session07/page04