WCudaMSE/Student_Cuda/src/cpp/core/02_ProduitScalaire/ProduitScalaire.cu

   1 #include "ProduitScalaire.h"
   2
   3 #include <iostream>
   4 #include <cmath>
   5 #include <stdio.h>
   6 using namespace std;
   7
   8 #include "Indice1D.h"
   9 #include "cudaTools.h"
  10 #include "Device.h"
  11 #include "Lock.h"
  12
  13 #define M_V 200
  14 #define M_W 200
  15
  16 #define VI 1.4422495703074083
  17 #define WI 0.739085133215160672293109200837
  18
  19 /*
  20  * Renvoie la valeur du ième élément du vecteur v.
  21  */
  22 __device__
  23 double v(long i)
  24     {
  25     double x = 1.5 + abs(cos(double(i)));
  26     for (long j = 1; j <= M_V; j++)
  27         {
  28         const double xCarre = x * x;
  29         x = x - (xCarre * x - 3) / (3 * xCarre);
  30         }
  31
  32     /* Debug afin d'ajuster VI
  33     if (Indice1D::tid() == 0)
  34         printf("x: %.30f, VI: %.30f, x / VI: %.30f\n", x, VI, x / VI);
  35         */
  36
  37     return (x / VI) * sqrt(double(i)); // x / VI doit être égal à 1.
  38     }
  39
  40 /*
  41  * Renvoie la valeur du ième élément du vecteur w.
  42  */
  43 __device__
  44 double w(long i)
  45     {
  46     double x = abs(cos(double(i)));
  47     for (long j = 1; j <= M_W; j++)
  48         x = x - (cos(x) - x) / (-sin(x) - 1);
  49
  50     /* Debug afin d'ajuster WI
  51     if (Indice1D::tid() == 0)
  52         printf("x: %.30f, WI: %.30f, x / WI: %.30f\n", x, WI, x / WI); */
  53
  54     return (x / WI) * sqrt(double(i));
  55     }
  56
  57 /*
  58  * 1) Chaque thread calcule un résultat intermediaire qu'il va ensuite placer en shared memory.
  59  * n: La taille des deux vecteurs.
  60  */
  61 __device__
  62 void reductionIntraThread(int n, double* tabSM)
  63     {
  64     const int NB_THREAD = Indice1D::nbThread();
  65     const int TID = Indice1D::tid();
  66     const int TID_LOCAL = Indice1D::tidLocal();
  67
  68     double threadResult = 0.0;
  69     int s = TID;
  70     while (s < n)
  71         {
  72         threadResult += v(s) * w(s);
  73         s += NB_THREAD;
  74         }
  75
  76     tabSM[TID_LOCAL] = threadResult;
  77     }
  78
  79 /*
  80  * Combine les résultats de 'tabSM' dans 'tabSM[0]'
  81  */
  82 __device__
  83 void combine(double* tabSM, int middle)
  84     {
  85     const int TID_LOCAL = Indice1D::tidLocal();
  86     const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
  87
  88     int s = TID_LOCAL;
  89     while (s < middle)
  90         {
  91         tabSM[s] = tabSM[s] + tabSM[s + middle];
  92         s += NB_THREAD_LOCAL;
  93         }
  94     }
  95
  96 /*
  97  * 2) La shared memory est réduite, le résultat est placé dans 'tabSM[0]'.
  98  */
  99 __device__
 100 void reductionIntraBlock(double* tabSM)
 101     {
 102     const int TAB_SIZE = blockDim.x;
 103     int middle = TAB_SIZE / 2;
 104
 105     while (middle > 0)
 106         {
 107         combine(tabSM, middle);
 108         middle /= 2;
 109         __syncthreads();
 110         }
 111     }
 112
 113 __device__
 114 int mutexReductionInterBlock = 0;
 115
 116 /*
 117  * 3) Le 'tabSM[0]' de chaque bloc est reduit dans 'ptrResult' qui se trouve en global memory.
 118  */
 119 __device__
 120 void reductionInterBlock(double* tabSM, double* ptrResult)
 121     {
 122     const int TID_LOCAL = Indice1D::tidLocal();
 123     if (TID_LOCAL == 0)
 124         {
 125         Lock lock(&mutexReductionInterBlock);
 126         lock.lock();
 127         (*ptrResult) += tabSM[0];
 128         lock.unlock();
 129
 130         // Si on travail en float (pas besoin de mutex) :
 131         // atomicAdd(ptrResult, float(tabSM[0]));
 132         }
 133     }
 134
 135 /**
 136  * La taille de la shared memory (en terme de # de sizeof(double)) doit
 137  * être égal à la taille des blocs.
 138  * n: La taille des deux vecteurs.
 139  * ptrResult: Le resultat du produit scalaire.
 140  */
 141 __global__
 142 void produitScalaire(int n, double* ptrResult)
 143     {
 144     extern __shared__ double tabSM[]; // Dynamic shared memory.;
 145
 146     // 1) Réduction intra-thread.
 147     reductionIntraThread(n, tabSM);
 148
 149     __syncthreads();
 150
 151     // 2) Réduction intra-block.
 152     reductionIntraBlock(tabSM);
 153
 154     // 3) Réduction inter-block.
 155     reductionInterBlock(tabSM, ptrResult);
 156     }
 157
 158 double resultatTheorique(long n)
 159     {
 160     n -= 1;
 161     return (n / 2.0) * (n + 1);
 162     }
 163
 164 bool produitScalaire()
 165     {
 166     const int N = 10000000; // Taille des deux vecteurs : 10 * 10^6.
 167
 168     // Allocation coté GPU en global memory (GM).
 169     double* ptrDevResult = 0;
 170     HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(double)));
 171     HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(double)));
 172
 173     // Paramètre de l'appel de la fonction sur le device.
 174     const dim3 dg(256, 1, 1);
 175     const dim3 db(256, 1, 1);
 176     Device::assertDim(dg, db);
 177     const size_t SMSize = db.x * sizeof(double); // 256 doubles;
 178
 179     produitScalaire<<<dg, db, SMSize>>>(N, ptrDevResult);
 180
 181     cudaDeviceSynchronize(); // Utilisé pour flusher les prints sur le stdout à partir du device.
 182
 183     double res;
 184     // Barrière implicite de synchronisation ('cudaMemCpy').
 185     HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(double), cudaMemcpyDeviceToHost));
 186
 187     const double resTheo = resultatTheorique(N);
 188
 189     cout.precision(30);
 190     cout << "Résultat : " << res << endl;
 191     cout << "Résultat théorique : " << resTheo << endl;
 192     cout << "Différence absolue : " << resTheo - res << endl;
 193     cout << "Différence relatif : " << 100 * (resTheo - res) / (resTheo + res) << " %" << endl;
 194
 195     return true;
 196     }