Neu in Mathematica 8 – Mathematica und Cuda

Mathematica und Cuda

Über das integrierte CUDALink Paket liefert Mathematica Unterstützung für die Grafikprozessorprogrammierung für grafikprozessorbeschleunigte lineare Algebra, diskrete Fourier-Transformationen und Bildverarbeitungsalgorithmen. Bei Bedarf können Sie CUDALink Module auch ganz einfach selbst schreiben.

Das CUDALink Paket von Mathematica liefert folgende Funktionen ganz ohne zusätzliche Kosten:

  • Zugriff auf die automatischen Interface-Builder, Import/Export-Funktionen und Visualisierungsfunktionen von Mathematica
  • Zugriff auf die betreuten Mathematica Datenbanken
  • Unterstützung für Arithmetik mit einfacher und doppelter Genauigkeit
  • Möglichkeit, benutzerdefinierte CUDA Programme in Mathematica zu laden
  • Skalierbarkeit für mehrere Geräte
  • Integration mit existierenden Mathematica Technologien wie Remote Sessions

Mehr Infos zu High Performance Computing mit Wolfram Mathematica.

Weitere Ressourcen:

Einleitung

Die "Common Unified Device Architecture" (CUDA) wurde im Jahre 2007 von NVIDIA entwickelt. Auch wenn es die Programmierung der GPU schon seit Jahren gibt, haben Hardwarerestriktionen das Schreiben von Software für die Grafikkarte nicht leicht gemacht. Eine Aufgabe von CUDALink von Mathematica ist es, das Programmieren und die Portierung der eigenen Software zu vereinfachen. Wenn man Mathematica benutzt, muss man sich keine Gedanken über die vielen Hürden machen. Somit bleibt das Hauptaugenmerk auf dem Schreiben der CUDA Kernel. Das wird durch durch Nutzen aller NVIDIA Architektur Stacks bewerkstelligt.

Architektur

Hier der schematische Aufbau der CUDA Architektur.

Grid und Blocks

CUDA-Programmierung basiert auf dem Daten-Parallelisierungs-Modell. Berechnungen werden auf viele hundert - wenn nicht gar tausend Threads geteilt. Wenn Folgendes Ihre Berechnung wäre:

OutputData = Table[fun[InputData[[i, j]]], {i, 10000}, {j, 10000}]

, dann wäre dies hier Ihr CUDA-Äquivalent:

CUDALaunch[fun, InputData, OutputData, {10000, 10000}]

wobei "fun" Ihre selbst definierte Funktion ist. Das Beispiel startet 1000 x 1000 Threads und übergibt jedem deren Indizes. Das Ergebnis wird InputData zugewiesen und in OutputData gespeichert.

Der Grund, weshalb CUDA Tausende von Threads starten können, liegt an der Hardwarearchitektur. Im folgenden Abschnitt erörtern wir nun, wie Threads geteilt werden.

CUDA teilt programmatische Probleme in ein-, zwei- oder dreidimensionale Grids auf. Folgendes Bild zeigt eine typische zweidimensionale CUDA-Thread-Konfiguration.

Jedes Grid beinhaltet mehrere Blöcke, und jeder Block besteht aus Threads. Ausgehend von dem obigen Bild, kann man sich die Grids, Blocks und Threads wie folgt vorstellen.

Ob man nun eine ein-, zwei- oder dreidimensionale Thread-Konfiguration verwendet, ist abhängig von dem Problem. Im Fall von Bildbearbeitung mappen wir die Threads wie in dem Bild und weisen jedem Pixel eine Funktion zu.

Der eindimensionale zelluläre Automat kann auf ein eindimensionales Grid gemappt werden.

Program Cycle

Das Wesentliche an der CUDA-Programmierung sind folgende drei Schritte:

  1. Kopieren der Daten beim Starten der Threads
  2. Warten, bis die GPU die Berechnung durchgeführt hat
  3. Kopieren des Ergebnisses von der GPU zum Host
Das obere Bild zeigt den typischen Ablauf eines CUDA-Programms im Einzelnen.
  1. Zuteilen des Arbeitsspeichers der GPU: GPU- und CPU-Speicher sind physikalisch getrennt, und der Programmierer muss die zugewiesenen Kopien verwalten.
  2. Der Speicher wird von der CPU zur GPU kopiert.
  3. Konfigurieren der Thread-Konfiguration: Der korrekte Block und die Grid-Dimension für das Problem werden ausgewählt.
  4. Die konfigurierten Threads werden gestartet.
  5. Durch die Synchronisation der CUDA-Threads wird sicher gestellt, dass die GPU alle Aufgaben beendet hat, um dann weitere Operationen auf dem Speicher der GPU auszuführen.
  6. Wenn die Operationen ausgeführt wurden, wird der Speicher von der GPU auf die CPU wieder zurückgeschrieben.
  7. Der GPU-Speicher wird geleert.

Wenn Mathematica genutzt wird, muss man sich über die vielen verschiedenen Schritte keine Gedanken mehr machen, sondern kann sich auf das Schreiben des CUDA Kernels konzentrieren.

Speicher-Hierarchie

CUDA-Speicher wird in verschiedene Bereiche getrennt; jeder dieser Speicherbereiche hat seine Vorteile und Beschränkungen. Das folgende Bild zeigt alle Typen von CUDA-Speicher.

Global Memory

Der Arbeitsspeicher beansprucht den größten Raum auf einer Grafikkarte und ist dabei auch am langsamsten. Alle Threads können Speicher aus dem "Global Memory" nehmen. Dies sollte man aber aus Gründen der Performance, wenn möglich, vermeiden.

Texture Memory

Das "Texture Memory" residiert in demselben Bereich wie das "Global Memory", ist aber nur lesbar. Dafür ist es sehr viel performanter. Dieser Speicher unterstützt nur die Werte string, char und int unterstützt.

Constant Memory

Hierbei handelt es sich um einen schnellen Speicher, der für jeden Thread aus dem Grid verfügbar ist. Der Speicher wird gecached, ist aber global auf 64 kB limitiert.

Shared Memory

Auch dies ist ein schneller Speicher, der mit jeweils einem Block verbunden ist. Unter aktueller Hardware ist dieser Speicher auf 16kB pro Block limitiert.

Local Memory

Dieser Speicher ist für jeden Thread lokal, residiert aber im "Global Memory", bis der Kompiler die Variablen in das Register schreibt. Auch hier sollte man die Zugriffe auf ein Minimum beschränken, selbst wenn dieser Speicher nicht ganz so langsam wie das "Global Memory" ist.

Compute Capability

Die "Compute Capability" gibt an, was die verwendete Hardware für Operationen unterstützt. Zurzeit sind folgende verfügbar:

Compute Capability:
Extra Features:
1.0 base implementation
1.1 atomic operations
1.2 shared atomic operations and warp vote functions
1.3 support for double-precision operations
2.0 double precision, L2 cache, and concurrent kernels

Sollten Sie CUDALink noch nicht importiert haben, holen Sie dies bitte jetzt nach:

Needs["CUDALink`"]

Folgender Befehl zeigt alle CUDA-fähigen Komponenten Ihres Systems an:

TabView[Table[
CUDAInformation[ii, "Name"] ->
CUDAInformation[ii,
"Compute Capabilities"], {ii, $CUDADeviceCount}]]

CUDA Devices

Wenn die Systemhardware es erlaubt, lässt CUDA den Anwender entscheiden, auf welcher GPU die Berechnungen durchgeführt werden sollen. Standardgemäß wird die schnellste genutzt, aber diese Einstellung kann von dem Benutzer außer Kraft gesetzt werden. Sobald eine GPU ausgewählt wurde, kann diese während der Kernel Session nicht mehr verändert werden.

 

Der erste Kernel

CUDA Kernel sind atomare Funktionen, die die Berechnung auf jedes Element einer Eingabeliste durchführen. Unser erster Kernel wird 2 Elemente miteinander addieren.

__global__ void addTwo_kernel(int * arry, int len) {
int index = threadIdx.x + blockIdx * blockDim.x;
if (index >= len) return;
arry[index] += 2;
}

__global __  ist ein Funktionsqualifizierer, der dem Compiler mitteilt, dass diese Funktion auf der GPU ausgeführt werden soll. __global__ functions kann direkt aus C heraus aufgerufen werden. Die andere Funktion ist __device__, welche Funktionen definiert, die von anderen Funktionen wie __global__ oder __device__ aufgerufen werden können, aber nicht von C.

Ein CUDA Kernel muss "void"-Ausgabe besitzen; um ein Ergebnis zu erhalten, müssen daher die Input Pointer überschrieben werden. In diesem Fall übergeben wir ein Array und überschreiben die Elemente (hier kann man sich vorstellen, dass das Array die Eingabe/Ausgabeparameter besitzt)

int index = threadIdx.x+blockIdx*blockDim.x;

Mit dem obigen Source Code erhalten wir den Index. Der CUDA Kernel stellt uns folgende Variablen zur Verfügung, die abhängig von dem auszuführenden Grid und der Blockgröße gesetzt werden.

  • threadIdx - Index des aktuellen Threads; der Threadindex liegt zwischen 0 und blockDim -1
  • blockIdx - Index des aktuellen Blocks; der Blockindex liegt zwischen 0 und gridDim -1
  • blockDim - die Dimension der Blockgröße
  • gridDim - die Dimension der Gridgröße

Diese Parameter werden von CUDA automatisch gesetzt in Abhängigkeit der Parameter, die beim Kernel-Start übergeben werden (die Block- und Griddimension). Höhere Dimensionen werden automatisch auf 1 gesetzt, wenn man eine Berechnung mit niedrigeren Dimensionen startet. Wenn man somit eine eindimensionale Berechnung ausführt, werden .y und .z auf 1 gesetzt.

In den meisten Fällen benutzt man in einem 1D-Grid threadIdx.x + blockIdx.x * blockDim.x, um das globale Offset zu finden. Um die globale Position in 2D zu erhalten, benutzt man threadIdx.x + blockIdx.x * blockDim.x + (threadIdx.y + blockIdx.y * blockDim.y)*width.

if (index >= len) return;

Da Threads im Mehrfachen der Blockdimension ausgeführt werden, muss sich der Anwender vergewissern, die Grenzen nicht zu überschreiben, sofern die Eingabedimension nicht das Mehrfache der Blockdimension ist. Dies wird sichergestellt durch:

in[index] += 2;

Diese Funktion wird auf jedes Element der Liste angewandt; hier wird 2 zu jedem Element addiert.

Der zweite Kernel

Der zweite Kernel implementiert eine finite Differenziermethode

__device__ float f(float x) {
return tanf(x);
}

__global__ void secondKernel(float * diff, float h, int listSize) {
int index = threadIdx.x + blockIdx.x * blockDim.x;
float f_n = f(((float) index) / h);
float f_n1 = f((index + 1.0) / h);
if( index < listSize) {
diff[index] = (f_n1 - f_n) / h;
}
}

__device__ float f(float x) {

Dies ist eine Device Funktion. __device__ Funktionen können von der __global__ Funktion aufgerufen werden, sind aber nicht aus Mathematica direkt ausführbar.

return tanf(x);

Durchführung des Tangens auf eine float-Eingabe.

__global__ void secondKernel(float*diff, float h, int listSize) {

Funktion, die die Liste diff ausgibt und die skalaren Größen h und listSize entgegen nimmt.

int index=threadIdx.x+blockIdx.x*blockDim.x;

Liefert den globalen Thread Offset.

float f_n=f(((float) index)/h);

Evaluiert die Funktion bei Index/h. Da Index und h vom Typ int sind, müssen wir das Ergebnis nach float konvertieren, um einen Abbruch zu vermeiden.

float f_n1=f((index+1.0)/h);

Evaluiert die Funktion bei (index+1)/h.

if(index<listSize) {

Verhindert das Überschreiben des Ausgabebuffers.

diff[index]=(f_n1-f_n)/h;

Berechnet die Differenz und speichert das Ergebnis in diff.

Kernel in Mathematica laden

Sollten Sie den CUDALink noch nicht geladen haben, holen Sie dies bitte jetzt nach.

Needs["CUDALink`"]

Nun speichern wir den erstellten Code in einen string.

secondKernelCode = "
__device__ float f(float x) {
return tanf(x);
}

__global__ void secondKernel(float * diff, float h, int listSize) {

int index = threadIdx.x + blockIdx.x * blockDim.x;
float f_n = f(((float) index) / h);
float f_n1 = f((index + 1.0) / h);
if( index < listSize) {
diff[index] = (f_n1 - f_n) / h;
}
}";

Dies lädt die CUDA-Funktion. Wir übergeben den Kernel-Code als erstes Argument, den Kernel-Namen, welcher ausgeführt wird ("secondKernel"), als zweites Argument, die Funktionsparameter als drittes Argument und die Blockgröße als viertes Argument. Das Ergebnis wird in secondKernel gespeichert.

secondKernel =
CUDAFunctionLoad[secondKernelCode,
"secondKernel", {{"Float"}, "Float", _Integer}, 16]

9

Im Hintergrund passieren folgende Dinge:

  • CUDA-Verfügbarkeit wird geprüft.
  • Der CUDA-Code wird zu einem Binary kompiliert und für die jeweilige GPU optimiert.
  • Der CUDA-Code wird gecached, um später Kompilierung zu verhindern.
  • Es wird geprüft, ob der Kernel aus dem CUDA-Code existiert.

Einmal geladen, kann die CUDA-Funktion wie eine Mathematica-Funktion benutzt werden. Bevor wir starten können, benötigen wir einen Speicher, in dem das Ergebnis abgelegt wird.

buffer = CUDAMemoryAllocate["Float", 1024]

Hier führen wir den zweiten CUDA Kernel aus.

secondKernel[buffer, 100.0, 1024]

Mit dem folgenden Befehl erhalten wir die ersten 30 Elemente aus dem Ausgabebuffer.

CUDAMemoryGet[buffer][[;; 30]]

o_5

Dies zeichnet unser Ergebnis:

ListPlot@CUDAMemoryGet[buffer]

o_6

Der Buffer muss vom Benutzer wieder gelöscht werden.

CUDAMemoryUnload[buffer];

Systemanforderung

Zur Verwendung von Mathematica's CUDALink ist erforderlich:
  • Betriebssystem: Windows, Linux oder Mac OS X, 32-Bit oder 64-Bit Architektur
  • NVIDIA CUDA-fähige Produkte
  • Mathematica 8.0 oder besser

Die leistungsfähigen GPU Computing Funktionen von Mathematica wurden auf Tesla und Quadro GPU Computing Produkten entwickelt und erfordern aktuelle CUDA-fähige NVIDIA Grafikprozessoren.

Tesla und Quadro GPU Computing Produkte werden speziell für maximale Rechenleistung und höchste numerische Präzision entwickelt und werden von führenden Herstellern professioneller Systeme angeboten und unterstützt.

Empfohlene Konfigurationen

Highend Workstation

  • Zwei Tesla C2050 oder C2070 Grafikprozessoren
  • Quadro NVS 295
  • zwei Quad-Core CPUs
  • 12 GB Systemspeicher
tesla-badge
Mid-Range Workstation
  • Tesla C2050 oder C2070 Grafikprozessor
  • Quad-core CPU
  • 8 GB Systemspeicher
tesla-badge
Einstiegs-Workstation
  • Quadro 4000 Grafikprozessor
  • Quad-Core CPU
  • 4 GB Systemspeicher
quadro-badge

© ADDITIVE GmbH. Alle Rechte, Irrtümer und Änderungen vorbehalten.

Impressum  /  Datenschutz  /  AGB