Massiv parallele Prozessoren: GPUs#

Graphics Processing Units (GPUs) stellen sehr viele parallele Rechenwerke zur Verfügung. Heutzutage werden diese Rechenwerke für vielfältige Berechnungen außerhalb von Computergraphik verwendet, z.B. zum Trainieren von Neuronalen Netzen oder für wissenschaftliche Berechnungen. Man spricht von GPGPU: General Purpose computation on GPUs.

GPUs wurden als massive parallele Prozessoren entwickelt, da für die Darstellung einer 3D Graphik sehr viele jedoch unabhängige Berechnungen angestellt werden müssen, z.B.:

  • 3D-Objekte werden je nach Komplexität aus vielen tausend Dreiecken zusammengesetzt. Für die Darstellung auf einem Bildschirm muss jeder Punkt von jedem Dreieck in Abhängigkeit von der Position des Betrachters in eine 2D-Koordinate umgerechnet werden

  • Jedem Pixel eines Dreiecks muss in Abhängigkeit von der Oberfläche (Texture) und den Lichtverhältnissen ein Farbwert zugewiesen werden

Die Stärke von GPUs liegt also darin die gleiche Berechnung (3D zu 2D Koordinaten, Shading) parallel auf sehr vielen Elementen (Dreiecke, Pixel) durchzuführen. Wenn man vom Graphikkontext verallgemeinert, dann entspricht dies dem Konzept von Datenparallelismus.

Die ersten Anwendungen von GPGPU ab ca. 2002 mussten die Probleme noch in 3D-Graphikberechnungen umwandeln und über entsprechende Graphikframeworks ausführen. Mittlerweile gibt es Frameworks, die diesen Umweg unnötig machen und auch die GPU-Entwickler unterstützen Features, die nicht Graphiken notwendig sind aber im GPGPU Kontext nützlich. Die Technologien, Sprachen, Bibliotheken für GPGPU sind derzeit (2021) noch von Inkompatibilitäten, Vendor Lockins und stark unterschiedlichen Abstraktionsleveln geprägt. Am beliebtesten ist CUDA - eine C-ähnliche Programmiersprache mit reichhaltigem Tooling. Der Nachteil ist, dass CUDA nur mit den GPUs des Herstellers NVIDIA funktioniert. OpenCL ist ein offener Standard, der von den meisten GPUs und auch CPUs unterstützt wird. Die Verbreitung ist aber geringer als bei CUDA. Die online und in Büchern verfügbaren Dokumentationen und Codebeispiele sind auch weniger umfangreich. Apple als einer der Initiatoren von OpenCL zieht die Unterstützung zurück und forciert eigene Ansätze wie Metal. Viele der Konzepte sind gut übertragbar zwischen den Standards - die Bezeichnungen sind aber oft unterschiedlich.

Massiv mehr Rechenkerne#

Bei der Programmierung von CPUs finden wir aktuell typischerweise zwischen 4 und 64 Kernen vor. Die Überlegung ist wie wir ein Problem auf ebensoviele Threads aufteilen können, um alle Kerne auszulasten. Im GPU-Szenario haben wir sehr viele Threads, z.B.: eine Farbmaske kann unabhängig für jeden Pixel in einer 1920 x 1080 Pixel großen Szene in einem eigenen Thread berechnet werden - hier haben wir also ca. 2 Millionen Threads. Das heißt wir können sehr hohe Zahlen von Rechenkerne auslasten. Aktuelle Graphikkarten im Jahr 2021 besitzen im Bereich von 1.000 bis 10.000 Rechenkerne, wobei hier die Begriffe und Zählweisen sehr unterschiedlich sind:

  • Apple M1 mit 8 GPU Core - jeder Core entspricht aber 128 Recheneinheiten, also in Summe 1024

  • NVIDIA RTX 3090 mit ca. 10.000 CUDA Cores - es sind auch um Faktor 32 höhere Zahlen bzgl. der möglichen Threads zu lesen - diese Threads können zwar gleichzeitig in den Cores geladen sein aber nicht für alle Threads können parallel Berechnungen durchgeführt werden. Dies hängt u.a. von den verfügbaren Rechenwerken pro Core ab (z.B. FP32, FP64, Integer, Tensor)

Wie auch immer gezählt wird: es befinden sich deutlich mehr Rechenkerne auf einer GPU als auf einer CPU. Dies ist möglich aufgrund der Eigenschaften von Graphikberechnungen, die es ermöglichen einige Einschränkungen zu akzeptieren, die im Gegenzug kleinere und damit mehr Kerne erlauben:

  1. Es gibt mehr Threads als Kerne und wir müssen warten bis alle Threads fertig sind. Throughput/Bandbreite wird wichtiger als Latenz. Wir verzichten auf große Caches (Latenz veringern). Wenn ein Thread auf Daten aus dem Speicher wartet wird so lange ein anderer Thread ausgeführt

  2. Die einzelnen Berechnungen sind relativ klein und haben wenige Verzweigungen. Wir können auf umfangreiche Komponenten für Instruction Level Parallelism verzichten

  3. Die Berechnungen für unterschiedliche Datenelemente sind die gleichen. Mehrere Recheneinheiten können sich die Logik für Instruktionsdekodierung teilen (ähnlich zu SIMD).

Verzweigungen sind ungeachtet der letzten beiden Punkte möglich über Maskierung. Dabei werden zwar auf allen Rechenkernen alle Instruktionen ausgeführt aber je nach Maske/Ergebnis einer If-Abfrage werden die nicht benötigten Ergebnisse verworfen. Dies kann im Extremfall zu Ineffizienzen führen, z.B. wenn wir auf 8 Recheneinheiten mit gemeinsamer Instruktionspipeline das folgende (sinnlose) Programm laufen lassen:

int i = 0;
switch(thread_id) {
    case 0: i += 5; break;
    case 1: i += 3; break;
    case 2: i += 2; break;
    case 3: i += 9; break;
    case 4: i += 7; break;
    case 5: i += 8; break;
    case 6: i += 17; break;
    case 7: i += 1; break;
}

Alle 8 Recheneinheiten müssen alle 8 Additionen nacheinander ausführen wovon jeweils 7 verworfen werden.

INFO
Illustrationen zur Architektur von GPUs auf Folien 8 - 27 Beispiele in Lectures 7, CS267 der UC Berkeley

Architektur von GPUs#

Die genaue Architektur von GPUs unterscheidet sich je nach Hersteller und die Namen für die einzelnen Teile je nach Hersteller und Programmiermodell. Wir orientieren uns im Folgenden an OpenCL. OpenCL hat als Ziel nicht nur GPUs sondern allgemeiner “Compute Devices” zu unterstützen dazu zählen u.a. auch CPUs. Im Architekturmodell gibt es einen Host (“Computer”) mit mehreren Compute Devices. Diese können mit clinfo angezeigt werden, z.B.:

$ clinfo | grep "Device Name"
  Device Name                                     Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz
  Device Name                                     Intel(R) UHD Graphics 630
  Device Name                                     AMD Radeon Pro 5300M Compute Engine
  ...

Hier gibt es Compute Devices für den Prozessor, die integrierte Graphikkarte und die dedizierte Graphikkarte. Eine Visualierung des obersten Modells und weiterer Ebenen ist in der (OpenCL Spezifikation)[https://www.khronos.org/registry/OpenCL/specs/opencl-1.2.pdf] zu finden. Ein Vorteil von OpenCL ist, dass ein mal geschriebener Code auf allen 3 Devices läuft, obwohl alle unterschiedliche Befehlssätze hat. Dies ist vergleichbar mit der JVM, die Java Binaries (.class-Dateien) auf unterschiedlichen Prozessorarchitekturen ausführen kann.

Ein Compute Device enthält Compute Units (bei CUDA “Streaming multiprocessor (SM)”) und einen globalen Speicher, der von allen Compute Units genutzt werden kann. Jede Compute Unit enthält wiederum Processing Elements (bei CUDA: “CUDA Core”) und einen lokalen Speicher, der von allen Processing Elements genutzt werden kann. Jedes Processing Element enthält darüber hinaus noch einen private Speicher (auch als Register bezeichnet).

Zu Beachten ist, dass der globale Speicher des Compute Devices sich vom Hauptspeicher des Hosts unterscheidet und Daten explizit hin und her übertragen werden müssen.

Arbeit für die GPU: Kernel-Funktionen#

Das Grundprinzip um hohe Parallelität in der GPU-Programmierung zu erreichen: die Definition einer Kernel-Funktion, die für jeden Datenpunkt des zu lösenden Programms ausgeführt wird. Dann kann für jeden Datenpunkt eine Instanz des Kernels ausgeführt werden. Diese Instanzen heißen in OpenCL work-items und können parallel ausgeführt werden. Die Kernel-Funktion entspricht oftmals dem innersten Codeblock einer (Reihe von) Schleifen. Der Unterschied ist, dass Schleifendurchläufe eine sequenzielle Semantik haben und nur parallelisiert werden können, wenn eine Abhängigkeit zwischen den Durchläufen ausgeschlossen ist.

Z.B. Vektor-Addition mit For-Schleife:

void vector_add(const int n, const float *a, const float *b, float *c) 
{
    for(int i=0; i<n; i++)
        c[i] = a[i] + b[i];
}

OpenCL Kernel-Funktion:

__kernel void vector_add(__global const float *a, __global const float *b, __global float *c)
{
    int i = get_global_id(0);
    c[i] = a[i] + b[i];
}

Die Kernel-Funktionen werden in OpenCL C geschrieben - es gibt einige Abweichungen zu Standard C (siehe z.B. __kernel und __global) auf die wir später eingehen werden.

Ausführung von OpenCL-Programmen#

OpenCL ist sehr flexibel bezüglich der Anzahl und Art der Compute Devices und übersetzt Kernel-Funktionen zur Laufzeit für die verschiedenen Geräte. Der Preis dafür ist ein relativ hoher Aufwand, um alle benötigten Datenstrukturen aufzusetzen bevor der erste Kernel ausgeführt werden kann. Glücklicherweise ist der Code dafür für die meisten Programme immer wieder der gleiche. Das prinzipielle Vorgehen im einfachsten Fall ist:

  • Herausfinden welche Plattformen und zugehörigen Compute Devices es gibt und welche Eigenschaften diese haben

  • Selektieren von einem Compute Devices und erstellen eines Contexts mit diesem Device

  • Erstellen einer Command Queue über die Aufgaben an das Device übergeben werden können

Woher kommen diese Aufgaben? Aus Programm-Dateien: Sammlungen von Kernel-Funktionen. Diese können z.B. als Strings im Host-Programm kodiert sein oder in einer extra Datei liegen. Ein einfaches Vorgehen ist:

  • Einlesen einer .cl-Datei das das Programm enthält in einen C-String

  • Programm im entsprechenden Context erstellen und compilen

  • Für jede Kernel-Funktion im Programm einen cl_kernel erstellen

Nun zur Ausführung:

  • Erstellen von Buffern für die Ein- und Ausgaben des Kernels im Speicher des Compute Device

  • Ggfs. Kopieren von Daten vom Host zum Compute Device

  • Übergeben des Kernels an die Command Queue

  • Warten auf Ausführung

  • Ggfs. Kopieren von Daten vom Compute Device zum Host

INFO
Gemeinsame Besprechung des Beispiel/Skelett-Programms:
opencl-skeleton.c
opencl-program.cl
Ebenso ein simplistischer Benchmark: opencl-skeleton.c

Daten und Speicher#

Wie wir schon im vector_add-Beispiel gesehen haben, können wir einem Kernel Argumente übergeben. Dazu gibt es die Funktion clSetKernelArg(cl_kernel_id kernel, cl_uint index, size_t size, const void *value) Über den Index wird bestimmt an welchen Parameter das Argument übergeben werden soll, es wird beginnend bei 0 gezählt. Der value ist entweder ein Pointer auf einen primitiven Datentyp, z.B. ein Integer, oder auf ein cl_mem Speicherobjekt.

Speicherobjekte werden mit der Funktion clCreateBuffer(cl_context context, cl_mem_flags options, size_t size, void *host_ptr, cl_int *error) erstellt. Die Optionen beinhalten Angaben wie das Compute Device auf die Daten zugreifen kann: CL_MEM_READ_WRITE, CL_MEM_READ_ONLY, CL_MEM_WRITE_ONLY. Darüber hinaus gibt es noch Optionen wie OpenCL mit dem Speicherbereich auf dem Host (host_ptr) umgeht. Sicher aber ggfs. ineffizient ist die Option CL_MEM_COPY_HOST_PTR - dabei werden die Daten vom host_ptr in einen separaten Speicherbereich kopiert. Dadurch kommt es nicht zu Problemen, wenn sich am Host Pointer etwas ändert bis die Datenübertragung zum Device stattfindet - dies kann aufgrund der Asynchronität der Command Queue zeitlich versetzt passieren. Wir nutzen meist die Alternative einen Buffer mit einem host_ptr = NULL zu erstellen und dann mit clEnqueReadBuffer und clEnqueWriteBuffer vom Device zum Host bzw. zurück zu übertragen. Folgendes Beispiel (ohne Fehlerbehandlung) erzeugt ein Float-Array auf dem Host, befüllt es, legt ein entsprechenden Buffer auf dem Device an und kopiert den Inhalt vom vom Host zum Device:

float *h_data = malloc(N * sizeof(float));
for(size_t i=0; i<N;i++) 
   h_data[i] = 0.0;

// Create array in device memory
cl_mem d_data = clCreateBuffer(context, CL_MEM_READ_WRITE, N * sizeof(float), NULL, NULL);
clEnqueueWriteBuffer(commands, d_data, CL_TRUE, 0, N * sizeof(float), h_data, 0, NULL, NULL);

Der Speicher muss natürlich auch wieder freigegeben werden:

clReleaseMemObject(d_data);
free(h_data);

Zur Unterscheidung zwischen Host und Device Speicher bietet es sich an h_ und d_ als Prefixe zu verwenden.

Arbeitsaufteilung und Speicherarten#

OpenCL-Kernel werden im Normalfall auf großen Arrays angewandt - für jedes Element wird ein Work Item also eine Instanz des Kernels ausgeführt. Woher weiß das Work Item nun aber auf welchem Element es arbeiten soll? Dazu gibt es spezielle Funktionen, die im Kernel aufgerufen werden können und das spezifische Work Item identifizieren. Wir haben bereits get_global_id(0) kennengelernt - dies liefert den Index des Work Items entlang der ersten (und im bisherigen Beispiel einzigen) Dimension. Wenn wir eine doppelt verschachtelte For-Schleife in einen Kernel umwandeln wollen, würden wir wie folgt vorgehen:

for(int y = 0; y < Y; y++) {
    for(int x = 0; x < X; x++) {
        calc(array[y][x]);
    }
}

wird zum Kernel:

int y = get_global_id(0);
int x = get_global_id(1);
calc(array[y][x]);

Um diese Dimensionen von Indizes zu setzen müssen wir unseren Kernel mit der Funktion

clEnqueueNDRangeKernel(cl_command_queue queue, cl_kernel kernel, cl_uint work_dims,
                       const size_t *global_work_offset, const size_t *global_work_size,
                       const size_t *local_work_size, cl_uint num_events, 
                       const cl_event *wait_list, cl_event *event)

starten. Die Dimensionen werden spezifiziert über:

  • work_dims: die Anzahl der Dimensionen im obigen Beispiel: 2

  • global_work_size: ein Array mit work_dims Einträgen - den Größen entlang jeder Dimension. Im obigen Beispiel wäre das ein Array, das Y und X enthält

Entlang der Dimensionen kann die Größe einer Work Group im Parameter local_work_size gesetzt werden. Work Groups bestimmen zusammengehörende Work Items entlang der folgenden wichtigen Kriterien:

  • Ein OpenCL Compute Device besteht aus mehreren Compute Units (siehe oben): jede Work Group wird von einer einzelnen Compute Unit ausgeführt und jede Compute Unit führt nur eine Work Group zur gleichen Zeit aus

  • Work Items in einer Work Group können auf einen gemeinsamen schnellen lokalen Speicher zugreifen

  • Work Items in einer Work Group können synchronisiert werden - auf die Ausführungsreihenfolge verschiedener Work Groups haben wir keinen Einfluss

Der erste Punkt hat weitreichende Konsequenzen für die Performance. Betrachten wir die beispielhafte Laufzeit für Vector-Add von 2 Vektoren mit je 64 Millionen Floats auf einer Radeon Pro 5300M Graphikkarte:

local_work_size

GFLOP/s

1

1.0

2

2.0

4

3.8

8

7.1

16

11.9

NULL

12.7

Wenn wir einen NULL-Pointer übergeben, dann wählt OpenCL selbst eine Größe aus. Was wir sehen ist, dass bei kleinen local_work_sizes eine Vergrößerung einen fast linearen Leistungsanstieg bringt. Ab einer gewissen Größe stagniert die Leistung - dies liegt wie so oft daran, dass nun die Speicherbandbreite der limitierende Faktor ist und nicht mehr die Rechenleistung. Bei 12 GFLOP/s und 12 Byte Speicherzugriffen pro FLOP (2 Reads, 1 Write, jeweils 4 Byte Float) ergibt sich eine Speicherbandbreite von 144 GB/s - nahe an der theoretisch maximalen Bandbreite von 192 GB/s, wenn auch noch um Faktor 250 von der maximalen Rechenleistung von 3200 GFLOP/s (FP32) (siehe Radeon Pro 5300M Spezifikation bei TechPowerUp GPU Database).

Folgende Funktionen stehen zur Verfügung, um rauszufinden welches Work Item bearbeitet werden soll, bzw. wo in/in welcher Work Group sich das Work Item befindet:

  • uint get_work_dim(): Die Anzahl der Dimensionen

  • size_t get_global_id(uint dim): Die ID des Work Items in der Dimension dim

  • size_t get_global_size(uint dim): Die Anzahl der Work Items in der Dimension dim

  • size_t get_group_id(uint dim): Die ID der Work Group in der Dimension dim

  • size_t get_local_id(uint dim): Die ID des Work Items in der Work Group in der Dimension dim

  • size_t get_local_size (uint dim): Die Anzahl der Work Items in der Work Group in Dimension dim

  • size_t get_num_groups(uint dim): Die Anzahl der Work Groups in Dimension dim

Um Speicher in den unterschiedlichen Bereichen zu nutzen gibt es folgende Stufen:

  • Private Memory: nur für Work Item verfügbar. Sehr wenig verfügbar. Ein genaues Limit ist abhängig vom Gerät - meist gibt es eine maximale Speichergröße pro Compute Unit. Wenn pro Work Item zu viel privater Speicher genutzt wird, dann wird die Zahl der Work Items pro Compute Unit dadurch begrenzt oder ein Teil des privaten Speichers in den Local oder Global Bereich ausgelagert, der um mehrere Größenordnungen langsamer sein kann. Dies ist der default für Variablen in Kernels ohne Modifier

  • Local Memory: geteilt von allen Work Items in einer Work Group. Dieser Speicher ist ca. eine Größenordnung größer und langsamer als Private Memory. Wird für Variablen mit dem Modifier __local gewählt.

  • Global Memory: geteilt vom gesamten Kernel. Dieser Speicher ist ca. eine Größenordnung größer und langsamer als Local Memory. Modifier __global.

Das folgende Beispiel illustriert die Zwischenspeicherung eines Werts aus dem globalen Speicher im private Memory:

__kernel void lot_of_flops(__global  float *data)
{
    int i = get_global_id(0);
    float a = data[i];
    for(int j=0;j<64;j++)
    {
        a = a + 3.0;
        a = a * 8.0;
        a = a + 1.0;
        a = a * 4.0;
    }
    data[i] = a;
}
ÜBUNG: Matrix Multiplikation

Erstellen Sie eine nicht optimierte Matrix Multiplikation mit OpenCL. Optimieren Sie dann nicht benötigte Global Memory Zugriffe weg: Vorlage .

Die Nutzung von Local Memory ist etwas komplizierter, da wir nun über mehrere Work Items hinweg synchronisieren müssen. Im folgenden illustrieren wir ein mögliches Vorgehen anhand eines Beispiels: wir wollen für ein 1D Array src für jedes Element den Durchschnitt aus dem Element und den 7 folgenden Elementen berechnen und im Array dst abspeichern. Unser Kernel hierzu:

__kernel void averages(const size_t n, __global const float *src, __global float *dst)
{
    size_t x = get_global_id(0);
    float avg = 0.0;
    for(int i=0; i<8; i++)
    {
        if(x + i < n)
            avg += src[x + i];
    }
    dst[i] = avg;
}

Wir beobachten, dass z.B. der Wert von src[7] in den Berechnungen von dst[0] bis dst[7] also insgesamt 8 mal verwendet wird. Um globale Speicherzugriffe einzusparen wollen wir immer einen Teil des src Arrays im lokalen Speicher ablegen, so dass mehrere dst-Berechnungen es wiederverwenden können. Dazu gehen wir wie folgt vor:

  • Erweitern der Deklaration um einen lokalen Parameter: __kernel void averages(const size_t n, __global const float *src, __global float *dst, __local float *buf)

  • Setzen eines leeren Puffers als Argument für den lokalen Parameter: clSetKernelArg(kernel, 3, sizeof(float) * BLOCK_SIZE, NULL);

  • Im Kernel befüllen wir nun den Puffer, indem jedes Workitem sich selbst in den Puffer lädt:

size_t x = get_global_id(0);
size_t local_id = get_local_id(0);
size_t group_id = get_groupd_id(0);
buf[local_id] = src[group_id * BLOCK_SIZE + local_id];
// Alternativ: buf[local_id] = src[x];

  • Die Ladevorgänge werden in den Workitems unabhängig ausgeführt - bevor wir weitermachen müssen wir warten bis alle Workitems fertig sind: barrier(CLK_LOCAL_MEM_FENCE);

  • Nun unterscheiden wir beim Zugriff innerhalb der For-Schleife, ob der Wert im lokalen Puffer ist oder aus dem globalen Speicher gelesen werden muss:

if(i + local_id < group_size)
    avg += buf[local_id + i];
else if(x + i < n)
    avg += src[x + i];

  • Nun setzen wir beim Aufruf des Kernels noch die Größe einer Workgroup:

size_t global_size[] = {N};
size_t local_size[] = {BLOCK_SIZE};
clEnqueueNDRangeKernel(commands, kernel, 1, NULL, global_size, local_size, 0, NULL, NULL);

Hier die Laufzeiten auf einer Radeon Pro 5300M für N = 1024*1024 und verschiedene local_sizes:

local_size

Laufzeit in s

ohne

0.3554

1

0.4555

2

0.2546

4

0.1525

8

0.0783

16

0.0406

32

0.0220

64

0.0196

128

0.0149

256

0.0127

Wir können sehen, dass die zusätzlichen Conditionals und das Abspeichern bei local_size = 1 (also kein wiederverwendbares Caching) zu einer Mehrlaufzeit führt. Danach sehen wir jedoch große Gewinne. Diese werden später kleiner, da nur noch ein kleiner Anteil von Zugriffen ungecached bleibt. Beim Vergleich der Gruppengröße 256 zum Fall ohne lokales Caching haben wir eine Geschwindkeitssteigerung >25. Der Code für dieses Beispiel steht zur Verfügung/

Wann ist ein Programm für GPUs geeignet?#

Ein Programm ist dann geeignet, wenn es a) einen hohen Grad an Parallelität vorweist und b) keine Dinge tut, die GPUs Probleme bereiten:

  • Unregelmäßige Speicherzugriffe, die sich nicht gut im privaten oder lokalen Speicher cachen lassen

  • Viele Fallunterscheidungen - wir haben gesehen, dass GPUs alle Fälle ausführen müssen und ggfs. Ergebnisse wegwerfen

  • Dynamische Speicherallokationen - Speicher wird von der CPU aus gesteuert, der Transfer ist teuer

  • Rekursion oder verschachtelte Funktionsaufrufe: der Speicherbereich für Rücksprungadressen (Stack) ist sehr begrenzt