--- /dev/null
+#include "Saucisson.h"
+
+#include <iostream>
+#include <cmath>
+#include <stdio.h>
+using namespace std;
+
+#include "Indice1D.h"
+#include "cudaTools.h"
+#include "Device.h"
+#include "Lock.h"
+
+/*
+ * 1) Chaque thread calcule un résultat intermediaire qu'il va ensuite placer en shared memory.
+ * n: Nombre d'échantillon.
+ */
+__device__
+void reductionIntraThread(int n, float deltaX, float* tabSM)
+ {
+ const int NB_THREAD = Indice1D::nbThread();
+ const int TID = Indice1D::tid();
+ const int TID_LOCAL = Indice1D::tidLocal();
+
+ float threadResult = 0.0;
+ int s = TID;
+ while (s < n)
+ {
+ const float i = s + 1;
+ const float xi = -1 + i * deltaX;
+ threadResult += sqrtf(1-xi*xi);
+
+ s += NB_THREAD;
+ }
+
+ tabSM[TID_LOCAL] = threadResult;
+ }
+
+/*
+ * Combine les résultats de 'tabSM' dans 'tabSM[0]'
+ */
+__device__
+void combine(float* tabSM, int middle)
+ {
+ const int TID_LOCAL = Indice1D::tidLocal();
+ const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
+
+ int s = TID_LOCAL;
+ while (s < middle)
+ {
+ tabSM[s] += tabSM[s + middle];
+ s += NB_THREAD_LOCAL;
+ }
+ }
+
+/*
+ * 2) La shared memory est réduite, le résultat est placé dans 'tabSM[0]'.
+ */
+__device__
+void reductionIntraBlock(float* tabSM)
+ {
+ const int TAB_SIZE = blockDim.x;
+ int middle = TAB_SIZE / 2;
+
+ while (middle > 0)
+ {
+ combine(tabSM, middle);
+ middle /= 2;
+ __syncthreads(); // Synchronisation des threads au niveau du bloc.
+ }
+ }
+
+/*
+ * 3) Le 'tabSM[0]' de chaque bloc est reduit dans 'ptrResult' qui se trouve en global memory.
+ */
+__device__
+void reductionInterBlock(float* tabSM, float* ptrResult)
+ {
+ const int TID_LOCAL = Indice1D::tidLocal();
+ if (TID_LOCAL == 0)
+ {
+ atomicAdd(ptrResult, float(tabSM[0]));
+ }
+ }
+
+/**
+ * La taille de la shared memory (en terme de # de sizeof(float)) doit
+ * être égal à la taille des blocs.
+ * n: le nombre d'échantillon
+ * ptrResult: Le resultat du calcul de pi.
+ */
+__global__
+void saucisson(int n, float deltaX, float* ptrResult)
+ {
+ extern __shared__ float tabSM[]; // Dynamic shared memory.
+
+ // 1) Réduction intra-thread.
+ reductionIntraThread(n, deltaX, tabSM);
+
+ __syncthreads();
+
+ // 2) Réduction intra-block.
+ reductionIntraBlock(tabSM);
+
+ // 3) Réduction inter-block.
+ reductionInterBlock(tabSM, ptrResult);
+ }
+
+bool saucisson()
+ {
+ // Nombre d'échantillon. Au-delà, la qualité du résultat n'est pas meilleure. Il faudrait employé des doubles à la place de floats.
+ const int N = 100000;
+
+ // Allocation coté GPU en global memory (GM).
+ float* ptrDevResult = 0;
+ HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(float)));
+ HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(float)));
+
+ // Paramètre de l'appel de la fonction sur le device.
+ const dim3 dg(256, 1, 1);
+ const dim3 db(256, 1, 1);
+ Device::assertDim(dg, db);
+ const size_t SMSize = db.x * sizeof(float); // 256 floats;
+
+ const float deltaX = 2.0f / N;
+ saucisson<<<dg, db, SMSize>>>(N, deltaX, ptrDevResult);
+
+ // cudaDeviceSynchronize(); // Utilisé pour flusher les prints sur le stdout à partir du device (debug).
+
+ float pi;
+ // Barrière implicite de synchronisation ('cudaMemCpy').
+ HANDLE_ERROR(cudaMemcpy(&pi, ptrDevResult, sizeof(float), cudaMemcpyDeviceToHost));
+ pi *= 2 * deltaX;
+
+ cout.precision(20);
+ cout << "Approximation de PI : " << pi << endl;
+
+ return true;
+ }