=====================================
Template-ifizierung: `gemm_blocked.h`
=====================================

Bei der Verallgemeinerung der geblockten GEMM-Implementierung ist folgendes
zu beachten:

- Die Blockgößen `MC`, `NC`, `KC`, `MR` und `NR` sind vom Element-Typ abhängig.
  Diese Konstanten mittels Macros festzulegen ist zwar prinzipiell möglich aber
  extrem unschön.
- Die Puffer sind ebenfalls vom Element-Typ abhängig.  Bisher wurden diese
  für Elemente vom Typ `double` hart-codiert auf dem Datensegment angelegt.
  Damit ist auch die Größe hart-codiert.  Zudem ist diese Umsetzung auch nicht
  Thread-Safe (dazu in späteren Vorlesungen mehr).

Wir verwenden folgende Techniken, um dies zu beheben:

- Die Blockgrößen werden durch Traits festgelegt.  Für wichtige Speziallfälle
  von Element-Typen wie `float`, `double`, `std::complex<float>` und
  `std::complex<double>` können Größen individuell festgelegt werden.  Bei
  anderen Element-Typen erhält man einen Default-Wert.
- Die Puffer werden im Frame-Algorithmus bei jedem Aufruf dynamisch erzeugt
  und am Ende wieder freigegeben.  Auf manchen Systemen wird dies allerdings
  zu signifikanten Laufzeiteinbußen führen.  Dies werden wir später beheben.


Erster Kontakt mit Traits
=========================
Kurz gesagt ermöglichen es die sogenannten _Traits_, dass zur Compilezeit
bekannte Typen und Konstanten auf andere Typen und Konstanten abgebildet
werden.

Typen auf Typen abbilden
------------------------
Betrachten wir folgendes Beispiel, das das Maximum von zwei Integer-Zahlen
zurückgibt:

---- CODE (type=cc) ------------------------------------------------------------
int
max(int a, int b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

Für Zahlen vom Typ `double` würde man dies analog programmieren:

---- CODE (type=cc) ------------------------------------------------------------
double
max(double a, double b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

Und schnell hat man die Idee, dass sich dies mit Templates wie folgt
verallgemeinern lässt:

---- CODE (type=cc) ------------------------------------------------------------
template <typename T>
T
max(T a, T b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

Allerdings setzt dies voraus, dass beide Argumente den gleichen Typ besitzen.
Das heisst der Funktionsaufruf

---- CODE (type=cc) ------------------------------------------------------------
double m = max(1, 2.0);
--------------------------------------------------------------------------------

erzeugt einen Compiler-Fehler, denn `1` hat den Typ `int` und `2.0` den Typ
`double`. Dies hat aber nichts mit den Templates zu tun.  Auch dann, wenn nur
überladene Funktionen benutzt werden müsste man dazu eine Variante

---- CODE (type=cc) ------------------------------------------------------------
double
max(int a, double b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

bereitstellen.  Und für einen Aufruf `max(2.0, 1)` eine weitere.  Deshalb wäre
es natürlich interessant, eine Template-Funktion `max` zu implementieren, bei
der die Argumente unterschiedliche Typen besitzen können.  Die Frage ist aber,
wie man hier den Rückgabe-Typ festlegen soll:

---- CODE (type=cc) ------------------------------------------------------------
template <typename T, typename S>
/* ????? */
max(T a, S b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

Sind die Paramter `T` und `S` vom Typ `int` oder `double`, dann sollte der
Rückgabe-Wert in Abhängigkeit von `T` und `S` wie folgt festgelegt werden:

   +--------+---------+--------------+
   | T      | S       | Rückgabe-Typ |
   +--------+---------+--------------+
   |`int`   |`int`    |`int`         |
   +--------+---------+--------------+
   |`double`|`int`    | `double`     |
   +--------+---------+--------------+
   |`int`   |`double` | `double`     |
   +--------+---------+--------------+
   |`double`|`double` | `double`     |
   +--------+---------+--------------+

Dies ist offensichtlich eine Verknüpfungstabelle, die `T` und `S` auf
den Rückgabe-Typen abbilden.  Diese soll durch eine Template-Klasse `Decl`
realisiert werden, die wir dann wie folgt verwenden können:

---- CODE (type=cc) ------------------------------------------------------------
template <typename T, typename S>
typename Decl<T,S>::Type
max(T a, S b)
{
    return (a>b) ? a : b;
}
--------------------------------------------------------------------------------

Erster Wurf
~~~~~~~~~~~
In einem ersten Wurf gehen wir völlig unelegant vor und codieren die
vollständige Tabelle mit C++-Mitteln.  Der Trick besteht darin, eine
Template-Klasse für alle aufgelisteten Fälle zu spezialieren und dabei
den Rückgabe-Typen durch einen `typedef` festzulegen:

---- CODE (file=session10/example4/version1/decl.h) ----------------------------
template <typename T, typename S>
struct Decl
{
};

template <>
struct Decl<int, int>
{
    typedef int Type;
};

template <>
struct Decl<double, int>
{
    typedef double Type;
};

template <>
struct Decl<int, double>
{
    typedef double Type;
};

template <>
struct Decl<double, double>
{
    typedef double Type;
};
--------------------------------------------------------------------------------

Um zu testen, ob diese Abbildung tatsächlich zur Compile-Zeit durchgeführt wird,
kann folgendes Test-Programm benutzt werden:

---- CODE (file=session10/example4/test_decl.cc) -------------------------------
#include <cstdio>
#include <decl.h>

template <typename T>
void
checkType(T)
{
    printf("Unknown type\n");
}

void
checkType(int)
{
    printf("int\n");
}

void
checkType(double)
{
    printf("double\n");
}

template <typename T, typename S>
typename Decl<T,S>::Type
max(T a, S b)
{
    return (a>b) ? a : b;
}

int
main()
{
    checkType(max(1,   1.0));
    checkType(max(1.0, 1  ));
    checkType(max(1,   1  ));
    checkType(max(1.0, 1.0));
}
--------------------------------------------------------------------------------

---- SHELL (path=session10/example4) -------------------------------------------
g++ -Wall -I version1 -o test_decl test_decl.cc
./test_decl
--------------------------------------------------------------------------------


Zweiter Wurf
~~~~~~~~~~~~
Zunächst stört uns, dass wir offensichtliche Fälle explizit codieren müssen
und das ganze Konstrukt nur für `int` und `double` funktioniert.
Sind beide Typen gleich, dann ist auch genau dies der Rückgabe-Typ.
Dies können wir durch eine teilweise Spezialisierung umsetzen:

---- CODE (file=session10/example4/version2/decl.h) ----------------------------
template <typename T, typename S>
struct Decl
{
};

// case where both types are the same
template <typename T>
struct Decl<T,T>
{
    typedef T Type;
};


template <>
struct Decl<double, int>
{
    typedef double Type;
};

template <>
struct Decl<int, double>
{
    typedef double Type;
};
--------------------------------------------------------------------------------

---- SHELL (path=session10/example4) -------------------------------------------
g++ -Wall -I version2 -o test_decl test_decl.cc
./test_decl
--------------------------------------------------------------------------------

Dritter Wurf
~~~~~~~~~~~~
Natürlich ist uns die Symmetrie der Tabelle nicht entgangen:  `Decl<T,S>` und
`Decl<S,T>` sollten immer den gleichen Typ definieren.  Auch das können wir
realisieren:

---- CODE (file=session10/example4/version3/decl.h) ----------------------------
template <typename T, typename S>
struct Decl
{
    typedef typename Decl<S,T>::Type Type;
};

// case where both types are the same
template <typename T>
struct Decl<T,T>
{
    typedef T Type;
};

template <>
struct Decl<int, double>
{
    typedef double Type;
};
--------------------------------------------------------------------------------

---- SHELL (path=session10/example4) -------------------------------------------
g++ -Wall -I version3 -o test_decl test_decl.cc
./test_decl
--------------------------------------------------------------------------------

Dadurch wurde der Programmieraufwand deutlich reduziert. Dennoch muss man diese
Tabelle für alle anderen in Frage kommenenden Typen (`long`, `float`, ...) und
Kombinationen mit viel Fleiß ergänzen.  Aber seit C++11 wurde dies von anderen
erledigt!

Ab C++11
~~~~~~~~
In `<type_traits>` sind eine Vielzahl an nützlichen Trait-Klassen definiert
(__Link__).  Eine davon ist `std::common_type<T,S>::type`.  Dies liefert den
Typ, in den sowohl Instanzen von `T` als auch `S` konvertiert werden können:

---- CODE (file=session10/example4/test_decl2.cc) -------------------------------
#include <cstdio>
#include <type_traits>

template <typename T>
void
checkType(T)
{
    printf("Unknown type\n");
}

void
checkType(int)
{
    printf("int\n");
}

void
checkType(double)
{
    printf("double\n");
}

template <typename T, typename S>
typename std::common_type<T,S>::type
max(T a, S b)
{
    return (a>b) ? a : b;
}

int
main()
{
    checkType(max(1,   1.0));
    checkType(max(1.0, 1  ));
    checkType(max(1,   1  ));
    checkType(max(1.0, 1.0));
}
--------------------------------------------------------------------------------

---- SHELL (path=session10/example4) -------------------------------------------
g++ -std=c++11 -Wall -I version1 -o test_decl2 test_decl2.cc
./test_decl2
--------------------------------------------------------------------------------

:links: Link -> http://en.cppreference.com/w/cpp/header/type_traits

Typen auf Konstanten abbilden
-----------------------------
Hier haben wir ein konkretes Beispiel.  Der Element-Typ einer Matrix soll
auf die Block-Dimensionen `MC`, `NC`, `KC`, `MR` und `NR` abgebildet werden.
Beispielsweise:

+-------------------------+---------+-------+---------+-------+-----------+
| Datentyp                |  `MC`   | `KC`  | `NC`    | `MR`  |   `NR`    |
+-------------------------+---------+-------+---------+-------+-----------+
| `float`                 |  256    | 512   | 4096    | 8     |  8        |
+-------------------------+---------+-------+---------+-------+-----------+
| `double`                |  256    | 256   | 4096    | 4     |  8        |
+-------------------------+---------+-------+---------+-------+-----------+
| `std::complex<float>`   |  256    | 256   | 4096    | 4     |  8        |
+-------------------------+---------+-------+---------+-------+-----------+
| `std::complex<double>`  |  256    | 128   | 4096    | 4     |  4        |
+-------------------------+---------+-------+---------+-------+-----------+
| Sonst                   |  64     | 64    | 256     | 2     |  2        |
+-------------------------+---------+-------+---------+-------+-----------+

Mit Traits werden wir folgendes ermöglichen:

---- CODE (type=cc) ------------------------------------------------------------
template <typename T>
void
foo(/* ...*/)
{
    int MC = BlockSize<T>::MC;

    /* ... */
}
--------------------------------------------------------------------------------

Dazu definieren wir eine Klasse `BlockSize` mit einem Template-Parameter für
den Element-Typ.  Dieser enthält in der allgemeinen Form die Default-Werte
als statische Konstanten:

---- CODE (type=cc) ------------------------------------------------------------
template <typename T>
struct BlockSize
{
    static const int MC = 64;
    static const int KC = 64;
    static const int NC = 256;
    static const int MR = 2;
    static const int NR = 2;
};
--------------------------------------------------------------------------------

Durch Spezialisierungen kann obige Tabelle codiert werden:

---- CODE (type=cc) ------------------------------------------------------------
template <>
struct BlockSize<float>
{
    static const int MC = 256;
    static const int KC = 512;
    static const int NC = 4096;
    static const int MR = 8;
    static const int NR = 8;
};

template <>
struct BlockSize<double>
{
    static const int MC = 256;
    static const int KC = 256;
    static const int NC = 4096;
    static const int MR = 4;
    static const int NR = 8;
};

/* ... */
--------------------------------------------------------------------------------

Damit haben wir vorerst alle Werkzeuge, um die geblockte GEMM-Implementierung
zu verallgemeinern.


Vorlage für `gemm_blocked.h`
============================

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

#include <complex>
#include <type_traits>
#include "ulmblas.h"

namespace blocked {

template <typename T>
struct BlockSize
{
    /* ... */
};

template <typename T, typename Index>
void
pack_A(Index mc, Index kc,
       const T *A, Index incRowA, Index incColA,
       T *p)
{
    /* ... */
}

template <typename T, typename Index>
void
pack_B(Index kc, Index nc,
       const T *B, Index incRowB, Index incColB,
       T *p)
{
    /* ... */
}

template <typename T, typename Index>
void
ugemm(Index kc, T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    /* ... */
}

template <typename T, typename Index>
void
mgemm(Index mc, Index nc, Index kc,
      T alpha,
      const T *A, const T *B,
      T beta,
      T *C, Index incRowC, Index incColC)
{
    /* ... */
}

template <typename T, typename Index>
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)
{
    /* ... */
}

} // namespace blocked

#endif // HPC_GEMM_BLOCKED_H
--------------------------------------------------------------------------------

:navigate: up    -> doc:index
           back  -> doc:session10/page05
           next  -> doc:session10/page07