1 #include "ProduitScalaire.h"
16 #define VI 1.4422495703074083
17 #define WI 0.739085133215160672293109200837
20 * Renvoie la valeur du ième élément du vecteur v.
25 double x = 1.5 + abs(cos(double(i)));
26 for (long j = 1; j <= M_V; j++)
28 const double xCarre = x * x;
29 x = x - (xCarre * x - 3) / (3 * xCarre);
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);
37 return (x / VI) * sqrt(double(i)); // x / VI doit être égal à 1.
41 * Renvoie la valeur du ième élément du vecteur w.
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);
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); */
54 return (x / WI) * sqrt(double(i));
58 * 1) Chaque thread calcule un résultat intermediaire qu'il va ensuite placer en shared memory.
59 * n: La taille des deux vecteurs.
62 void reductionIntraThread(int n, double* tabSM)
64 const int NB_THREAD = Indice1D::nbThread();
65 const int TID = Indice1D::tid();
66 const int TID_LOCAL = Indice1D::tidLocal();
68 double threadResult = 0.0;
72 threadResult += v(s) * w(s);
76 tabSM[TID_LOCAL] = threadResult;
80 * Combine les résultats de 'tabSM' dans 'tabSM[0]'
83 void combine(double* tabSM, int middle)
85 const int TID_LOCAL = Indice1D::tidLocal();
86 const int NB_THREAD_LOCAL = Indice1D::nbThreadBlock();
91 tabSM[s] = tabSM[s] + tabSM[s + middle];
97 * 2) La shared memory est réduite, le résultat est placé dans 'tabSM[0]'.
100 void reductionIntraBlock(double* tabSM)
102 const int TAB_SIZE = blockDim.x;
103 int middle = TAB_SIZE / 2;
107 combine(tabSM, middle);
114 int mutexReductionInterBlock = 0;
117 * 3) Le 'tabSM[0]' de chaque bloc est reduit dans 'ptrResult' qui se trouve en global memory.
120 void reductionInterBlock(double* tabSM, double* ptrResult)
122 const int TID_LOCAL = Indice1D::tidLocal();
125 Lock lock(&mutexReductionInterBlock);
127 (*ptrResult) += tabSM[0];
130 // Si on travail en float (pas besoin de mutex) :
131 // atomicAdd(ptrResult, float(tabSM[0]));
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.
142 void produitScalaire(int n, double* ptrResult)
144 extern __shared__ double tabSM[]; // Dynamic shared memory.;
146 // 1) Réduction intra-thread.
147 reductionIntraThread(n, tabSM);
151 // 2) Réduction intra-block.
152 reductionIntraBlock(tabSM);
154 // 3) Réduction inter-block.
155 reductionInterBlock(tabSM, ptrResult);
158 double resultatTheorique(long n)
161 return (n / 2.0) * (n + 1);
164 bool produitScalaire()
166 const int N = 10000000; // Taille des deux vecteurs : 10 * 10^6.
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)));
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;
179 produitScalaire<<<dg, db, SMSize>>>(N, ptrDevResult);
181 cudaDeviceSynchronize(); // Utilisé pour flusher les prints sur le stdout à partir du device.
184 // Barrière implicite de synchronisation ('cudaMemCpy').
185 HANDLE_ERROR(cudaMemcpy(&res, ptrDevResult, sizeof(double), cudaMemcpyDeviceToHost));
187 const double resTheo = resultatTheorique(N);
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;