Introduzione

Per la maggior parte della mia carriera nell'ottica, ho simulato l'immagine di modulatori luminosi spaziali, principalmente SLM illuminati con sorgenti parzialmente coerenti. In genere, l'FFT svolge questo compito piuttosto rapidamente, almeno per un'unica immagine. A volte, però, bisogna effettuare molte simulazioni. Un esempio di questo è quando si lavora sui requisiti ottici di un sistema ottico. Un approccio che si è rivelato utile è quello di utilizzare Campionamento del cubo ipersferico in latino. Per un sistema ottico, probabilmente si includerebbero aberrazioni fino al 5° ordine, forse un po’ di più e questo diventa una simulazione piuttosto ingente. Forse 1000 o 2000 immagini. Quindi, forse si analizzano i dati. Forse lo spazio (gamma) analizzato era troppo piccolo o troppo grande e potrebbe essere necessario ripetere queste simulazioni un paio di volte.

Sebbene sia il computer a farlo, non è sempre facile programmare altre attività e, anche se si riesce a farlo, passare da una attività all’altra con intervalli di mezz’ora o un’ora (o forse anche di più) non è molto produttivo. Pertanto, per poter calcolare rapidamente FFT bidimensionali è necessario fare molta attenzione e essere molto concentrati.

La struttura di una simulazione dell'immagine

L'approccio alle simulazioni di immagini parzialmente coerenti dipende solitamente dal fatto che si intenda calcolare un gran numero di immagini utilizzando una specifica distribuzione di sorgenti; in questo caso l'approccio più appropriato potrebbe essere SOCS (Sum Of Coherent Systems), oppure se si desidera ricercare una buona distribuzione dell'intensità, l'overhead associato a SOCS non giustifica l'utilizzo di questo approccio e si consiglia di calcolare semplicemente le trasformazioni inverse delle ampiezze riflesse dalla maschera, oppure, come nel mio caso, SLM. Ogni di queste trasformazioni inverse rappresenta un'immagine la cui intensità deve essere sommata per poter determinare l'esposizione completa.

A seconda delle dimensioni del dominio periodico e del fattore di riempimento della pupilla dell’illuminazione (spesso chiamato sigma), potremmo facilmente finire per dover effettuare diverse centinaia di trasformazioni inverse. In questo caso, si può decidere di effettuare un campionamento più fine della sorgente luminosa, ma non si dovrebbe scendere significativamente al di sotto dei 100 punti di illuminazione. Se si scende a 25 (per esempio), ci saranno sicuramente delle variazioni visibili dovute alla periodicità della funzione di coerenza. Anche un valore di 50 è un numero accettabile; in alcuni casi la densità di campionamento viene fornita da un omogeneizzatore, quindi non c’è molto altro da aggiungere.

Per riassumere, l’immagine di base che parlerò qui consiste in una trasformazione in avanti per calcolare l’ampiezza che illumina la pupilla del sistema, e una serie di trasformazioni inverse che consistono semplicemente in traslazioni triviali di una trasformazione per calcolare le immagini parziali. Tra queste ci sono alcune operazioni di filtraggio della pupilla, ma questo è un procedimento molto economico computazionalmente.

API CUDA

In un primo momento, il API CUDA Sembra un po’ intimidante, e anche dopo aver passato un po’ di tempo a farlo, ero un po’ preoccupato su quanto dolore avrei dovuto sopportare prima di superare la fase di apprendimento. Ciò che complicava le cose era il fatto che volevo un codice che potesse essere utilizzato insieme a Octave o Matlab e che fosse scritto in C. Questo codice doveva restituire risultati identici anche quando veniva compilato per utenti che non avessero accesso a CUDA.

A mia grande sorpresa, l’accesso a CUDA fu molto più semplice di quanto avessi previsto, e i problemi che mi preoccupavano, come l’impossibilità di riavviare il sistema più volte a causa del blocco della scheda grafica, non si sono mai verificati.

Qualche scelta strategica

Il primo punto di decisione riguardava l’utilizzo della libreria CUDA FFTW-wrapper o l’utilizzo della libreria cuFFT direttamente. Se si dispone di un’applicazione standalone basata sulla libreria FFTW, allora potrebbe essere una scelta, anche se ci sono alcune cose ovvie che questa libreria non può fare, come ad esempio allocare sia i buffer per i dati sia quelli per il dispositivo e per il processore host. Nel linguaggio CUDA, la memoria “device” è presente sulla GPU e la memoria “host” è associata al processore host.

Con Octave/Matlab c’è un altro problema: l’utilizzo della libreria CUDA FFTW-wrapper può creare conflitti di nomi tra FFTW e le funzioni con lo stesso nome nella libreria wrapper. Anche se decidessimo di fornire implementazioni per tutte le varianti FFT di Octave/Matlab, non sapremmo se non ci sia un binario presente in qualche altro posto che stia aspettando di ricevere un programma da fftw_plan_dft_2d. Inoltre, sebbene la libreria wrapper fornisca una tale funzione, le librerie non sono compatibili con il formato binario. Pertanto, utilizzare la guida di conversione per cuFFT è una buona idea.

La soluzione più semplice

Perciò, iniziamo con qualcosa di più semplice. Qual è il minimo indispensabile per eseguire FFT sulla GPU? In realtà, non è molto. Bisogna

• Allocare la memoria della GPU – device_mem = cudaMalloc(n);
• Realizzare un piano – cufftPlan2D(&plan, n1, n2, CUFFT_Z2Z);
• Copiare i dati dell'host sulla GPU – cudaMemcpy(device_mem, host_mem, n, cudaMemcopyDefault);
• Eseguire il piano: cufftExecZ2Z(piano, dispositivo_mem_in, dispositivo_mem_out, CUFFT_FORWARD);
• Ritorna il risultato – cudaMemcpy(host_mem, device_mem, n, cudaMemcpyDefault);

Ovviamente sono piuttosto sciatto in questo ambito, ma procedere in questo modo ci permetterebbe di eseguire l’FFT sulla GPU, ma ciò sarebbe comunque un po’ una delusione. La maggior parte del tempo sarebbe speso per copiare i dati e formare il dispositivo. C’è una buona probabilità (se “buona” è la parola giusta anche per questo) che non ci accorgessimo che la GPU stesse effettivamente svolgendo qualsiasi attività, sebbene in realtà stesse calcolando il risultato corretto.

Così, cosa dovevamo fare con i dati che abbiamo copiato indietro dalla GPU? L'unica cosa che dovevamo fare era calcolare il prodotto della trasformazione e del suo coniugato per ogni campione e accumulare la somma finché non avessimo completato la scansione delle fonti di illuminazione. Se avessimo potuto farlo sulla GPU, almeno non avremmo dovuto rileggere il risultato più di una volta. Ciò significa che dobbiamo imparare a eseguire un certo codice banale sulla GPU.

Un codice che potrebbe farlo è:

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

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

Il prototipo di addIntensity2Df deve essere esportato in un header e può essere utilizzato da codice C normale. Dopo aver completato il ciclo principale sulle sorgenti, dovremo rileggere l'intensità accumulata e liberare un sacco di buffer.

Una dimensione del blocco di 256 è risultata abbastanza buona, ma questa caratteristica è regolabile. È una buona idea inserire in questo file .cu dei pezzi di codice come questo e compilarlo con nvcc. Questo ci ha permesso di compilare il resto del codice C standard con le consuete ottimizzazioni utilizzando GCC o Clang. E non ci sono problemi nel collegare tutto a un file .mex e usarlo con Octave.

L'API CUDA offre molte modalità per affrontare questa simulazione al fine di ottenere una migliore utilizzazione della GPU, come ad esempio i flussi, ma anche con questa semplice implementazione possiamo ottenere un'utilizzazione del 70% della GPU per grandi domini di simulazione (ad esempio 2048×2048 punti di campionamento). Tuttavia, non si ottiene molto implementando un'alternativa a FFTW per un singolo trasformato. L'overhead consuma quasi completamente il vantaggio derivante dal fatto di eseguire il calcolo sulla GPU.

Vale la pena l'impegno?

Quindi, valeva la pena? La mia GPU principale è un Ryzen 7 8700G. Il codice che utilizzo solitamente è stato ottimizzato per questo tipo di CPU da un po’ di tempo per sfruttare al meglio i 8 core. La GPU con cui volevo confrontarlo è un RTX 3070Ti, non un processore da gaming anche considerando gli standard dei tempi in cui è stato rilasciato. Tuttavia, i risultati dimostrano che calcola immagini parzialmente coerenti circa 10 volte più velocemente rispetto a questo processore relativamente recente. È un bel miglioramento rispetto al primo tentativo. Non è così facile trovare un processore 10 volte più veloce ad un prezzo ragionevole.

Per mantenere questa velocità elevata, è assolutamente necessario ridurre la quantità di dati trasferiti per ciclo. Poiché questo era un intervento minimo basato su un codice ottimizzato per una CPU, ero piuttosto riluttante a convertire troppo codice per essere eseguito sulla GPU. Fortunatamente, è sufficiente copiare la parte del trasformato che attraversa la pupilla, che è una quantità di dati molto limitata da copiare. Forse in una futura revisione potrò spostare più codice sulla GPU e quindi tutte le copie per iterazione verranno eliminate. Potrebbe essere possibile ottenere un ulteriore aumento della velocità di 60%-70%, e una volta eliminati tutti i piccoli ostacoli, un 5070Ti potrebbe davvero rivelarsi eccezionale.

Guardandoci indietro, la computazione di immagini parzialmente coerenti è perfetta per una GPU. Il gran numero di trasformazioni da effettuare può essere gestito dalla GPU e l'accumularsi dell'intensità è davvero semplice, e la GPU svolge un lavoro fantastico in questo senso. Dal mio punto di vista, non riesco a trovare un'applicazione che apporti così tanti miglioramenti con così poco sforzo.