ajout correction TP2 et ajout TP3

exercices/02-correction-gamma/correction/gamma.cu

+#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

exercices/03-equation-chaleur/README.md

+# 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`$:
+
+![dom](fig/domain.png)
+
+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:
+
+![sten](fig/stencil.png)
+
+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.
+
+
+