Labo "Produit Scalaire".

[GPU.git] / WCudaMSE / Student_Cuda / src / cpp / core / 02_ProduitScalaire / ProduitScalaire.cu

diff --git a/WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu b/WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu
index c96562a..27c3a55 100644 (file)
--- a/WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu
+++ b/WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu
@@ -1,16 +1,20 @@
+#include "ProduitScalaire.h"
+
 #include <iostream>
 #include <cmath>
+#include <stdio.h>
 using namespace std;
 
 #include "Indice1D.h"
 #include "cudaTools.h"
 #include "Device.h"
+#include "Lock.h"
 
 #define M_V 200
 #define M_W 200
 
 #define VI 1.4422495703074083
-#define WI 0.7390850782394409
+#define WI 0.739085133215160672293109200837
 
 /*
  * Renvoie la valeur du ième élément du vecteur v.
@@ -24,7 +28,13 @@ double v(long i)
         const double xCarre = x * x;
         x = x - (xCarre * x - 3) / (3 * xCarre);
         }
-    return (x / VI) * sqrt(double(i));
+
+    /* Debug afin d'ajuster VI
+    if (Indice1D::tid() == 0)
+        printf("x: %.30f, VI: %.30f, x / VI: %.30f\n", x, VI, x / VI);
+        */
+
+    return (x / VI) * sqrt(double(i)); // x / VI doit être égal à 1.
     }
 
 /*
@@ -36,13 +46,20 @@ double w(long i)
     double x = abs(cos(double(i)));
     for (long j = 1; j <= M_W; j++)
         x = x - (cos(x) - x) / (-sin(x) - 1);
+
+    /* Debug afin d'ajuster WI
+    if (Indice1D::tid() == 0)
+        printf("x: %.30f, WI: %.30f, x / WI: %.30f\n", x, WI, x / WI); */
+
     return (x / WI) * sqrt(double(i));
     }
 
 /*
+ * 1) Chaque thread calcule un résultat intermediaire qu'il va ensuite placer en shared memory.
  * n: La taille des deux vecteurs.
  */
-__device__ void reductionIntraThread(int n, double* tabResultSM)
+__device__
+void reductionIntraThread(int n, double* tabSM)
     {
     const int NB_THREAD = Indice1D::nbThread();
     const int TID = Indice1D::tid();
@@ -55,44 +72,63 @@ __device__ void reductionIntraThread(int n, double* tabResultSM)
         threadResult += v(s) * w(s);
         s += NB_THREAD;
         }
-    tabResultSM[TID_LOCAL] = threadResult;
+
+    tabSM[TID_LOCAL] = threadResult;
     }
 
 /*
- * Combine les résultats de 'tabResultSM' dans 'tabResulSM[0]'
+ * Combine les résultats de 'tabSM' dans 'tabSM[0]'
  */
-__device__ void combine(double* tabResultSM, int middle)
+__device__
+void combine(double* tabSM, int middle)
     {
-        const int TID_LOCAL = Indice1D::tidLocal();
-        const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
-
-        int s = TID_LOCAL;
-        while (s < middle)
-            {
-            tabResultSM[s] = tabResultSM[s] + tabResultSM[s + middle];
-            s += NB_THREAD_LOCAL;
-            }
+    const int TID_LOCAL = Indice1D::tidLocal();
+    const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
+
+    int s = TID_LOCAL;
+    while (s < middle)
+        {
+        tabSM[s] = tabSM[s] + tabSM[s + middle];
+        s += NB_THREAD_LOCAL;
+        }
     }
 
-__device__ void reductionIntraBlock(double* tabResultSM)
+/*
+ * 2) La shared memory est réduite, le résultat est placé dans 'tabSM[0]'.
+ */
+__device__
+void reductionIntraBlock(double* tabSM)
     {
     const int TAB_SIZE = blockDim.x;
     int middle = TAB_SIZE / 2;
 
     while (middle > 0)
         {
-        combine(tabResultSM, middle);
+        combine(tabSM, middle);
         middle /= 2;
         __syncthreads();
         }
     }
 
-__device__ void reductionInterBlock(double* tabResultSM, float* ptrResult)
+__device__
+int mutexReductionInterBlock = 0;
+
+/*
+ * 3) Le 'tabSM[0]' de chaque bloc est reduit dans 'ptrResult' qui se trouve en global memory.
+ */
+__device__
+void reductionInterBlock(double* tabSM, double* ptrResult)
     {
     const int TID_LOCAL = Indice1D::tidLocal();
     if (TID_LOCAL == 0)
         {
-        atomicAdd(ptrResult, float(tabResultSM[0]));
+        Lock lock(&mutexReductionInterBlock);
+        lock.lock();
+        (*ptrResult) += tabSM[0];
+        lock.unlock();
+
+        // Si on travail en float (pas besoin de mutex) :
+        // atomicAdd(ptrResult, float(tabSM[0]));
         }
     }
 
@@ -103,36 +139,36 @@ __device__ void reductionInterBlock(double* tabResultSM, float* ptrResult)
  * ptrResult: Le resultat du produit scalaire.
  */
 __global__
-void produitScalaire(int n, float* ptrResult)
+void produitScalaire(int n, double* ptrResult)
     {
-    extern __shared__ double tabResultSM[]; // Shared memory.
+    extern __shared__ double tabSM[]; // Dynamic shared memory.;
 
     // 1) Réduction intra-thread.
-    reductionIntraThread(n, tabResultSM);
+    reductionIntraThread(n, tabSM);
 
     __syncthreads();
 
     // 2) Réduction intra-block.
-    reductionIntraBlock(tabResultSM);
+    reductionIntraBlock(tabSM);
 
     // 3) Réduction inter-block.
-    reductionInterBlock(tabResultSM, ptrResult);
+    reductionInterBlock(tabSM, ptrResult);
     }
 
 double resultatTheorique(long n)
-{
+    {
     n -= 1;
-    return (n / 2.0) * (n+1);
-}
+    return (n / 2.0) * (n + 1);
+    }
 
 bool produitScalaire()
     {
-    const int N = 100000000; // Taille des deux vecteurs.
+    const int N = 10000000; // Taille des deux vecteurs : 10 * 10^6.
 
     // Allocation coté GPU en global memory (GM).
-    float* ptrDevResult = 0;
-    HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(float)));
-    HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(float)));
+    double* ptrDevResult = 0;
+    HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(double)));
+    HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(double)));
 
     // Paramètre de l'appel de la fonction sur le device.
     const dim3 dg(256, 1, 1);
@@ -142,13 +178,15 @@ bool produitScalaire()
 
     produitScalaire<<<dg, db, SMSize>>>(N, ptrDevResult);
 
-    float res;
+    cudaDeviceSynchronize(); // Utilisé pour flusher les prints sur le stdout à partir du device.
+
+    double res;
     // Barrière implicite de synchronisation ('cudaMemCpy').
-    HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(float), cudaMemcpyDeviceToHost));
+    HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(double), cudaMemcpyDeviceToHost));
 
-    double resTheo = resultatTheorique(N);
+    const double resTheo = resultatTheorique(N);
 
-    cout.precision(10);
+    cout.precision(30);
     cout << "Résultat : " << res << endl;
     cout << "Résultat théorique : " << resTheo << endl;
     cout << "Différence absolue : " << resTheo - res << endl;