====================================
GEMM für unterschiedliche Datentypen                                        [TOC]
====================================

Die Referenz-Implementierung der GEMM Operation zu verallgemeinern ist relativ
einfach.  Dies geschieht analog zu den vorigen Beispielen `gecopy`, `geaxpy`
und `gescal`.  Lediglich die Anzahl der Template-Parameter erhöht sich:

---- CODE (type=cc) ------------------------------------------------------------
namespace refblas {

template <typename Index, typename Alpha, typename TA, typename TB,
          typename Beta, typename TC>
void
gemm(Index m, Index n, Index k,
     Alpha alpha,
     const TA *A, Index incRowA, Index incColA,
     const TB *B, Index incRowB, Index incColB,
     Beta beta,
     TC *C, Index incRowC, Index incColC)
{
    hpc::ulmblas::gescal(m, n, beta, C, incRowC, incColC);
    for (Index j=0; j<n; ++j) {
        for (Index l=0; l<k; ++l) {
            for (Index i=0; i<m; ++i) {
                C[i*incRowC+j*incColC] += alpha*A[i*incRowA+l*incColA]
                                               *B[l*incRowB+j*incColB];
            }
        }
    }
}
--------------------------------------------------------------------------------

Etwas anspruchsvoller ist es, die geblockte GEMM Operation aufzubohren ohne die
Effizienz zu verlieren.  Insbesondere die optimierten Assembler Micro-Kernel
setzen voraus, dass die Elemente von Paneelen aus `A` und `B` den gleichen Typ
besitzen.

Andererseits ist es doch recht einfach, denn die Matrix-Multiplikation wird
blockweise durchgeführt.  Und bevor zwei Blöcke miteinander multipliziert werden,
müssen diese in Puffer gepackt werden.  Den Typ der beiden Puffer legen wir als
den aus __Session 10, Page 6__ bekannten `common_type` für `A`, `B` und `alpha`
fest.  Wie bei `gecopy` wird beim Packen also eventuell ein upcast
durchgeführt, wenn sich der Typ der Matrix-Elemente vom Typ des Puffer
unterscheiden.  Für die gepackten Blöcke kann also ein Micro-Kernel verwendet
werden, der vorraussetzt, dass Paneele homogene Datentypen enthalten.

Zu beachten ist allerdings noch, dass sich die Elemente von `C` noch im Typ
unterscheiden können.  Auf diesen Fall gehen wir weiter unten ein.  Wir
betrachten die Anpassung des geblockten GEMM Algorithmus in folgender
Reihenfolge:

- Festlegung von Blockgrößen durch Macros
- Festlegung von Blockgrößen durch Traits
- Funktionen zum Packen: `hpc::ulmblas::pack_A` und `hpc::ulmblas::pack_B`
- Macro-Kernel: `hpc::ulmblas::mgemm`
- Micro-Kernel: `hpc::ulmblas::ugemm`
- Frame-Algorithmus: `hpc::ulmblas::gemm`

Diese einzelnen Komponenten pflegen wir dabei in unsere neue
Verzeichnisstruktur ein.

:links: Session 10, Page 6 -> doc:session10/page06

Festlegung von Blockgrößen durch Macros
=======================================
Wir verwenden Macros, um für Blockgößen Default-Werte festzulegen, die wir aber
beim Übersetzen eventuell abändern können.  Dies geschieht nach dem Muster

---- CODE (type=cc) ------------------------------------------------------------
//
// double precision(double)
//
#ifndef D_BLOCKSIZE_MR
#define D_BLOCKSIZE_MR 4
#endif

#ifndef D_BLOCKSIZE_NR
#define D_BLOCKSIZE_NR 4
#endif

#ifndef D_BLOCKSIZE_MC
#define D_BLOCKSIZE_MC 256
#endif

#ifndef D_BLOCKSIZE_KC
#define D_BLOCKSIZE_KC 256
#endif

#ifndef D_BLOCKSIZE_NC
#define D_BLOCKSIZE_NC 4096
#endif
--------------------------------------------------------------------------------

Dabei steht `D` für `double`.  Analog benutzen wir gemäß der BLAS Konvention
Macros mit Prefix `S` für `float`, `C` für `std::complex<float>` und `Z` für
`std::complex<double>`.  Für andere Datentypen legen wir gemeinsame
Default-Werte mit Prefix `DEFAULT` fest.

Aufgabe
-------
Übernehmt das fertige __config.h__ und kopiert es in `hpc/ulmblas`.

:links: config.h -> file:session12/page01/hpc/ulmblas/config.h


Festlegung von Blockgrößen durch Traits
=======================================
Die oben definierten Macro verwenden wir wieder indirekt mit Hilfe einer
Traits-Klasse.  Damit kann die Blockgröße in Abhängigkeit eines
Template-Parameter zur Compile-Zeit festgelegt werden:

---- CODE (type=cc) ------------------------------------------------------------
template <typename Index, typename T>
void
foo(Index m, Index k, T *A)
{
    Index MR = hpc::ulmblas::BlockSize<T>::MR;
    /*
       ...
    */
}
--------------------------------------------------------------------------------


Aufgabe
-------
Die Traits-Klasse `hpc::ulmblas::BlockSize` packen wir in
`hpc/ulmblas/blocksize.h`.  Dabei sollt ihr folgende Vorlage um
Spezialisierungen für `float`, `double`, `std::complex<float>` und
`std::complex<double>` ergänzen:

---- CODE (type=cc) ------------------------------------------------------------
#ifndef HPC_ULMBLAS_BLOCKSIZE_H
#define HPC_ULMBLAS_BLOCKSIZE_H 1

#include <complex>
#include <hpc/ulmblas/config.h>

namespace hpc { namespace ulmblas {

template <typename T>
struct BlockSize
{
    static const int MC = DEFAULT_BLOCKSIZE_MC;
    static const int KC = DEFAULT_BLOCKSIZE_KC;
    static const int NC = DEFAULT_BLOCKSIZE_NC;
    static const int MR = DEFAULT_BLOCKSIZE_MR;
    static const int NR = DEFAULT_BLOCKSIZE_NR;


    static_assert(MC>0 && KC>0 && NC>0 && MR>0 && NR>0, "Invalid block size.");
    static_assert(MC % MR == 0, "MC must be a multiple of MR.");
    static_assert(NC % NR == 0, "NC must be a multiple of NR.");
};

/*

   Special cases of BlockSize for float, double, ...

 */

} } // namespace ulmblas, hpc

#endif // HPC_ULMBLAS_BLOCKSIZE_H
--------------------------------------------------------------------------------


Funktionen zum Packen von Blöcken
=================================
Als Vorlage zum Packen von Blöcken dient:

---- CODE (type=cc) ------------------------------------------------------------
template <typename Index, typename T>
void
pack_A(Index mc, Index kc,
       const T *A, Index incRowA, Index incColA,
       T *p)
{
    Index MR = BlockSize<T>::MR;
    Index mp = (mc+MR-1) / MR;

    for (Index j=0; j<kc; ++j) {
        for (Index l=0; l<mp; ++l) {
            for (Index i0=0; i0<MR; ++i0) {
                Index i  = l*MR + i0;
                Index nu = l*MR*kc + j*MR + i0;
                p[nu]   = (i<mc) ? A[i*incRowA+j*incColA]
                                 : T(0);
            }
        }
    }
}

template <typename Index, typename T>
void
pack_B(Index kc, Index nc,
       const TB *B, Index incRowB, Index incColB,
       T *p)
{
    Index NR = BlockSize<T>::NR;
    Index np = (nc+NR-1) / NR;

    for (Index l=0; l<np; ++l) {
        for (Index j0=0; j0<NR; ++j0) {
            for (Index i=0; i<kc; ++i) {
                Index j  = l*NR+j0;
                Index nu = l*NR*kc + i*NR + j0;
                p[nu]   = (j<nc) ? B[i*incRowB+j*incColB]
                                 : T(0);
            }
        }
    }
}
--------------------------------------------------------------------------------

Aufgaben
--------
- Die Funktionen `hpc::ulmblas::pack_A` und `hpc::ulmblas::pack_B` schieben
  wir gemeinsam in `hpc/ulmblas/pack.h`.
- Festlegen der Include-Guards und Namensräume sollte klar sein ;-)
- Welche `#include` Direktiven benötigt werden sollte ebenfalls nicht weiter
  schwer sein.
- Die Funktionen sollen so angepasst werden, dass sich Matrix-Elemente und
  Puffer-Elemente im Typ unterscheiden können.  Dies geschiet analog zu
  `gecopy` usw.

Etwas mehr Überblick und Weitsicht ist hier notwendig:

- Hängen die Blockgrößen `MR` und `NR` vom Typ der Matrix-Elemente
  oder vom Typ der Puffer-Elemente ab?

Macro-Kernel
============
Bei folgender Vorlage wurde ein Bug beseitigt.  Da der Puffer `C_` auf dem
Stack liegt kann diese Werte wie `NaN`und `inf` beihalten.  In diesem Fall
wird der Inhalt beim Aufruf des Micro-Kernel mit `beta=0` nicht mit Null
Initialisiert.  Mit `std::fill_n(C_, MR*NR, T(0))` wird der Puffer vor
jedem Aufruf des Micro-Kernel mit Null überschrieben.  Mehr zu `std::fill`
findet ihr __hier__.

:links: hier -> http://www.cplusplus.com/reference/algorithm/fill_n/

---- CODE (type=cc) ------------------------------------------------------------
template <typename Index, typename T>
void
mgemm(Index mc, Index nc, Index kc,
      T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    Index MR = BlockSize<T>::MR;
    Index NR = BlockSize<T>::NR;
    T C_[BlockSize<T>::MR*BlockSize<T>::NR];

    Index mp  = (mc+MR-1) / MR;
    Index np  = (nc+NR-1) / NR;
    Index mr_ = mc % MR;
    Index nr_ = nc % NR;

    for (Index j=0; j<np; ++j) {
        Index nr = (j!=np-1 || nr_==0) ? NR : nr_;

        for (Index i=0; i<mp; ++i) {
            Index mr = (i!=mp-1 || mr_==0) ? MR : mr_;

            if (mr==MR && nr==NR) {
                ugemm(kc, alpha,
                      &A[i*kc*MR], &B[j*kc*NR],
                      beta,
                      &C[i*MR*incRowC+j*NR*incColC],
                      incRowC, incColC);
            } else {
                // Fill buffer with zeros in case of infs and nans
                std::fill_n(C_, MR*NR, T(0));
                ugemm(kc, alpha,
                      &A[i*kc*MR], &B[j*kc*NR],
                      T(0),
                      C_, Index(1), MR);
                gescal(mr, nr, beta,
                       &C[i*MR*incRowC+j*NR*incColC],
                       incRowC, incColC);
                geaxpy(mr, nr, T(1), C_, Index(1), MR,
                       &C[i*MR*incRowC+j*NR*incColC],
                       incRowC, incColC);
            }
        }
    }
}
--------------------------------------------------------------------------------

Aufgabe
-------
Der Macro-Kernel erhält mit `A` und `B` gepackte Blöcke.  Das Packen geschieht
im Frame-Algorithmus und zwar so, dass `A` und `B` den gleichen Typ besitzen.
Allerdings kann sich dieser vom Typ für `beta` und `C` unterscheiden.  Ändert
die Implementierung so ab, dass dies möglich ist.  Überlegt euch dabei bereits
welche Signatur vom Micro-Kernel im Allgemeinen unterstützt werden muss.

Der Macro-Kernel soll in `hpc/ulmblas/mgemm.h` mit entsprechenden Include-Guards
und Namesräumen implementiert werden.


Micro-Kernel
============
Für den Fall, dass der Typ der Puffer-Elemente mit dem Typ für `beta` und
Puffer-Elementen identisch ist (homogene Typen), kann die bisherige
Implementierung verwendet werden:

---- CODE (type=cpp) -----------------------------------------------------------
template <typename Index, typename T>
void
ugemm(Index kc, T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    const Index MR = BlockSize<T>::MR;
    const Index NR = BlockSize<T>::NR;
    T P[BlockSize<T>::MR*BlockSize<T>::NR];

    for (Index l=0; l<MR*NR; ++l) {
        P[l] = 0;
    }
    for (Index l=0; l<kc; ++l) {
        for (Index j=0; j<NR; ++j) {
            for (Index i=0; i<MR; ++i) {
                P[i+j*MR] += A[i+l*MR]*B[l*NR+j];
            }
        }
    }
    for (Index j=0; j<NR; ++j) {
        for (Index i=0; i<MR; ++i) {
            C[i*incRowC+j*incColC] *= beta;
            C[i*incRowC+j*incColC] += alpha*P[i+j*MR];
        }
    }
}
--------------------------------------------------------------------------------

Für spezielle Element-Typen und Blöckgrößen `MR` und `NR` stellen wir zudem
optimierte Assembler Implementierungen zur Verfügung.

Fälle mit inhomogenen Datenen Typen (bei denen `beta` oder Elemente von `C`
einen anderen Typ als `alpha`, `A` oder  `B` besitzen) führen wir deshalb auf
den Fall von homogenen Typen zurück.  Diese Fälle werden von einer
verallgemeinerten `hpc::ulmblas::ugemm` Variante behandelt, die einen Puffer
`P` der Dimension `MR x NR` auf dem Stack anlegt.  Dieser ist vom gleichen Typ
wie die Elemente von `A` und `B`.  Die Operation wird dann wie folgt
durchgeführt:

- $P \leftarrow \alpha\,A B$  (Micro-Kernel mit homogenen Typen)
- $C \leftarrow \beta C$ (gescal)
- $C \leftarrow C + P$ (geaxpy)

Entscheidend ist dabei, dass hier eventuell ein optimierter Micro-Kernel
verwendet werden kann.

Aufgabe
-------
Die Micro-Kernel landen im Verzeichnis `hpc/ulmblas/kernels/`:

- `hpc/ulmblas/kernels/ugemm_ref.h` enthält obige Template-Funktion für den
  Fall homogener Typen.
- `hpc/ulmblas/kernels/ugemm_buf.h` enthält eine Verallgemeinerung für den
  Fall inhomogener Typen.  Dieser geht wie oben geschildert vor.  Zu beachten
  ist, dass der Puffer mit Null initialisiert wird ;-)
- Die Verwendung von Include-Guards und Namensräume ist ab jetzt
  selbstverständlich.

Die Liste an Micro-Kernel wird sich regelmäßig erweitern.  Deshalb soll es
möglich sein, dass mit einem einzigen `#include` alle verfügbaren Micro-Kernel
eingebunden werden.  Dazu dient folgendes Include-File:

:import: session12/page01/hpc/ulmblas/kernels/ugemm.h

Frame-Algorithmus
=================
Als Vorlage für den Frame-Algorithmus verwenden wir:

---- CODE (type=cpp) -----------------------------------------------------------
template <typename Index, typename T>
void
gemm(Index m, Index n, Index k,
     T alpha,
     const T *A, Index incRowA, Index incColA,
     const T *B, Index incRowB, Index incColB,
     T beta,
     T *C, Index incRowC, Index incColC)
{
    Index MC = BlockSize<T>::MC;
    Index NC = BlockSize<T>::NC;
    Index KC = BlockSize<T>::KC;

    Index mb = (m+MC-1) / MC;
    Index nb = (n+NC-1) / NC;
    Index kb = (k+KC-1) / KC;

    Index mc_ = m % MC;
    Index nc_ = n % NC;
    Index kc_ = k % KC;

    T *A_ = new T[MC*KC];
    T *B_ = new T[KC*NC];

    if (alpha==T(0) || k==0) {
        gescal(m, n, beta, C, incRowC, incColC);
        return;
    }

    for (Index j=0; j<nb; ++j) {
        Index nc = (j!=nb-1 || nc_==0) ? NC : nc_;

        for (Index l=0; l<kb; ++l) {
            Index   kc    = (l!=kb-1 || kc_==0) ? KC : kc_;
            T beta_ = (l==0) ? beta : T(1);

            pack_B(kc, nc,
                   &B[l*KC*incRowB+j*NC*incColB],
                   incRowB, incColB,
                   B_);

            for (Index i=0; i<mb; ++i) {
                Index mc = (i!=mb-1 || mc_==0) ? MC : mc_;

                pack_A(mc, kc,
                       &A[i*MC*incRowA+l*KC*incColA],
                       incRowA, incColA,
                       A_);

                mgemm(mc, nc, kc,
                      alpha, A_, B_, beta_,
                      &C[i*MC*incRowC+j*NC*incColC],
                      incRowC, incColC);
            }
        }
    }
    delete [] A_;
    delete [] B_;
}
--------------------------------------------------------------------------------

Aufgaben
--------
- Der Frame-Algorithmus soll nun wie die obige verallgemeinerte
  Referenz-Implementierung (ganz oben! Ja, genau, am Seiten Anfang) den
  allgemeinen Fall behandeln.  Damit ist bereits die Signatur festgelegt.
- Der Typ für den Puffer soll als `common_type` der Typen für `alpha`,  `A`
  und `B` festgelegt werden.  Dabei ist es schön, dass hier mehr als ein Typ
  verwendet werden kann:

  ---- CODE (type=cpp) ---------------------------------------------------------
  typedef typename std::common_type<float, int, double>::type  FooType;
  ------------------------------------------------------------------------------

  Wie dies realisiert wird ist allerdings ein eigenes Thema ...

- Das Ganze packen wir in `hpc/ulmblas/gemm.h`.


Test
====
Da wir die Matrix-Klassen noch noch angepasst haben und auch andere Hilfsmittel
noch übertragen werden müssen, sieht der low-level Benchmark wieder etwas
umständlicher aus:

:import: session12/page01/hpc/tests/gemm.cc

Aufgaben
--------
- Testen, ob alles tut.
- Das Matrix-Produkt für verschiedene Matrix-Element Typen sowie verschieden
  Typen für die Skalare `alpha` und `beta` testen: `float`, `double`,
  `std::complex<float>`, `std::complex<double>`, aber auch `int` oder `long`.
  Manche Kombinationen werden nicht übersetzen.  Hier sollte man sich jeweils
  Gedanken über den Grund machen.  Beispiel: Wo geht es schief, wenn man als
  Matrix-Element Typ `int` festlegt?



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