Physikalische Optik auf dem Nvidia-GPU

Intro

Die meiste Zeit meiner Laufbahn in der Optik habe ich mich mit der Simulation der Abbildung von räumlichen Lichtmodulatoren beschäftigt, vor allem mit SLMs, die mit teilkohärenten Quellen beleuchtet werden. Normalerweise geht das mit der FFT recht schnell, zumindest für ein einziges Bild. Manchmal muss man jedoch eine ganze Reihe von Simulationen durchführen. Ein solches Beispiel ist die Arbeit an den optischen Anforderungen für ein optisches System. Ein Ansatz, der sich als nützlich erwiesen hat, ist die Verwendung von Latin Hypercube Stichprobenverfahren. Bei einem optischen System würde man wahrscheinlich Abbildungsfehler bis zur 5. Ordnung einbeziehen, vielleicht auch etwas mehr, und das wird eine ziemlich große Simulation. Vielleicht 1000 oder 2000 Bilder. Dann kann man sich die Daten ansehen. Vielleicht war der untersuchte Raum (Bereich) zu klein oder zu groß, und man sollte diese Simulationen ein paar Mal wiederholen.

Auch wenn es der Computer ist, der diese Aufgabe erledigt, ist es nicht immer einfach, eine andere Aufgabe einzuplanen, und selbst wenn dies möglich ist, ist ein halbstündiger oder stündlicher (oder vielleicht auch mehr) Wechsel der Aufgabe nicht sehr produktiv. Daher ist es wichtig, zweidimensionale FFTs schnell berechnen zu können, und zwar sehr wichtig.

Die Struktur einer Bildsimulation

Die Herangehensweise an teilkohärente Bildsimulationen hängt in der Regel davon ab, ob Sie vorhaben, viele Bilder mit einer bestimmten Verteilung der Quellen zu berechnen, dann wäre Ihr Ansatz vielleicht SOCS (Sum Of Coherent Systems), oder wenn Sie nach einer guten Intensitätsverteilung suchen wollen, lohnt sich der Aufwand mit SOCS nicht und der Ansatz wird sein, nur die inversen Transformationen der von Ihrer Maske reflektierten Amplituden zu berechnen oder, wie es bei mir meist der Fall war, SLM. Jede dieser inversen Transformationen stellt ein Bild dar, dessen Intensität aufsummiert werden muss, um die Gesamtbelichtung zu ermitteln.

Je nach Größe des periodischen Bereichs und des Pupillenfüllfaktors der Beleuchtung (in der Regel Sigma genannt) können wir leicht mehrere hundert inverse Transformationen benötigen. An diesem Punkt kann man beschließen, eine gröbere Abtastung der Lichtquelle vorzunehmen, aber ich würde nicht wesentlich unter 100 Beleuchtungspunkte gehen. Wenn wir auf 25 (zum Beispiel) heruntergehen, wird es definitiv sichtbare Beiträge der Periodizität der Kohärenzfunktion geben. Selbst 50 sind zu viel, so dass 100 als Faustregel keine schlechte Zahl ist. Manchmal ist die Abtastdichte durch einen Homogenisator gegeben, so dass wir nicht viel mehr hinzufügen können.

Zusammenfassend lässt sich also sagen, dass die grundlegende Bildsimulation, über die ich hier spreche, aus einer Vorwärtstransformation zur Berechnung der Amplitude besteht, die die Pupille des Systems beleuchtet, sowie aus einer Reihe von Rücktransformationen, die nur triviale Übersetzungen einer Transformation sind, um die Teilbilder zu berechnen. Dazwischen gibt es eine Pupillenfilterung, aber das ist ein rechnerisch sehr billiger Vorgang.

CUDA-API

Zunächst war die CUDA-API mag ein wenig einschüchternd wirken, und selbst mit meiner beträchtlichen Zeit, die ich mit diesen Dingen verbracht habe, war ich ein wenig besorgt, wie viel Schmerz ich mir antun müsste, bevor ich die Lernkurve verlasse. Erschwerend kam hinzu, dass ich einen Code wollte, der zusammen mit Octave oder Matlab verwendet werden kann, und zwar in C. Dieser Code muss auch identische Ergebnisse liefern, wenn er für Benutzer kompiliert wird, die keinen Zugang zu CUDA haben.

Zu meiner Überraschung war der Einstieg in CUDA viel einfacher, als ich erwartet hatte, und Dinge, die ich befürchtet hatte, wie z. B. endlos oft neu starten zu müssen, weil ich es geschafft hatte, meine Grafikkarte zu sperren, sind nie passiert.

Einige strategische Entscheidungen

Der erste Entscheidungspunkt in Bezug auf die API war, ob ich die CUDA FFTW-Wrapper-Bibliothek verwenden sollte oder mich ein wenig zurückziehen und direkt die cuFFT-Bibliothek verwenden sollte. Wenn Sie eine eigenständige Anwendung haben, die auf der FFTW-Bibliothek basiert, dann ist das vielleicht eine gute Wahl, auch wenn es einige offensichtliche Dinge gibt, die dieser Wrapper einfach nicht kann, wie z. B. die Zuweisung von Geräte- und Hostpuffern für Ihre Daten. Im CUDA-Jargon befindet sich der "Geräte"-Speicher auf der GPU und der Host-Speicher ist an Ihre CPU gebunden.

Bei Octave/Matlab kommt noch hinzu, dass die Verwendung der CUDA FFTW-Wrapper-Bibliothek zu Namenskonflikten zwischen FFTW und gleichnamigen Funktionen in der Wrapper-Bibliothek führen wird. Selbst wenn wir uns entscheiden, Implementierungen für alle Octave/Matlab-FFT-Varianten bereitzustellen, werden wir immer noch nicht wissen, ob es nicht irgendwo eine Binärdatei gibt, die erwartet, einen Plan von fftw_plan_dft_2d zu erhalten, und während die Wrapper-Bibliothek eine solche Funktion bereitstellt, sind die Bibliotheken nicht binärkompatibel. Daher ist es eine gute Idee, den Konvertierungsleitfaden zu cuFFT zu verwenden.

Die einfachste Lösung

Fangen wir also ganz einfach an. Was ist das absolute Minimum, das wir tun müssen, um FFTs auf der GPU auszuführen? Eigentlich ist dafür nicht viel erforderlich. Wir müssen

• Zuweisung von GPU-Speicher - device_mem = cudaMalloc(n);
• Einen Plan erstellen - cufftPlan2D(&plan, n1, n2, CUFFT_Z2Z);
• Kopieren von Host-Daten auf die GPU - cudaMemcpy(device_mem, host_mem, n, cudaMemcopyDefault);
• Führen Sie den Plan aus - cufftExecZ2Z(plan, device_mem_in, device_mem_out, CUFFT_FORWARD);
• Das Ergebnis zurückbringen - cudaMemcpy(host_mem, device_mem, n, cudaMemcpyDefault);

Offensichtlich bin ich hier ziemlich schlampig, aber wenn wir so vorgehen, können wir zwar die FFT auf der GPU durchführen, aber es wäre auch eine kleine Enttäuschung. Die meiste Zeit würde für das Kopieren von Daten zum und vom Gerät verwendet werden. Die Chancen stehen gut (wenn "gut" überhaupt das richtige Wort ist), dass wir gar nicht bemerken würden, dass die GPU überhaupt etwas tut, außer dass sie tatsächlich das richtige Ergebnis berechnet.

Was mussten wir also mit den Daten tun, die wir von der GPU zurückkopiert hatten? Das Einzige, was wir tun mussten, war, das Produkt aus der Transformation und ihrer Konjugierten für jede Probe zu berechnen und die Summe zu akkumulieren, bis wir mit dem Scannen der Beleuchtungsquellen fertig sind. Wenn wir das auf der GPU machen könnten, müssten wir das Ergebnis wenigstens nicht mehr als einmal zurücklesen. Das bedeutet, dass wir lernen müssen, wie man einen trivialen Code auf dem Grafikprozessor ausführt.

Ein Stück Code, das das tun würde, ist:

__global__ void _launchAddIntensity2Df(float *img, const cufftComplex *amp, float factor, int n)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        float re = amp[idx].x;
        float im = amp[idx].y;
        img[idx] += (re * re + im * im) * Faktor;
    }
}

extern "C" void addIntensity2Df(float *img, const cufftComplex *amp, float factor, int n) {
    int blockSize = 256;
    int gridSize = (n + blockSize - 1) / blockSize;
    _launchAddIntensity2Df<<>>(img, amp, factor, n);
}

Der Prototyp für addIntensity2Df muss in einen Header exportiert werden und kann von normalem C-Code aus verwendet werden. Nachdem die Hauptschleife über die Quellen beendet ist, müssen wir die akkumulierte Intensität zurücklesen und eine ganze Menge Puffer freigeben.

Eine Blockgröße von 256 funktionierte ganz gut, aber das ist einstellbar. Es ist eine gute Idee, Teile des Codes wie diesen in eine eigene .cu-Datei zu packen und mit nvcc zu kompilieren. Dies ermöglichte es uns, den restlichen Standard-C-Code mit den üblichen Optimierungen mit GCC oder Clang zu kompilieren. Und es ist kein Problem, das Ganze dann in eine .mex-Datei zu verlinken und mit Octave zu verwenden.

Die CUDA-API bietet viele Möglichkeiten, diese Simulation anzugehen, um eine bessere GPU-Auslastung zu erreichen, wie z. B. Streams, aber selbst mit dieser einfachen Implementierung können wir eine GPU-Auslastung von 70% für große Simulationsdomänen (z. B. 2048×2048 Abtastpunkte) erreichen. Die Implementierung einer Alternative zu FFTW für eine einzelne Transformation ist jedoch wenig sinnvoll. Der Overhead macht den Vorteil der Berechnung auf dem Grafikprozessor fast vollständig zunichte.

Lohnt sich der Aufwand?

Also, war es das wert? Meine Haupt-GPU ist ein Ryzen 7 8700G. Der Code, den ich normalerweise ausführe, wurde schon seit einiger Zeit für diese Art von CPU optimiert, um das meiste aus den 8 Kernen herauszuholen. Die GPU, mit der ich sie vergleichen wollte, ist eine RTX 3070Ti, kein Geschwindigkeitsdämon, selbst nach den Standards der Tage, als sie veröffentlicht wurde. Dennoch zeigen die Ergebnisse, dass sie teilweise kohärente Bilder etwa 10 Mal schneller berechnet als diese etwas neuere CPU. Das ist ein netter Geschwindigkeitszuwachs für einen ersten Versuch. Es ist nicht so einfach, eine 10-mal schnellere CPU zu einem beliebigen Preis zu finden.

Um dies mit angemessener Geschwindigkeit zu bewerkstelligen, muss man die Menge der pro Schleife übertragenen Daten in Grenzen halten. Da dies ein minimaler Aufwand war, der auf einem für eine CPU optimierten Code basierte, war ich eher zurückhaltend, zu viel Code zu konvertieren, damit er auf der GPU läuft. Glücklicherweise muss man nur den Teil der Transformation kopieren, der durch die Pupille geht, was eine sehr begrenzte Menge an zu kopierenden Daten darstellt. Vielleicht werde ich bei einer künftigen Überarbeitung mehr Code auf den Grafikprozessor verlagern, so dass alle Kopien pro Iteration entfallen. Es könnte ein weiterer Geschwindigkeitsgewinn mit 60%-70% möglich sein, und wenn alle kleinen Engpässe beseitigt sind, könnte eine 5070Ti wirklich glänzen.

Im Nachhinein betrachtet, ist die Berechnung von teilweise kohärenten Bildern perfekt für einen Grafikprozessor geeignet. Die große Anzahl von Transformationen, die durchgeführt werden müssen, können auf dem Grafikprozessor gespeichert werden, und die Akkumulation der Intensität ist wirklich einfach, und der Grafikprozessor leistet dabei fantastische Arbeit. Ich kann beim besten Willen keine Anwendung finden, die mit so wenig Aufwand so viele Verbesserungen bringt.