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

   1 #include <iostream>
   2 #include <cmath>
   3 using namespace std;
   4
   5 #include "Indice1D.h"
   6 #include "cudaTools.h"
   7 #include "Device.h"
   8
   9 #define M_V 200
  10 #define M_W 200
  11
  12 #define VI 1.4422495703074083
  13 #define WI 0.7390850782394409
  14
  15 /*
  16  * Renvoie la valeur du ième élément du vecteur v.
  17  */
  18 __device__
  19 double v(long i)
  20     {
  21     double x = 1.5 + abs(cos(double(i)));
  22     for (long j = 1; j <= M_V; j++)
  23         {
  24         const double xCarre = x * x;
  25         x = x - (xCarre * x - 3) / (3 * xCarre);
  26         }
  27     return (x / VI) * sqrt(double(i));
  28     }
  29
  30 /*
  31  * Renvoie la valeur du ième élément du vecteur w.
  32  */
  33 __device__
  34 double w(long i)
  35     {
  36     double x = abs(cos(double(i)));
  37     for (long j = 1; j <= M_W; j++)
  38         x = x - (cos(x) - x) / (-sin(x) - 1);
  39     return (x / WI) * sqrt(double(i));
  40     }
  41
  42 /*
  43  * n: La taille des deux vecteurs.
  44  */
  45 __device__ void reductionIntraThread(int n, double* tabResultSM)
  46     {
  47     const int NB_THREAD = Indice1D::nbThread();
  48     const int TID = Indice1D::tid();
  49     const int TID_LOCAL = Indice1D::tidLocal();
  50
  51     double threadResult = 0.0;
  52     int s = TID;
  53     while (s < n)
  54         {
  55         threadResult += v(s) * w(s);
  56         s += NB_THREAD;
  57         }
  58     tabResultSM[TID_LOCAL] = threadResult;
  59     }
  60
  61 /*
  62  * Combine les résultats de 'tabResultSM' dans 'tabResulSM[0]'
  63  */
  64 __device__ void combine(double* tabResultSM, int middle)
  65     {
  66         const int TID_LOCAL = Indice1D::tidLocal();
  67         const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
  68
  69         int s = TID_LOCAL;
  70         while (s < middle)
  71             {
  72             tabResultSM[s] = tabResultSM[s] + tabResultSM[s + middle];
  73             s += NB_THREAD_LOCAL;
  74             }
  75     }
  76
  77 __device__ void reductionIntraBlock(double* tabResultSM)
  78     {
  79     const int TAB_SIZE = blockDim.x;
  80     int middle = TAB_SIZE / 2;
  81
  82     while (middle > 0)
  83         {
  84         combine(tabResultSM, middle);
  85         middle /= 2;
  86         __syncthreads();
  87         }
  88     }
  89
  90 __device__ void reductionInterBlock(double* tabResultSM, float* ptrResult)
  91     {
  92     const int TID_LOCAL = Indice1D::tidLocal();
  93     if (TID_LOCAL == 0)
  94         {
  95         atomicAdd(ptrResult, float(tabResultSM[0]));
  96         }
  97     }
  98
  99 /**
 100  * La taille de la shared memory (en terme de # de sizeof(double)) doit
 101  * être égal à la taille des blocs.
 102  * n: La taille des deux vecteurs.
 103  * ptrResult: Le resultat du produit scalaire.
 104  */
 105 __global__
 106 void produitScalaire(int n, float* ptrResult)
 107     {
 108     extern __shared__ double tabResultSM[]; // Shared memory.
 109
 110     // 1) Réduction intra-thread.
 111     reductionIntraThread(n, tabResultSM);
 112
 113     __syncthreads();
 114
 115     // 2) Réduction intra-block.
 116     reductionIntraBlock(tabResultSM);
 117
 118     // 3) Réduction inter-block.
 119     reductionInterBlock(tabResultSM, ptrResult);
 120     }
 121
 122 double resultatTheorique(long n)
 123 {
 124     n -= 1;
 125     return (n / 2.0) * (n+1);
 126 }
 127
 128 bool produitScalaire()
 129     {
 130     const int N = 100000000; // Taille des deux vecteurs.
 131
 132     // Allocation coté GPU en global memory (GM).
 133     float* ptrDevResult = 0;
 134     HANDLE_ERROR(cudaMalloc(&ptrDevResult, sizeof(float)));
 135     HANDLE_ERROR(cudaMemset(ptrDevResult, 0, sizeof(float)));
 136
 137     // Paramètre de l'appel de la fonction sur le device.
 138     const dim3 dg(256, 1, 1);
 139     const dim3 db(256, 1, 1);
 140     Device::assertDim(dg, db);
 141     const size_t SMSize = db.x * sizeof(double); // 256 doubles;
 142
 143     produitScalaire<<<dg, db, SMSize>>>(N, ptrDevResult);
 144
 145     float res;
 146     // Barrière implicite de synchronisation ('cudaMemCpy').
 147     HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(float), cudaMemcpyDeviceToHost));
 148
 149     double resTheo = resultatTheorique(N);
 150
 151     cout.precision(10);
 152     cout << "Résultat : " << res << endl;
 153     cout << "Résultat théorique : " << resTheo << endl;
 154     cout << "Différence absolue : " << resTheo - res << endl;
 155     cout << "Différence relatif : " << 100 * (resTheo - res) / (resTheo + res) << " %" << endl;
 156
 157     return true;
 158     }