diff --git a/exercices/02-correction-gamma/correction/gamma.cu b/exercices/02-correction-gamma/correction/gamma.cu new file mode 100644 index 0000000000000000000000000000000000000000..3a86275bf38c832dc606f1d6f4edb1e4259dff1e --- /dev/null +++ b/exercices/02-correction-gamma/correction/gamma.cu @@ -0,0 +1,89 @@ +#include<stdio.h> + +__device__ unsigned int t_gamma(unsigned int pixel, double gamma) { + double pixel_norm = ((double) pixel)/255.0; + double new_pixel = 255.0*pow(pixel_norm, gamma); + return rint(new_pixel); +} +__global__ void alpha_ker(unsigned int* img, int nRows, int nCols, double gamma) { + int i = blockIdx.x * blockDim.x + threadIdx.x; + int j = blockIdx.y * blockDim.y + threadIdx.y; + int k = (i*nCols) + j; + if (i >= nRows || j >= nCols) + return; + img[k] = t_gamma(img[k], gamma); +} + +void print_ary(unsigned int* img, int nRows, int nCols) { + for (int i = 0; i < nRows; i++) { + for (int j = 0; j < nCols; j++) { + printf("%d, ", img[i*nCols + j]); + } + printf("\n"); + } +} + +void h_gamma(unsigned int* img, int nRows, int nCols, double gamma){ + for(int i=0; i<nRows; i++){ + for(int j=0; j<nCols; j++){ + double pixel_norm = ((double) img[i*nCols+j])/255.0; + double new_pixel = 255.0*pow(pixel_norm, gamma); + img[i*nCols+j] = rint(new_pixel); + } + } +} + +int main() { + double gamma = 10.0; + int NUM_ROWS = 20; + int NUM_COLS = 20; + int imgSize = NUM_ROWS*NUM_COLS*sizeof(unsigned int); + + unsigned int* h_img = (unsigned int*) malloc(imgSize); + unsigned int* h_img_res = (unsigned int*) malloc(imgSize); + unsigned int* d_img; + + //cudaError_t cuError = cudaSuccess; + srand(10); + for (int i = 0; i < NUM_ROWS; i++) { + for (int j = 0; j < NUM_COLS; j++) { + h_img[i*NUM_COLS + j] = ((double)rand()/RAND_MAX)*255; + } + } + memcpy(h_img_res, h_img, imgSize); + + print_ary(h_img_res, NUM_ROWS, NUM_COLS); + + cudaMalloc(&d_img, imgSize); + cudaMemcpy(d_img, h_img_res, imgSize, cudaMemcpyHostToDevice); + + //cuError = cudaMalloc((void**)&da, iSize); + //if (cudaSuccess != cuError) { + // printf("Failed to allocate memory\n"); + // return 1; + //} + //cuError = cudaMemcpy(da, ha, iSize, cudaMemcpyHostToDevice); + + dim3 dimGrid (2, 2, 1); + dim3 dimBlock (10, 10, 1); + + + alpha_ker<<<dimGrid, dimBlock>>>(d_img, NUM_ROWS, NUM_COLS, gamma); + + //cuError = cudaGetLastError(); + cudaMemcpy(h_img_res, d_img, imgSize, cudaMemcpyDeviceToHost); + cudaFree(d_img) ; + + printf("======================================\n"); + + print_ary(h_img_res, NUM_ROWS, NUM_COLS); + + h_gamma(h_img, NUM_ROWS, NUM_COLS, gamma); + + if(memcmp(h_img, h_img_res, imgSize)==0){ + printf("tout va bien\n"); + return 0; + } + printf("argl!\n"); + return 1; +} \ No newline at end of file diff --git a/exercices/03-equation-chaleur/README.md b/exercices/03-equation-chaleur/README.md new file mode 100644 index 0000000000000000000000000000000000000000..9a3fed8cca27d847629ba7579c886b6825aa9fb8 --- /dev/null +++ b/exercices/03-equation-chaleur/README.md @@ -0,0 +1,204 @@ +# TP 3: l'équation de la chaleur + +Le but de ce TP est de procéder à une implémentation simple d'un [code stencil](https://en.wikipedia.org/wiki/Stencil_code). +Ceci dans le but de résoudre l'équation de la chaleur sur un domaine carré en deux dimensions. +L'implémentation parallèle se fera avec CUDA C. + +L'algorithme décrit dans cet énoncé est admis sans démonstration. +Il n'est donc pas nécessaire de comprendre les mathématiques associées à la résolution de l'équation de la chaleur, ni de comprendre l'équation elle-même. +En effet, il faudrait une ou deux leçons complètes (c-à-d de deux à quatre heure) pour en introduire tous les concepts. +Néanmoins, la "traduction" de l'approche en algorithme devra être claire pour l'étudiant et son fonctionnement devrait être intuitif. + +Toutefois, les étudiants curieux peuvent se référer à: [équation de la chaleur](https://en.wikipedia.org/wiki/Heat_equation), ainsi qu'aux concepts suivants pour comprendre la description du problème: +- description du [flux](https://en.wikipedia.org/wiki/Flux), +- le [théorème de la divergence](https://en.wikipedia.org/wiki/Divergence_theorem), +- et l'[opérateur de Laplace](https://en.wikipedia.org/wiki/Laplace_operator), +- ou encore l'[équation de Laplace](https://en.wikipedia.org/wiki/Laplace%27s_equation). + +Ensuite pour la résolution numérique du problème: +- la [série de Taylor](https://en.wikipedia.org/wiki/Taylor_series), +- les [différences fines](https://en.wikipedia.org/wiki/Finite_difference). + +Les articles sont proposés en anglais, car je trouve que les articles Wikipedia anglophones couvrant les mathématiques sont plus claires que ceux écrits en français. +Les articles francophones étant souvent très formels sans tenter de donner une intuition sur ce qui est expliqué. +Il s'agit d'un avis personel, libre à vous de lire les articles dans la langue de votre choix. + +## Diffusion de la chaleur sur un plaque rectangulaire + +Le problème que nous cherchons à résoudre est le suivant: +- on a une plaque rectangulaire, +- on ne s'intéresse pas à la température initiale de la plaque, +- on applique une source de chaleur constante sur les bords de la plaque, +- la chaleur va se propager (diffuser) vers/depuis les bords de la plaque (du chaud vers le froid, principe de la thermodynamique), +- au bout d'un certain temps, la plaque va atteindre un état stable (c-à-d. plus de propagation, plus de changement de température au cours du temps). + +Ce phénomène se décrit bien mathématiquement par une équation aux dérivées partielles dont une solution analytique a été proposée par [Jean-Baptiste Joseph Fourier](https://en.wikipedia.org/wiki/Joseph_Fourier). +Mais dans notre cas, nous n'allons pas chercher à résoudre analytiquement l'équation de la chaleur, mais nous allons simplement appliquer un schéma numérique (un algorithme) pour approximer la solution. + +Tout d'abord, on ne peut pas considérer un domaine (c-à-d. notre plaque) continu. +Nous somme obligé de considérer un domaine discret, dénoté par $`U`$. +Ceci est inhérent à l'utilisation d'un ordinateur pour représenter le domaine. + +Le domaine s'apparente donc à une grille et nous sommes intéressés uniquement à approximer les valeurs de chaleur (la température) aux coordonnées $`(i,j)`$ rouges de la grille $`U`$: + + + +Les points aux coordonnées bleues sont les valeurs aux bords et ce sont ces points qui restent constants au court du temps. +Il n'est donc pas nécessaire de les "résoudre", vu qu'on les connait (ils sont constants). +Néanmoins ce sont eux qui vont définir la répartition finale de la température sur la plaque. + +Il est important de noter que les points de grilles (dénotés $`u_{ij}`$) sont équidistants les uns des autres. +Ceci même si la figure précédente ne l'illustre pas parfaitement. +Ce qui implique qu'on peut facilement représenter ce domaine à l'aide d'un tableau en deux dimensions composé de valeurs réelles. + +Une fois le domaine discrétisé et les valeurs aux bords imposées, la résolution se fait à l'aide d'une simple moyenne. +En effet, une solution au problème consiste à calculer l'évolution de la température d'un point $`u_{ij}`$ au temps $`t`$, comme étant la température moyenne de ses voisins au temps $`t-1`$. +Ce calcul s'applique du temps initial, noté $`t_0`$, jusqu'au temps où le système converge, noté $`t_n`$. +C'est-à-dire que les températures des tous les points $`u_{ij}`$ du domaine ne changent plus en appliquant la moyenne. + +Faire une telle opération se nomme une application de stencil. +Un stencil est un pattern ou motif fix qui décrit quels voisins sont utilisés pour calculer la valeur d'un point $`u_{ij}`$. +Il est appliqué en même temps à tous les points de l'espace considéré (application synchrone). + +Le stencil qui nous intéresse peut être schématisé comme: + + + +donc dans notre cas nous utiliseront les quatre voisins (points blancs) pour calculer la valeur d'un point du domaine (point noir, notre $`u_{ij}`$). +Plus formellement, nous calculons la moyenne des quatre voisins à chaque application du stencil au cours du temps et pour chaque point du domaine: + +$` u_{i,j}^{t+1} = \frac{1}{4} \cdot (u_{i-1,j}^{t} + u_{i+1,j}^{t} + u_{i,j-1}^{t} + u_{i,j+1}^{t} )`$ + +où $`u_{i,j}^{t}`$ est la valeur (la chaleur) à la coordonnée $`(i,j)`$ du domaine $`U`$ à l'itération $`t`$. + +Note pour les étudiants ayant suivi le cours de VISNUM: +- vous avez déjà appliqué un stencil, p.ex. avec le filtre moyenneur (passe-bas), +- il s'agit du **même** filtre que celui appliqué pour détecter les contours d'un objet. D'ailleur ce fitre s'appelait le filtre de Laplace. + +Pour atteindre l'état stationnaire, on applique itérativement le stencil sur tous les points $`(i,j)`$ du domaine et on produit successivement les domaines: $`U^0`$ (état inital, au temps $`t_0`$), $`U^1, U^2, \ldots, U^n`$ (au temps $`t_n`$). +Ceci jusqu'à ce que les valeurs convergent. +Numériquement cela implique que toutes les points points $`(i,j)`$ de $`U`$ ne changent plus suffisamment d'une itération à l'autre. + +Mais pour des questions de confort, spécialement pour les mesures de performance du code, nous utiliserons un nombre d'itération maximum $`n`$ à appliquer. + +### Pseudo-code séquentiel des boucles temporelle et spatiale + +Le pseudo-code suivant illustre l'algorithme approximant la solution de l'équation de Laplace pour un nombre d'itérations fixe `max_T`, pour un domaine discret composé de `max_I` $`\times`$ `max_J` points: + +```c +U = init() +for (int t = 0; t < max_T; t++) +{ + for (int i = 0; i < max_I; i++) + { + for (int j = 0; j < max_J; j++) + { + U_temp[i][j] = 0.25 * ( + U[i-1][j] + U[i+1][j] + + U[i][j-1] + U[i][j+1] + ) + } + } + U = U_temp +} +``` + +Pour l'état initiale (i.e. `init()`) vous pouvez choisir 0 comme valeur à l'intérieur du domaine et 1 aux bords du domaine: + +```c +U = 0 +for (int i = 0; i < max_I; i++) +{ + U[i][0] = 1.0 + U[i][max_J - 1] = 1.0 +} + +for (int j = 0; i < max_J; j++) +{ + U[0][j] = 1.0 + U[max_i - 1][j] = 1.0 +} +``` + +Nous rappelons que les valeurs au bord du domaine ne changent pas au cours du temps. + +Voici un exemple de taille $`7 \times 7`$ d'un état initial $`U^0`$: + +``` +1.0 1.0 1.0 1.0 1.0 1.0 1.0 +1.0 0.0 0.0 0.0 0.0 0.0 1.0 +1.0 0.0 0.0 0.0 0.0 0.0 1.0 +1.0 0.0 0.0 0.0 0.0 0.0 1.0 +1.0 0.0 0.0 0.0 0.0 0.0 1.0 +1.0 0.0 0.0 0.0 0.0 0.0 1.0 +1.0 1.0 1.0 1.0 1.0 1.0 1.0 +``` + +et après application d'une itération nous obtenons $`U^1`$: + +``` +1.0 1.0 1.0 1.0 1.0 1.0 1.0 +1.0 0.5 0.25 0.25 0.25 0.5 1.0 +1.0 0.25 0.0 0.0 0.0 0.25 1.0 +1.0 0.25 0.0 0.0 0.0 0.25 1.0 +1.0 0.25 0.0 0.0 0.0 0.25 1.0 +1.0 0.5 0.25 0.25 0.25 0.5 1.0 +1.0 1.0 1.0 1.0 1.0 1.0 1.0 +``` + + + + + +## Travail à réaliser + +Votre travail consiste à accélérer ce code stencil en utilisant l'architecture du GPU et le modèle de programmation [CUDA](https://en.wikipedia.org/wiki/CUDA). +Vous devrez également procéder à des mesures de performances. + +### Partie implémentation + +On rappelle le pseudo-code séquentiel de la boucle temporelle/spatiale. +Le pseudo-code suivant illustre l'algorithme approximant la solution de l'équation de Laplace pour un nombre d'itérations fixe `T`, pour un domaine discret composé de `X` $`\times`$ `Y` points: + +```c +U = init() // U a squared domain +for (int t = 0; t < T; t++) { + for (int i = 0; i < X; i++) { + for (int j = 0; j < Y; j++) { + U_temp[i][j] = 0.25 * ( + U[i-1][j] + U[i+1][j] + + U[i][j-1] + U[i][j+1] + ) + } + } + U = U_temp +} +``` + +## Synchronisation des threads CUDA + +On procède à un bref rappel sur la synchronisation des kernels et threads avec CUDA. + +Premièrement, le lancement d'un kernel sur le device vis-à-vis de l'hôte est une opération asynchrone. +Il faut donc procéder à une opération de synchronisation pour que l'hôte attende sur le device. +L'utilisation de [cudaDeviceSynchronize](https://docs.nvidia.com/cuda/cuda-runtime-api/group__CUDART__DEVICE.html#group__CUDART__DEVICE_1g10e20b05a95f638a4071a655503df25d) utilisé sur le host pour qu'il attendre la completion des kernels sur le device. +On rappelle que certaines opérations snychronisent implicitement le device et le host, comme [cudaMemcpy](https://docs.nvidia.com/cuda/cuda-runtime-api/group__CUDART__MEMORY.html#group__CUDART__MEMORY_1gc263dbe6574220cc776b45438fc351e8). + +Ensuite, les threads ne sont synchrones que par warp, mais il est possible de facilement les synchroniser sur un même block avec la fonction [__syncthreads](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#synchronization-functions). + +Finalement, les threads entre les blocks d'une grille sont **asynchrones**. Pour ce travail, il est nécessaire de synchroniser les threads entre chaque itération temporelle. + +On rappelle également que pour un stream donné, l'ordre des kernels soumis est préservé à l'exécution, on peut donc se baser dessus pour effectuer la synchronisation des threads. + + +Il est également possible de synchroniser les threads en utilisant les groupes coopératifs. Les étudiants curieux ou les CUDA enthusiasts peuvent se renseigner sur le sujet. + + +### Partie mesures + +Pour les mesures vous devrez simplement faire varier la taille du problème et mesurer son temps d'exécution. +Vous pouvez également jouer avec la taille de la grille. +Pensez à fixer le nombre d'itérations et ne pas utiliser de critère de convergence pour mettre fin à votre exécution. + + + diff --git a/exercices/03-equation-chaleur/fig/domain.png b/exercices/03-equation-chaleur/fig/domain.png new file mode 100644 index 0000000000000000000000000000000000000000..57072569b9d4f0acf44c590b88776cb8960ecd76 Binary files /dev/null and b/exercices/03-equation-chaleur/fig/domain.png differ diff --git a/exercices/03-equation-chaleur/fig/paradom.png b/exercices/03-equation-chaleur/fig/paradom.png new file mode 100644 index 0000000000000000000000000000000000000000..95778f8a9a2a05771d637fe27cc0395be24c4bfb Binary files /dev/null and b/exercices/03-equation-chaleur/fig/paradom.png differ diff --git a/exercices/03-equation-chaleur/fig/stencil.png b/exercices/03-equation-chaleur/fig/stencil.png new file mode 100644 index 0000000000000000000000000000000000000000..7099597bb901ed15efb23e12e552cfd954e6e786 Binary files /dev/null and b/exercices/03-equation-chaleur/fig/stencil.png differ