Ajout des cas de tests pour les TP non-graphiques. Cleanage en tous genres.

[GPU.git] / WCudaMSE / Student_Cuda / src / cpp / core / 04_MonteCarlo / MonteCarlo.cu

diff --git a/WCudaMSE/Student_Cuda/src/cpp/core/04_MonteCarlo/MonteCarlo.cu b/WCudaMSE/Student_Cuda/src/cpp/core/04_MonteCarlo/MonteCarlo.cu
new file mode 100644 (file)
index 0000000..c51910e
--- /dev/null
+++ b/WCudaMSE/Student_Cuda/src/cpp/core/04_MonteCarlo/MonteCarlo.cu
@@ -0,0 +1,176 @@
+#include "MonteCarlo.h"
+
+#include <iostream>
+#include <cmath>
+#include <stdio.h>
+using namespace std;
+
+#include <curand_kernel.h>
+
+#include "Indice1D.h"
+#include "cudaTools.h"
+#include "Device.h"
+#include "Lock.h"
+
+// Paramètres pour le calcul de pi.
+const float X_MIN = 0;
+const float X_MAX = 1;
+const float M = 1;
+
+/*
+ * 1) Chaque thread calcule un résultat intermediaire qu'il va ensuite placer en shared memory.
+ * n: Nombre d'échantillon.
+ */
+__device__
+void mc_reductionIntraThread(int n, uint* tabSM, curandState* tabGeneratorThread)
+    {
+    const int NB_THREAD = Indice1D::nbThread();
+    const int TID = Indice1D::tid();
+    const int TID_LOCAL = Indice1D::tidLocal();
+
+    curandState& localState = tabGeneratorThread[TID];
+
+    uint threadResult = 0;
+    int s = TID;
+    while (s < n)
+        {
+        // Tire un point aléatoire entre (0,0) et (1,1), voir X_MIN, X_MAX et M.
+        const float x = curand_uniform(&localState);
+        const float y = curand_uniform(&localState);
+        const float fx = sqrtf(1 - x * x);
+
+        // Si y est sous la courbe alors on le comptabilise.
+        if (y < fx)
+            threadResult += 1;
+
+        s += NB_THREAD;
+        }
+
+    tabSM[TID_LOCAL] = threadResult;
+    }
+
+/*
+ * Combine les résultats de 'tabSM' dans 'tabSM[0]'
+ */
+__device__
+void mc_combine(uint* 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 mc_reductionIntraBlock(uint* tabSM)
+    {
+    const int TAB_SIZE = blockDim.x;
+    int middle = TAB_SIZE / 2;
+
+    while (middle > 0)
+        {
+        mc_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 mc_reductionInterBlock(uint* tabSM, uint* ptrResult)
+    {
+    const int TID_LOCAL = Indice1D::tidLocal();
+    if (TID_LOCAL == 0)
+        {
+        atomicAdd(ptrResult, tabSM[0]);
+        }
+    }
+
+/**
+ * La taille de la shared memory (en terme de # de sizeof(uint)) doit
+ * être égal à la taille des blocs.
+ * n: le nombre de point à tiré aléatoirement.
+ * ptrResult: Le resultat du calcul.
+ * tabGeneratorThread: les generateurs aléatoires
+ */
+__global__
+void monteCarlo(int n,  uint* ptrResult, curandState* tabGeneratorThread)
+    {
+    extern __shared__ uint tabSM[]; // Dynamic shared memory.
+
+    // 1) Réduction intra-thread.
+    mc_reductionIntraThread(n, tabSM, tabGeneratorThread);
+
+    __syncthreads();
+
+    // 2) Réduction intra-block.
+    mc_reductionIntraBlock(tabSM);
+
+    // 3) Réduction inter-block.
+    mc_reductionInterBlock(tabSM, ptrResult);
+    }
+
+__global__
+void setupKernelRand(curandState* tabGeneratorThread, int deviceId)
+    {
+    int tid = Indice1D::tid();
+    int deltaSeed = deviceId * INT_MAX;
+    int deltaSequence = deviceId * 100;
+    int deltaOffset = deviceId * 100;
+
+    int seed = 1234 + deltaSeed;
+    int sequenceNumber = tid + deltaSequence;
+    int offset = deltaOffset;
+    curand_init(seed, sequenceNumber, offset, &tabGeneratorThread[tid]);
+    }
+
+bool monteCarlo()
+    {
+    cout << endl << "monteCarlo() ..." << endl;
+
+    // Nombre de point total tiré aléatoirement.
+    const int N = 1000000;
+
+    // Allocation coté GPU en global memory (GM).
+    uint* ptrDevResult;
+    HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(uint)));
+    HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(uint)));
+
+    // 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(uint); // 256 uint;
+
+    // Gestion des générateurs aléatoires (un par thread).
+    curandState* tabGeneratorThread;
+    HANDLE_ERROR(cudaMalloc(&tabGeneratorThread, db.x * dg.x * sizeof(curandState)));
+    setupKernelRand<<<dg, db>>>(tabGeneratorThread, 0); // Device 0 (à utilise si utilisation en MP).
+    cudaDeviceSynchronize();
+
+    monteCarlo<<<dg, db, SMSize>>>(N, ptrDevResult, tabGeneratorThread);
+
+    // cudaDeviceSynchronize(); // Utilisé pour flusher les prints sur le stdout à partir du device (debug).
+
+    uint nb;
+    // Barrière implicite de synchronisation ('cudaMemCpy').
+    HANDLE_ERROR(cudaMemcpy(&nb, ptrDevResult, sizeof(uint), cudaMemcpyDeviceToHost));
+
+    double integrale = (double)nb / (double)N * (X_MAX - X_MIN) * M;
+    double pi = 4 * integrale;
+
+    cout.precision(20);
+    cout << "Approximation de PI : " << pi << endl;
+
+    return true;
+    }