Ajout du TP "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
new file mode 100644 (file)
index 0000000..c96562a
--- /dev/null
+++ b/WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu
@@ -0,0 +1,158 @@
+#include <iostream>
+#include <cmath>
+using namespace std;
+
+#include "Indice1D.h"
+#include "cudaTools.h"
+#include "Device.h"
+
+#define M_V 200
+#define M_W 200
+
+#define VI 1.4422495703074083
+#define WI 0.7390850782394409
+
+/*
+ * Renvoie la valeur du ième élément du vecteur v.
+ */
+__device__
+double v(long i)
+    {
+    double x = 1.5 + abs(cos(double(i)));
+    for (long j = 1; j <= M_V; j++)
+        {
+        const double xCarre = x * x;
+        x = x - (xCarre * x - 3) / (3 * xCarre);
+        }
+    return (x / VI) * sqrt(double(i));
+    }
+
+/*
+ * Renvoie la valeur du ième élément du vecteur w.
+ */
+__device__
+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);
+    return (x / WI) * sqrt(double(i));
+    }
+
+/*
+ * n: La taille des deux vecteurs.
+ */
+__device__ void reductionIntraThread(int n, double* tabResultSM)
+    {
+    const int NB_THREAD = Indice1D::nbThread();
+    const int TID = Indice1D::tid();
+    const int TID_LOCAL = Indice1D::tidLocal();
+
+    double threadResult = 0.0;
+    int s = TID;
+    while (s < n)
+        {
+        threadResult += v(s) * w(s);
+        s += NB_THREAD;
+        }
+    tabResultSM[TID_LOCAL] = threadResult;
+    }
+
+/*
+ * Combine les résultats de 'tabResultSM' dans 'tabResulSM[0]'
+ */
+__device__ void combine(double* tabResultSM, 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;
+            }
+    }
+
+__device__ void reductionIntraBlock(double* tabResultSM)
+    {
+    const int TAB_SIZE = blockDim.x;
+    int middle = TAB_SIZE / 2;
+
+    while (middle > 0)
+        {
+        combine(tabResultSM, middle);
+        middle /= 2;
+        __syncthreads();
+        }
+    }
+
+__device__ void reductionInterBlock(double* tabResultSM, float* ptrResult)
+    {
+    const int TID_LOCAL = Indice1D::tidLocal();
+    if (TID_LOCAL == 0)
+        {
+        atomicAdd(ptrResult, float(tabResultSM[0]));
+        }
+    }
+
+/**
+ * La taille de la shared memory (en terme de # de sizeof(double)) doit
+ * être égal à la taille des blocs.
+ * n: La taille des deux vecteurs.
+ * ptrResult: Le resultat du produit scalaire.
+ */
+__global__
+void produitScalaire(int n, float* ptrResult)
+    {
+    extern __shared__ double tabResultSM[]; // Shared memory.
+
+    // 1) Réduction intra-thread.
+    reductionIntraThread(n, tabResultSM);
+
+    __syncthreads();
+
+    // 2) Réduction intra-block.
+    reductionIntraBlock(tabResultSM);
+
+    // 3) Réduction inter-block.
+    reductionInterBlock(tabResultSM, ptrResult);
+    }
+
+double resultatTheorique(long n)
+{
+    n -= 1;
+    return (n / 2.0) * (n+1);
+}
+
+bool produitScalaire()
+    {
+    const int N = 100000000; // Taille des deux vecteurs.
+
+    // 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(double); // 256 doubles;
+
+    produitScalaire<<<dg, db, SMSize>>>(N, ptrDevResult);
+
+    float res;
+    // Barrière implicite de synchronisation ('cudaMemCpy').
+    HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(float), cudaMemcpyDeviceToHost));
+
+    double resTheo = resultatTheorique(N);
+
+    cout.precision(10);
+    cout << "Résultat : " << res << endl;
+    cout << "Résultat théorique : " << resTheo << endl;
+    cout << "Différence absolue : " << resTheo - res << endl;
+    cout << "Différence relatif : " << 100 * (resTheo - res) / (resTheo + res) << " %" << endl;
+
+    return true;
+    }