Reading Time: 11 minutes

Die GPU-Beschleunigung im wissenschaftlichen Rechnen beginnt häufig mit Bibliotheken auf hoher Ebene. Ein Forscher ersetzt Numpy durch Cupy, verwendet eine beschleunigte FFT oder ruft eine GPU-fähige lineare Algebra-Routine auf. Dies kann zu einer wesentlichen Verbesserung führen, ohne dass detaillierte Kenntnisse der GPU-Hardware erforderlich sind.

Die benutzerdefinierte GPU-Kernel-Programmierung geht eine Ebene tiefer. Anstatt eine vordefinierte Operation aufzurufen, schreibt der Entwickler die Funktion, die jeder GPU-Thread ausführt. Dies bietet Kontrolle über Thread-Indizierung, Datenlayout, Speicherübertragungen, Synchronisation, gemeinsam genutzten Speicher und die Anzahl der während jedes Kernels gestarteten Operationen.

Diese Steuerung ist nützlich, wenn ein physikalisches Modell Zugriffsmuster enthält, die nicht effizient durch gewöhnliche Array-Operationen ausgedrückt werden können. Partikelwechselwirkungen, Finite-Differenz-Schablonen, Kollisionsregeln, gitterbasierte Fluidmodelle, Materialaktualisierungen und Geometrieabfragen sind gängige Beispiele.

Ein benutzerdefinierter Kernel ist nicht automatisch schneller als eine optimierte Bibliothek. Es wird wertvoll, wenn der Algorithmus genügend parallele Arbeit aufzeigt und der Entwickler den Speicherverkehr reduzieren, Operationen kombinieren oder die Berechnung an die GPU-Architektur anpassen kann.

Was ist ein GPU-Kernel?

Ein GPU-Kernel ist eine Funktion, die von vielen leichten Threads ausgeführt wird. Jeder Thread behandelt normalerweise ein Teilchen, eine Gitterzelle, ein Element, eine Fläche, ein Pixel oder einen Eintrag in einem Array.

Das Host-Programm startet den Kernel mit einem Raster von Thread-Blöcken:

kernel[blocks, threads_per_block](arguments)

Jeder Thread bestimmt, welche Daten er aus seinen Block- und Thread-Indizes verarbeiten soll. Ein typischer eindimensionaler Index ist:

index =
    block_index * block_size
    + thread_index

Threads sind in Gruppen unterteilt, die als Blöcke bezeichnet werden. Threads in einem Block können Daten über gemeinsam genutzten Speicher synchronisieren und austauschen. Es wird allgemein erwartet, dass separate Blöcke unabhängig ausgeführt werden.

Wenn benutzerdefinierte Kernel es wert sind, geschrieben zu werden

GPU-Bibliotheken auf hoher Ebene sollten in der Regel die erste Option sein. Sie bieten bereits optimierte Matrixmultiplikation, FFTs, Reduktionen, spärliche Operationen, Zufallszahlengenerierung und viele elementweise Funktionen.

Ein benutzerdefinierter Kernel wird nützlich, wenn:

  • Die Simulation verwendet ein nicht standardmäßiges Speicherzugriffsmuster.
  • Mehrere kleine Operationen können zu einem Durchgang über Speicher verschmolzen werden.
  • Benachbarte Threads benötigen wiederholt dieselben lokalen Daten.
  • Eine spezielle Partikel-, Kollisions- oder Schablonenregel dominiert die Laufzeit.
  • Zwischen Arrays verbrauchen zu viel Gerätespeicher.
  • Die Simulation benötigt eine benutzerdefinierte differenzierbare Operation.
  • Eine vorhandene Bibliothek kann die erforderliche Grenz- oder Materiallogik nicht effizient ausdrücken.

Profilieren Sie das Programm, bevor Sie einen Kernel schreiben. Die Optimierung einer optisch komplexen Funktion, die nur einen kleinen Teil der Laufzeit ausmacht, führt nicht zu einer aussagekräftigen Gesamtbeschleunigung.

GPU-Programmiermodelle

Wissenschaftliche Software kann über mehrere Abstraktionsebenen auf GPUs zugreifen.

Sich nähern Beispiele Hauptstärke Hauptkompromiss
Direktives Offloading OpenMP, OpenACC Inkrementelle Beschleunigung des vorhandenen C, C++ oder FORTRAN-Codes Weniger direkte Kontrolle über generierte Kernel
Herstellerorientierte Kernel-Programmierung Cuda, Hüfte Detaillierte Kontrolle über die GPU-Ausführung Höhere Implementierungs- und Wartungskosten
Leistungstragbares C ++ Sycl, Kokkos, Alpaka Ein Quellmodell für mehrere Hardware-Backends Portabilität garantiert nicht überall gleiche Leistung
Python-GPU-Arrays und -Kernel Cupy, Numba-Cuda, Nvidia Warp Schnelle Entwicklung mit Zugriff auf kompilierten GPU-Code Framework-spezifische Einschränkungen und Zusammenstellungsverhalten

Enccs Einführung in die GPU-Programmierung Modelle bieten einen umfassenderen Vergleich von CUDA, HIP, OpenCL, Sycl, Kokkos, OpenMP und verwandten Ansätzen.

Cuda

CUDA bietet direkten Zugriff auf die NVIDIA-GPU-Programmierung über C++, FORTRAN und mehrere Sprachbindungen. Es stellt Gitter, Blöcke, Warps, gemeinsam genutzten Speicher, Streams, Ereignisse, Gerätebibliotheken und hardwarespezifische Optimierungstools bereit.

CUDA wird häufig ausgewählt, wenn maximale NVIDIA-spezifische Leistung und ausgereifte Tools wichtiger sind als die Portabilität der Hardware. Die Hauptkosten sind das manuelle Speichermanagement, ein Programmiermodell auf niedrigerer Ebene und die Abhängigkeit von der NVIDIA-Plattform.

Hüfte und ROCM

HIP bietet ein CUDA-ähnliches C ++ – Programmiermodell im ROCM-Ökosystem von AMD. Eine große Menge an CUDA-Code kann an HIP angepasst werden, aber die Portabilität der Quellen bedeutet nicht, dass eine Binärdatei auf jeder GPU unverändert ausgeführt wird.

Für jede Zielarchitektur ist möglicherweise weiterhin eine separate Kompilierung, Prüfung und Leistungsoptimierung erforderlich. Details finden Sie in der Official Hip Documentation .

HIP ist eine praktische Option, wenn AMD-GPUs unterstützt werden müssen oder wenn ein Projekt die Abhängigkeit von einem einzelnen Beschleunigerhersteller verringern möchte, während ein CUDA-ähnlicher Programmierstil beibehalten wird.

Cupy benutzerdefinierte Kernel

Cupy ist am bekanntesten für seine numpy-kompatiblen Arrays und beschleunigten numerischen Funktionen. Es unterstützt auch die benutzerdefinierte GPU-Programmierung.

Entwickler können verwenden:

  • ElementwiseKernel für benutzerdefinierte elementweise Ausdrücke
  • ReductionKernel für benutzerdefinierte Reduzierungen
  • RawKernel für in CUDA C++ geschriebene Kernel
  • cupyx.jit.rawkernel für JIT-Kernel-Definitionen im Python-Stil

Ein einfacher Cupy-Raw-Kernel kann wie folgt definiert werden:

import cupy as cp

scale_kernel = cp.RawKernel(
    r'''
    extern "C" __global__
    void scale_values(
        const float* input,
        float* output,
        const float factor,
        const int size
    ) {
        int index =
            blockDim.x * blockIdx.x
            + threadIdx.x;

        if (index < size) {
            output[index] =
                factor * input[index];
        }
    }
    ''',
    "scale_values"
)

size = 1_000_000
threads_per_block = 256
blocks = (
    size + threads_per_block - 1
) // threads_per_block

input_values = cp.arange(
    size,
    dtype=cp.float32
)

output_values = cp.empty_like(
    input_values
)

scale_kernel(
    (blocks,),
    (threads_per_block,),
    (
        input_values,
        output_values,
        cp.float32(2.0),
        cp.int32(size)
    )
)

Cupy ist nützlich, wenn die meisten Anwendungen bereits GPU-Arrays verwenden und nur ausgewählte Operationen benutzerdefinierte Kernel benötigen.

Numba-Cuda

Numba-Cuda kompiliert eine eingeschränkte Teilmenge von Python in GPU-Kernel. Das Programmiermodell folgt CUDA C. Die Entwickler definieren eine Funktion mit einem Kerneldekorator, berechnen einen Thread-Index und starten die Funktion mit einer Raster- und Blockkonfiguration.

Die aktive Projektdokumentation ist unter Numba-cuda verfügbar. Da sich der Entwicklungsstatus ändern kann, sollten Teams, die langlebige Software erstellen, die aktuellen Wartungs- und Migrationsrichtlinien überprüfen, bevor sie sich für das Framework engagieren.

ein Numba-Teilchenkern

Das folgende Lehrbeispiel berechnet die Gravitationsbeschleunigung mit direkten All-Pairs-Interaktionen:

import math
import numpy as np
from numba import cuda

@cuda.jit
def gravity_kernel(
    positions,
    masses,
    accelerations,
    gravitational_constant,
    softening_squared
):
    particle = cuda.grid(1)
    particle_count = positions.shape[0]

    if particle >= particle_count:
        return

    px = positions[particle, 0]
    py = positions[particle, 1]
    pz = positions[particle, 2]

    ax = 0.0
    ay = 0.0
    az = 0.0

    for other in range(particle_count):
        if other == particle:
            continue

        dx = positions[other, 0] - px
        dy = positions[other, 1] - py
        dz = positions[other, 2] - pz

        distance_squared = (
            dx * dx
            + dy * dy
            + dz * dz
            + softening_squared
        )

        inverse_distance = (
            1.0
            / math.sqrt(distance_squared)
        )

        inverse_distance_cubed = (
            inverse_distance
            * inverse_distance
            * inverse_distance
        )

        scale = (
            gravitational_constant
            * masses[other]
            * inverse_distance_cubed
        )

        ax += scale * dx
        ay += scale * dy
        az += scale * dz

    accelerations[particle, 0] = ax
    accelerations[particle, 1] = ay
    accelerations[particle, 2] = az


particle_count = 10_000
threads_per_block = 256

blocks = (
    particle_count
    + threads_per_block
    - 1
) // threads_per_block

host_positions = np.random.random(
    (particle_count, 3)
).astype(np.float32)

host_masses = np.ones(
    particle_count,
    dtype=np.float32
)

device_positions = cuda.to_device(
    host_positions
)

device_masses = cuda.to_device(
    host_masses
)

device_accelerations = cuda.device_array(
    (particle_count, 3),
    dtype=np.float32
)

gravity_kernel[
    blocks,
    threads_per_block
](
    device_positions,
    device_masses,
    device_accelerations,
    np.float32(1.0),
    np.float32(1e-4)
)

cuda.synchronize()

accelerations = (
    device_accelerations.copy_to_host()
)

Dieser Kernel hat rechnerische Komplexität von O(N²) . Es ist zur Erklärung der Thread-Mapping geeignet, aber es ist kein effizienter Produktionsalgorithmus für sehr große Partikelzahlen.

Große Gravitations- oder Partikelsysteme erfordern häufig eine Baummethode, eine schnelle Multipolmethode, räumliche Behälter oder Nachbarlisten. Eine GPU kann das Skalierungsproblem eines ineffizienten Algorithmus nicht beseitigen.

Nvidia Warp

NVIDIA Warp ermöglicht Entwicklern, stark typisierte Kernel mit Python-Syntax zu definieren. Warp kompiliert diese Funktionen in CPU- oder CUDA-Code und bietet spezielle Primitive für Geometrie, Simulation, spärliche Operationen, endliche Elemente, Optimierung und automatische Differenzierung.

Die aktuelle Nvidia Warp-Produktseite und Warp-GitHub-Repository Enthält Beispiele für Partikel, Flüssigkeiten, Netze, Optimierungen, differenzierbare Simulationen und Tile-based GPU-Programmierung.

der gleiche Partikelkern in Warp

import numpy as np
import warp as wp

wp.init()

@wp.kernel
def gravity_kernel(
    positions: wp.array(dtype=wp.vec3),
    masses: wp.array(dtype=wp.float32),
    accelerations: wp.array(dtype=wp.vec3),
    gravitational_constant: wp.float32,
    softening_squared: wp.float32
):
    particle = wp.tid()
    particle_count = positions.shape[0]

    position = positions[particle]
    acceleration = wp.vec3(0.0, 0.0, 0.0)

    for other in range(particle_count):
        if other != particle:
            displacement = (
                positions[other]
                - position
            )

            distance_squared = (
                wp.dot(
                    displacement,
                    displacement
                )
                + softening_squared
            )

            inverse_distance = (
                1.0
                / wp.sqrt(distance_squared)
            )

            inverse_distance_cubed = (
                inverse_distance
                * inverse_distance
                * inverse_distance
            )

            acceleration += (
                gravitational_constant
                * masses[other]
                * inverse_distance_cubed
                * displacement
            )

    accelerations[particle] = acceleration


particle_count = 10_000

host_positions = np.random.random(
    (particle_count, 3)
).astype(np.float32)

host_masses = np.ones(
    particle_count,
    dtype=np.float32
)

positions = wp.array(
    host_positions,
    dtype=wp.vec3,
    device="cuda"
)

masses = wp.array(
    host_masses,
    dtype=wp.float32,
    device="cuda"
)

accelerations = wp.zeros(
    particle_count,
    dtype=wp.vec3,
    device="cuda"
)

wp.launch(
    kernel=gravity_kernel,
    dim=particle_count,
    inputs=[
        positions,
        masses,
        accelerations,
        wp.float32(1.0),
        wp.float32(1e-4)
    ],
    device="cuda"
)

wp.synchronize()

Warp bestimmt intern eine geeignete Startkonfiguration für einen gewöhnlichen Kernel-Start. Fortgeschrittene Benutzer können bei Bedarf weiterhin mit Blockdimensionen, Befehlsdiagrammen, gekachelten Ausführungen und speziellen Geräteoperationen arbeiten.

Grid-Stride-Schleifen

Ein Kernel benötigt nicht für jedes Element einen fest zugewiesenen Thread. Eine Grid-Stride-Schleife ermöglicht es jedem Thread, mehrere Einträge zu verarbeiten:

from numba import cuda

@cuda.jit
def scale_with_stride(
    input_values,
    output_values,
    factor
):
    index = cuda.grid(1)
    stride = cuda.gridsize(1)

    for position in range(
        index,
        input_values.size,
        stride
    ):
        output_values[position] = (
            factor
            * input_values[position]
        )

Dieses Muster trennt die Anzahl der gestarteten Threads von der Gesamtdatengröße. Es ist nützlich, wenn sehr große Arrays verarbeitet oder eine feste Startkonfiguration wiederverwendet werden.

Eine praktische Einführung in dieses Muster finden Sie in Das Numba- und CUDA-Tutorial des Datenfroschs .

Gedächtnis verschmelzen

GPU-Kernel sind häufig eher durch Speicherbandbreite als durch den arithmetischen Durchsatz begrenzt. Threads in einem Warp sollten idealerweise auf die nahe gelegenen Speicheradressen zugreifen, damit die Hardware ihre Anfragen kombinieren kann.

Betrachten Sie Partikeldaten als gespeichert:

particle_0: x, y, z, mass
particle_1: x, y, z, mass
particle_2: x, y, z, mass

Dieses Array-of-Struktur-Layout kann für objektorientierten Code zweckmäßig sein. Ein Layout mit Struktur-of-Arrays speichert separate zusammenhängende Arrays:

x_positions[]
y_positions[]
z_positions[]
masses[]

Das zweite Layout kann eine bessere Koaleszenz liefern, wenn jeder Thread das gleiche Feld für ein anderes Teilchen liest. Das beste Layout hängt immer noch davon ab, auf welche Felder zusammen zugegriffen wird.

Vektortypen korrigieren nicht automatisch ein schlechtes Zugriffsmuster. Entwickler sollten die tatsächlichen Adressen überprüfen, die von benachbarten Threads angefordert werden.

Geteilter Speicher

Shared Memory ist ein kleiner Speicherbereich mit geringer Latenz, auf den Threads im selben Block zugegriffen werden kann. Es kann wiederholte Lesevorgänge aus dem globalen Speicher reduzieren.

Eine eindimensionale Schablone kann einen Werteblock und seinen Halo in den gemeinsamen Speicher laden:

import numpy as np
from numba import cuda, float32

BLOCK_SIZE = 256

@cuda.jit
def three_point_stencil(
    input_values,
    output_values
):
    shared = cuda.shared.array(
        shape=BLOCK_SIZE + 2,
        dtype=float32
    )

    local_index = cuda.threadIdx.x
    global_index = cuda.grid(1)
    size = input_values.size

    center = local_index + 1

    if global_index < size:
        shared[center] = (
            input_values[global_index]
        )
    else:
        shared[center] = 0.0

    if local_index == 0:
        left_index = global_index - 1

        shared[0] = (
            input_values[left_index]
            if left_index >= 0
            else 0.0
        )

    if local_index == BLOCK_SIZE - 1:
        right_index = global_index + 1

        shared[BLOCK_SIZE + 1] = (
            input_values[right_index]
            if right_index < size
            else 0.0
        )

    cuda.syncthreads()

    if global_index < size:
        output_values[global_index] = (
            shared[center - 1]
            + shared[center]
            + shared[center + 1]
        ) / 3.0

Der Aufruf von cuda.syncthreads() stellt sicher, dass jeder Thread seine Werte lädt, bevor die Schablonenberechnung beginnt.

Shared Memory sollte nicht automatisch verwendet werden. Übermäßige Zuweisung von Shared-Memory kann die Belegung verringern, die Synchronisierungskosten erhöhen und den Kernel langsamer machen. Profilierung erforderlich.

Kernfusion

Separate Array-Operationen erzeugen häufig mehrere Kernel-Starts und Zwischenarrays:

velocity += dt * acceleration
position += dt * velocity
energy = compute_energy(position, velocity)

Ein fusionierter Kernel kann alle drei Aktualisierungen berechnen, während die notwendigen Werte in Registern verbleiben. Dies reduziert den Start-Overhead und den globalen Speicherverkehr.

Fusion ist am vorteilhaftesten, wenn die Operationen einfach und speichergebunden sind. Zu viel Arbeit kann die Nutzung des Registers erhöhen, die Belegung verringern und die Wartung des Kernels erschweren.

Blockgröße wählen

Werte wie 128 oder 256 Threads pro Block sind vernünftige Ausgangspunkte, nicht universelle Optima.

Die beste Blockgröße hängt ab von:

  • Register pro Thread verwendet
  • Shared Memory pro Block verwendet
  • Zweigdivergenz
  • Anleitungsmix
  • Speicherzugriffsverhalten
  • Die Ziel-GPU-Architektur

Eine Regel wie das Starten von mindestens doppelt so vielen Blöcken wie Streaming-Multiprozessoren kann ein nützliches anfängliches Experiment sein, garantiert jedoch keine maximale Leistung. Belegungsrechner und Profilierungswerkzeuge sollten die endgültige Konfiguration leiten.

differenzierbare Physik

Die differenzierbare Simulation berechnet, wie sich eine Ausgabe in Bezug auf Eingaben wie Materialeigenschaften, Kräfte, Geometrie oder Anfangsbedingungen ändert.

Warp kann unterstützte Kerneloperationen aufzeichnen und die automatische Differenzierung des Reverse-Modus ausführen. Dies kann für inverse Probleme, Designoptimierung, Parameterschätzung und Integration in maschinelle Lern-Workflows verwendet werden.

Ein vereinfachtes Muster verwendet ein Band:

with wp.Tape() as tape:
    wp.launch(
        kernel=simulation_kernel,
        dim=element_count,
        inputs=[state, parameters],
        outputs=[result],
        device="cuda"
    )

    wp.launch(
        kernel=loss_kernel,
        dim=element_count,
        inputs=[result, target],
        outputs=[loss],
        device="cuda"
    )

tape.backward(loss)

Eine automatische Differenzierung ist nicht für jeden Kernel garantiert. In-Place-Overwrites, nichtdeterministische Atome, externer nativer Code, diskontinuierliche Logik und nicht unterstützte Operationen erfordern möglicherweise eine Neuformulierung oder benutzerdefinierte Gradienten.

Diese Fähigkeit verbindet sich direkt mit Anwendungen wie Adjoint Optimization und Physikinformierte neuronale Netze .

Wenn eine GPU möglicherweise nicht schneller ist

Es gibt keinen universellen Schwellenwert für die Problemgröße, bei dem eine GPU schneller wird als eine CPU. Der Crossover hängt von Hardware, Präzision, Datenbewegung, Algorithmenstruktur, Compilerqualität und der Häufigkeit der Wiederverwendung derselben geräteresidenten Daten ab.

Arbeitsanfall GPU-Eignung Hauptüberlegung
Große Partikelaktualisierung Oft stark viele ähnliche unabhängige Operationen
Große reguläre Schablone Oft stark Vorhersehbarer paralleler Speicherzugriff
Dichte lineare Algebra Stark bei der Verwendung optimierter Bibliotheken Hohe arithmetische Intensität
Kleine Simulation Hardwareabhängig Start und Transfer-Overhead können dominieren
Unregelmäßige Graphendurchquerung Gemischt Divergenz und unvorhersehbarer Speicherzugriff
I/O-gebundener Workflow in der Regel begrenzt Die GPU kann einen Speicher- oder Netzwerkengpass nicht entfernen
stark sequentieller Algorithmus normalerweise schwach Unzureichende unabhängige Arbeit

Der richtige Workflow besteht darin, die CPU-Version, die GPU-Bibliotheksversion und die benutzerdefinierten Kernel-Version mit repräsentativen Daten zu messen. Forscherhandbuch für wissenschaftliche Rechenleistung Erklärt, wie Sie den tatsächlichen Engpass identifizieren, bevor Sie ihn optimieren.

Vermeidung von Transfer-Overhead

Wiederholte Übertragungen zwischen Host- und Gerätespeicher können den Nutzen der GPU-Berechnung beseitigen.

Eine ineffiziente Schleife kann diesem Muster folgen:

  1. Daten auf die GPU kopieren.
  2. Führen Sie einen kleinen Kernel aus.
  3. Kopieren Sie das Ergebnis in die CPU.
  4. Ändern Sie es auf der CPU.
  5. Kopieren Sie es wieder auf die GPU.

Ein besseres Design behält den Simulationszustand über viele Zeitschritte auf dem Gerät und überträgt nur die Ausgabe, die für Visualisierung, Checkpointing oder Analyse erforderlich ist.

Asynchrone Streams, festgesteckter Host-Speicher, überlappende Kommunikation und Unified-Memory-Funktionen können helfen, sollten jedoch erst eingeführt werden, nachdem die normalen Übertragungskosten gemessen wurden.

Die XLB-Fallstudie

Das XLB-Projekt von Autodesk Research bietet ein nützliches Beispiel für die Python-native GPU-Simulation. XLB ist eine Open-Source-Lattice-Boltzmann-Bibliothek mit mehreren rechnerischen Backends, darunter Nvidia Warp.

In von Autodesk Research und NVIDIA gemeldeten Benchmark-Konfigurationen erreichte das Warp-Backend eine Leistung in der Nähe der verglichenen C ++ / OpenCL FluidX3D-Implementierung für einen bestimmten Lid-gesteuerten Hohlraumfall. Ein separater Vergleich ergab eine ungefähre achtfache Beschleunigung über das JAX-Backend von XLB bei ausgewählten Hardware und Konfigurationen.

Das Team zeigte auch einen Out-of-Core-Ansatz für einen GH200-Cluster mit acht Knoten mit einer Domäne von ungefähr 50 Milliarden Gitterzellen. Diese Ergebnisse gelten für den gemeldeten Solver, Problem, Hardware und Implementierungsoptionen. Sie sollten nicht als allgemeine Beschleunigungsgarantien für Python-Kernel behandelt werden.

Die vollständige Fallstudie ist in verfügbar. Autodesk-Forschung bringt die Warp-Geschwindigkeit in die rechnerische Fluiddynamik auf NVIDIA GH200 . Der aktuelle Projektcode ist im Autodesk XLB-Repository verfügbar.

Wie man einen Kernel richtig bewertet

GPU-Operationen sind normalerweise asynchron. Das Messen nur der Python-Funktionsdauer kann die Zeit angeben, die erforderlich ist, um den Kernel in die Warteschlange zu stellen, und nicht die Zeit, die für die Ausführung erforderlich ist.

Ein grundlegender Benchmark sollte:

  1. Führen Sie den Kernel mehrmals aus, um die Kompilierung und das Aufwärmen auszulösen.
  2. Synchronisieren, bevor Sie den Timer starten.
  3. Führen Sie mehrere gemessene Iterationen aus.
  4. Synchronisieren, bevor der Timer gestoppt wird.
  5. Berichten Sie über Durchschnitt und Variation über Wiederholungen.
  6. Trennen Sie die Datenübertragungszeit von der Kernel-Ausführungszeit.
  7. Stellen Sie sicher, dass CPU- und GPU-Versionen gleichwertige Ergebnisse liefern.
import time
from numba import cuda

# Warm-up and JIT compilation
kernel[blocks, threads](*arguments)
cuda.synchronize()

start = time.perf_counter()

for _ in range(100):
    kernel[blocks, threads](*arguments)

cuda.synchronize()

elapsed = time.perf_counter() - start
average = elapsed / 100

print("Average kernel time:", average)

Der Vergleich sollte eine realistische Simulationsgröße verwenden und den vollständigen Workflow einschließen, wenn die Kosten für die Übertragung oder Vorverarbeitung von Bedeutung sind.

Die Arbeitspräzisionsmethode kann auch die Laufzeit mit numerischen Fehlern vergleichen. Die SCIML-Benchmark-Dokumentation enthält Beispiele für diesen Bewertungsstil.

Ein praktischer Entwicklungsworkflow

  1. Implementieren und überprüfen Sie eine klare CPU-Referenzversion.
  2. Profil Die Anwendung, um die dominante Operation zu lokalisieren.
  3. Versuchen Sie es mit einer optimierten GPU-Bibliothek, bevor Sie einen Kernel schreiben.
  4. Bewahren Sie häufig wiederverwendete Daten auf dem Gerät auf.
  5. Schreiben Sie den einfachsten korrekten Kernel.
  6. Überprüfen Sie es gegen das CPU-Ergebnis.
  7. Messen Sie Speicherübertragungen und Ausführung separat.
  8. Überprüfen Sie die Koaleszenz, Belegung, Divergenz und Registernutzung.
  9. Testen Sie den gemeinsamen Speicher oder die Fusion nur, wenn die Profilerstellung sie unterstützt.
  10. Vergleichen Sie mehrere Problemgrößen und GPU-Architekturen.

Auswahl eines Frameworks

Anforderung Möglicher Ausgangspunkt
Vorhandener GPU-Workflow im Numpy-Stil Cupy mit integrierten Operationen oder benutzerdefinierten Kernel
Kleine Anzahl von Python-Kernel im Cuda-Stil NUMBA-CUDA, nachdem er seinen aktuellen Support-Status überprüft hat
Simulation, Geometrie und differenzierbare Kernel Nvidia Warp
Maximale NVIDIA-spezifische Steuerung CUDA C ++
AMD-GPU-Ziel mit CUDA-ähnlicher Quelle Hüfte und ROCM
Portable C++ über mehrere Backends SYCL, KOKKOS oder eine andere Performance-Portabilitätsschicht

Häufige Fehler beim Programmieren von Kernel

  • Verschieben von Daten zwischen CPU und GPU innerhalb jedes Zeitschritts
  • Ignorieren von Out-of-Bounds-Threads
  • Verwenden eines ineffizienten Algorithmus und Erwarten von Hardware, um seine Skalierung zu reparieren
  • Zugriff auf Speicher mit einem nicht koaleszierten Layout
  • Hinzufügen eines gemeinsamen Speichers, ohne zu messen, ob es hilft
  • Viele winzige Kernel starten, anstatt Fusion in Betracht zu ziehen
  • Verwenden übermäßiger Register oder gemeinsam genutzter Speicher pro Block
  • Benchmarking asynchroner Code ohne Synchronisation
  • Vergleichen von Ausgaben ohne numerische Genauigkeit zu überprüfen
  • Verwenden von festen Leistungsansprüchen über verschiedene GPUs und Workloads
  • Angenommen, die tragbare Quelle bietet tragbare Leistung
  • Anwenden der automatischen Differenzierung auf nicht unterstützte In-Place-Operationen

Verwandte Anleitungen

Schlussfolgerung

Benutzerdefinierte GPU-Kernel ermöglichen es wissenschaftlichen Entwicklern, Partikel, Gitterzellen, Materialregeln und andere physikalische Operationen direkt in massiv parallele Funktionen zu übersetzen. Sie bieten Kontrolle über Indexierung, Speicherzugriff, Synchronisation, Kernelfusion und gerätespezifische Optimierung.

Cupy bietet sowohl hochrangige Array-Operationen als auch benutzerdefinierte Kernel-Schnittstellen. NUMBA-CUDA bietet ein CUDA-ähnliches Modell in Python, während Nvidia Warp Simulationsprimitive, Kacheloperationen und automatische Differenzierung hinzufügt. CUDA und HIP bieten C++-Steuerung auf niedrigerer Ebene, und leistungstragbare Frameworks unterstützen breitere Hardware-Ziele.

Die größten Leistungssteigerungen ergeben sich selten allein aus einer sich ändernden Syntax. Sie stammen aus der Auswahl eines parallelen Algorithmus, der Speicherung von Daten auf dem Gerät, der Reduzierung des Speicherverkehrs, der Verwendung eines geeigneten Datenlayouts und der Eliminierung unnötiger Zwischenoperationen.

Ein benutzerdefinierter Kernel sollte erst entwickelt werden, nachdem die Profilerstellung einen echten Engpass identifiziert hat. Es sollte anhand einer Referenzlösung validiert und mit Synchronisation, repräsentativen Problemgrößen und vollständiger Abrechnung der Speicherübertragungskosten verglichen werden.

Wenn diese Bedingungen erfüllt sind, können Python-basierte Kernel-Tools ernsthafte Physiksimulationen unterstützen, ohne die Forscher zu zwingen, jeden Teil der Anwendung in C ++ auf niedriger Ebene zu verschieben.