Optique physique sur le GPU Nvidia

jarek

Intro

Pendant la majeure partie de ma carrière en optique, j'ai simulé l'imagerie de modulateurs spatiaux de lumière, et plus particulièrement de SLMs illuminés par des sources partiellement cohérentes. En général, la FFT fait un travail assez rapide, au moins pour une seule image. Parfois, cependant, il est nécessaire d'effectuer un grand nombre de simulations. C'est le cas par exemple lorsque l'on travaille sur les exigences optiques d'un système optique. Une approche qui s'est avérée utile est l'utilisation de Échantillonnage hypercube latin. Pour un système optique, on inclurait probablement jusqu'aux aberrations de cinquième ordre, voire un peu plus, et cela deviendrait une simulation assez importante. Peut-être 1000 ou 2000 images. Ensuite, vous pouvez examiner les données. Il se peut que l'espace sondé (la portée) soit trop petit ou trop grand, et que l'on veuille refaire ces simulations plusieurs fois.

Bien que ce soit l'ordinateur qui le fasse, il n'est pas toujours facile de programmer une autre tâche, et même si c'est le cas, il n'est pas très productif de changer de tâche toutes les demi-heures ou toutes les heures (voire plus). C'est pourquoi il est important de pouvoir calculer rapidement des FFT bidimensionnelles, et c'est même très important.

La structure d'une simulation d'image

L'approche des simulations d'images partiellement cohérentes dépend généralement du fait que vous prévoyez de calculer un grand nombre d'images en utilisant une distribution spécifique de sources, alors votre approche pourrait être SOCS (Sum Of Coherent Systems), ou si vous voulez rechercher une bonne distribution d'intensité, le surcoût de SOCS n'en vaut pas la peine et l'approche consistera simplement à calculer les transformées inverses des amplitudes réfléchies par votre masque ou, comme c'était généralement le cas pour moi, SLM. Chacune de ces transformées inverses représente une image dont l'intensité doit être additionnée afin de trouver l'exposition complète.

En fonction de la taille du domaine périodique et du facteur de remplissage de la pupille de l'illumination (généralement appelé sigma), nous pourrions facilement nous retrouver à devoir effectuer plusieurs centaines de transformations inverses. À ce stade, on peut décider de procéder à un échantillonnage plus grossier de la source lumineuse, mais je ne descendrais pas en dessous de 100 points d'illumination. Si nous descendons à 25 (par exemple), il y aura certainement des contributions visibles de la périodicité de la fonction de cohérence. Même 50 est un peu juste, alors en règle générale, 100 n'est pas un mauvais chiffre. Parfois, la densité d'échantillonnage est donnée par un homogénéisateur, il n'y a donc pas grand-chose à ajouter.

Pour résumer, la simulation d'image de base dont je parlerai ici consiste en une transformation directe pour calculer l'amplitude qui illumine la pupille du système, et une série de transformations inverses qui ne sont que des translations triviales d'une transformation, pour calculer les images partielles. Il y a un filtrage de la pupille entre les deux, mais il s'agit d'une opération très peu coûteuse sur le plan informatique.

API CUDA

Dans un premier temps, le API CUDA peut sembler un peu intimidant, et même avec le temps considérable que j'ai passé à faire ces choses, j'étais un peu inquiet de la douleur que j'allais devoir endurer avant de sortir de la courbe d'apprentissage. Ce qui a quelque peu compliqué les choses, c'est que je voulais un code qui puisse être utilisé avec Octave ou Matlab et je le voulais en C. Ce code doit également renvoyer des résultats identiques lorsqu'il est compilé pour un utilisateur qui n'a pas accès à CUDA.

À ma grande surprise, CUDA a été beaucoup plus facile à appréhender que je ne l'avais prévu, et les choses que je craignais, comme devoir redémarrer un nombre infini de fois parce que j'avais réussi à bloquer ma carte graphique, ne se sont jamais produites.

Quelques choix stratégiques

Le premier point de décision concernant l'API était de savoir si je devais utiliser la bibliothèque CUDA FFTW-wrapper ou faire un peu marche arrière et utiliser directement la bibliothèque cuFFT. Si vous avez une application autonome basée sur la bibliothèque FFTW, c'est peut-être un choix, même s'il y a des choses évidentes que ce wrapper ne peut pas faire, comme allouer des tampons à la fois au périphérique et à l'hôte pour vos données. Dans le jargon CUDA, la mémoire "device" se trouve sur le GPU et la mémoire "host" est liée à votre CPU.

Avec Octave/Matlab, il y a une chose de plus, c'est que l'utilisation de la bibliothèque CUDA FFTW-wrapper créera des conflits de noms entre FFTW et les fonctions de même nom dans la bibliothèque wrapper. Même si nous décidons de fournir des implémentations pour toutes les variantes FFT d'Octave/Matlab, nous ne saurons toujours pas s'il n'y a pas un binaire quelque part qui s'attend à obtenir un plan de fftw_plan_dft_2d, et bien que la bibliothèque wrapper fournisse une telle fonction, les bibliothèques ne sont pas compatibles au niveau binaire. Par conséquent, l'utilisation du guide de conversion vers cuFFT est une bonne idée.

La solution la plus simple

Commençons donc par le bas. Quel est le strict minimum à faire pour exécuter des FFT sur le GPU ? En fait, il n'y a pas grand-chose à faire. Nous devons

Allocation de la mémoire du GPU - device_mem = cudaMalloc(n) ;
Faire un plan - cufftPlan2D(&plan, n1, n2, CUFFT_Z2Z) ;
Copie des données de l'hôte vers le GPU - cudaMemcpy(device_mem, host_mem, n, cudaMemcopyDefault) ;
Exécuter le plan - cufftExecZ2Z(plan, device_mem_in, device_mem_out, CUFFT_FORWARD) ;
Ramener le résultat - cudaMemcpy(host_mem, device_mem, n, cudaMemcpyDefault) ;

Je suis évidemment un peu négligent, mais en suivant ces lignes, nous pourrons effectivement effectuer la FFT sur le GPU, mais ce sera aussi un peu décevant. La plupart du temps serait consacré à copier des données vers et depuis le périphérique. Il y a de bonnes chances (si le mot "bonnes" est approprié) que nous ne remarquions pas que le GPU fait quoi que ce soit, si ce n'est qu'il calcule le bon résultat.

Alors, que devons-nous faire avec les données que nous avons copiées depuis le GPU ? La seule chose que nous devions faire était de calculer le produit de la transformée et son conjugué pour chaque échantillon et d'accumuler la somme jusqu'à ce que nous ayons fini de scanner les sources d'illumination. Si nous pouvions le faire sur le GPU, nous n'aurions pas à relire le résultat plus d'une fois. Cela signifie que nous devons apprendre à exécuter un code trivial sur le GPU.

Voici un morceau de code qui permettrait de faire cela :


void __global__ _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) * facteur ;
    }
}

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) ;
}

Le prototype de addIntensity2Df doit être exporté dans un en-tête et peut être utilisé à partir d'un code C normal. Une fois la boucle principale sur les sources terminée, nous devrons relire l'intensité accumulée et libérer un grand nombre de tampons.

Une taille de bloc de 256 a très bien fonctionné, mais elle peut être ajustée. C'est une bonne idée de mettre des morceaux de code comme celui-ci dans son propre fichier .cu et de le compiler avec nvcc. Cela nous a permis de compiler le reste du code C standard avec ses optimisations habituelles en utilisant GCC ou Clang. Et il n'y a aucun problème pour lier le tout à un fichier .mex et l'utiliser avec Octave.

L'API CUDA offre de nombreuses façons d'aborder cette simulation afin d'obtenir une meilleure utilisation du GPU, comme par exemple les flux, mais même avec cette implémentation simple, nous pouvons obtenir une utilisation du GPU de 70% pour de grands domaines de simulation (par exemple 2048×2048 points d'échantillonnage). Il y a cependant peu à gagner à mettre en œuvre une alternative à la FFTW pour une seule transformation. Les frais généraux consomment presque entièrement le gain apporté par l'exécution du calcul sur le GPU.

Le jeu en vaut-il la chandelle ?

Alors, est-ce que cela en valait la peine ? Mon GPU principal est un Ryzen 7 8700G. Le code que j'exécute habituellement a été optimisé pour ce type de CPU depuis un certain temps afin de tirer le meilleur parti des 8 cœurs. Le GPU avec lequel je voulais le comparer est une RTX 3070Ti, qui n'est pas un démon de vitesse, même selon les standards de l'époque où elle est sortie. Néanmoins, les résultats montrent qu'il calcule des images partiellement cohérentes environ 10 fois plus vite que ce CPU assez récent. C'est une belle accélération pour une première tentative. Il n'est pas facile de trouver un processeur 10 fois plus rapide à n'importe quel prix.

Pour que cela fonctionne à une vitesse décente, il faut réduire la quantité de données transférées par boucle, ce qui est obsolète. Comme il s'agissait d'un effort minimal basé sur un code optimisé pour un CPU, j'étais plutôt réticent à convertir trop de code pour l'exécuter sur le GPU. Heureusement, il suffit de copier la partie de la transformation qui passe par la pupille, ce qui représente une quantité très limitée de données à copier. Peut-être que lors d'une prochaine révision, je transférerai une plus grande partie du code sur le GPU, ce qui permettra d'éliminer toutes les copies par itération. Il pourrait y avoir un autre gain de vitesse de 60%-70% et une fois que tous les petits goulots d'étranglement seront éliminés, un 5070Ti pourrait vraiment briller.

En y repensant, le calcul d'images partiellement cohérentes convient parfaitement à un GPU. Le grand nombre de transformations à effectuer peut être conservé sur le GPU et l'accumulation de l'intensité est vraiment simple, et le GPU fait un travail fantastique. De mémoire, je n'ai pas trouvé d'application qui apporte autant d'améliorations avec si peu d'efforts.