-
Notifications
You must be signed in to change notification settings - Fork 0
/
mcKerneLxCorrelated.cu
377 lines (281 loc) · 10.1 KB
/
mcKerneLxCorrelated.cu
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
/*
Progetto GPGPU
Luca Steccanella
W83_000009
Monte Carlo Hit or Miss
Generazione numeri casuali
*/
//Includes
#include "WarpCorrelated.cuh" //Libreria MWC64X, ci sono tre implementazioni diverse del warp
#include <unistd.h>
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <climits>
#include <vector>
#include <sstream>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <fstream>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <numeric>
//cuda incs
#include <cuda_runtime_api.h>
#include "device_launch_parameters.h"
#include <curand.h>
#include <curand_kernel.h>
#define _CRT_RAND_S
#define restrict __restrict__
using namespace std;
#define BLOCK_SIZE 1024
extern __shared__ unsigned rngShmem[];
__global__ void init(unsigned int seed, curandState_t *states)
{
int i = threadIdx.x + blockIdx.x * blockDim.x;
curand_init(seed, i, 0, states + i); //funzione proprietaria della curand: genera i semi
} //Kernel di inizializzazione: qui vengono generati i semi per curand
__global__ void mc_curand(int *hits, curandState_t *states, int c)
{
int i = threadIdx.x + blockIdx.x * blockDim.x;
int tHits = 0;
for (int j = 0; j < c; j++)
{
//tramite la funzione curand_uniform ottengo un valore casuale
float x = curand_uniform(states + i);
float y = curand_uniform(states + i);
//effettuo il passo dell'algoritmo mc_curand dove si verifica la somma dei quadrati
float z = (x * x + y * y);
if (z <= 1.0)
tHits++;
}
hits[i] = tHits;
} //Kernel che esegue l'algoritmo Monte Carlo Hit or Miss usando la curand
__global__ void mc_warp(int *hits, unsigned *state, int c)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
//inizializzo l'RNG
unsigned rngRegs[WarpCorrelated_REG_COUNT];
WarpCorrelated_LoadState(state, rngRegs, rngShmem);
//Ottengo i valori casuali ed effettuo i passi dell'algoritmo come suggerito nella documentazione
int tHits = 0;
for (int j = 0; j < c; j++)
{
unsigned long x = WarpCorrelated_Generate(rngRegs, rngShmem);
unsigned long y = WarpCorrelated_Generate(rngRegs, rngShmem);
x = (x * x) >> 3;
y = (y * y) >> 3;
if (x + y <= (1UL << 61))
{
tHits++;
}
}
hits[i] = tHits;
WarpCorrelated_SaveState(rngRegs, rngShmem, state);
} //Kernel che esegua l'algoritmo Monte Carlo Hit or Miss usando MWC64X
void check_error(cudaError_t err, const char *msg)
{
if (err != cudaSuccess)
{
fprintf(stderr, "%s : errore %d (%s)\n",
msg, err, cudaGetErrorString(err));
exit(err);
}
} //controllo errori
__global__ void red(int *vrid, const int *v1, int numels)
{
__shared__ int local[BLOCK_SIZE];
int i = threadIdx.x + blockIdx.x * blockDim.x;
int val = 0;
while (i < numels)
{
val += v1[i];
i += blockDim.x * gridDim.x;
}
local[threadIdx.x] = val;
for (int numthreads = BLOCK_SIZE / 2; numthreads > 0; numthreads /= 2)
{
__syncthreads();
if (threadIdx.x < numthreads)
{
local[threadIdx.x] += local[numthreads + threadIdx.x];
}
}
if (threadIdx.x == 0)
vrid[blockIdx.x] = local[0];
} //kernel che effettua la riduzione, come visto a lezione
cudaEvent_t evt[2]; //array di eventi ad uso generale
void do_rid(int *vrid, const int *v1, int N, int numBlocks)
{
cudaEventRecord(evt[0]);
red<<<numBlocks, BLOCK_SIZE>>>(vrid, v1, N);
cudaEventRecord(evt[1]);
cudaError_t error = cudaEventSynchronize(evt[1]);
check_error(error, "rid sync");
} //kernel che gesisce la riduzione
int main(int argc, char **argv)
{
//controlli avvio
if (argc != 3)
{
fprintf(stderr, "Inserire numero di campioni e il metodo di randomizazione; mcCUDA [campioni] [metodo]\n0: curand\n1: MWC64X (Warp Standard)\n");
return 1;
}
int N = atoi(argv[1]);
int t = atoi(argv[2]);
if (t < 0 || t > 1)
{
fprintf(stderr, "Metodo non valido\n");
return 1;
}
//allocazione di eventi, array hits, etc...
cudaError_t errore;
int devId = -1;
cudaDeviceProp devProps;
cudaGetDevice(&devId);
cudaGetDeviceProperties(&devProps, devId);
unsigned gridSize = devProps.multiProcessorCount;
unsigned totalThreads = BLOCK_SIZE * gridSize;
int totThrds = (int)totalThreads;
if (N < totThrds)
fprintf(stderr, "Usa almeno: %d campioni\n", totThrds);
int *d_hits;
size_t numbytes = totThrds * sizeof(int);
errore = cudaMalloc((void **)&d_hits, numbytes);
check_error(errore, "alloc hits");
errore = cudaMemset(d_hits, 0, numbytes);
check_error(errore, "memset hits");
cudaEvent_t prima_init, prima_mc;
cudaEvent_t dopo_init, dopo_mc;
errore = cudaEventCreate(&prima_init);
check_error(errore, "create prima init");
errore = cudaEventCreate(&prima_mc);
check_error(errore, "create prima mc");
errore = cudaEventCreate(&dopo_init);
check_error(errore, "create dopo init");
errore = cudaEventCreate(&dopo_mc);
check_error(errore, "create dopo mc");
float runtime;
size_t iobytes;
int numBlocks = devProps.multiProcessorCount * 6;
int Nt = N / totThrds;
printf("Nt: %d Threads: %d\n", Nt, totThrds);
if (t == 1)
{
//printf("Prima init\n");
/*
Di seguito inizializzo le variabili necessarie al funzionamento della libreria
e genero i semi che verranno usati dall'RNG; MWC64X non genera i semi nella GPU come CuRAND
*/
void *seedDevice = 0;
unsigned totalRngs = totalThreads / WarpCorrelated_K;
unsigned rngsPerBlock = BLOCK_SIZE / WarpCorrelated_K;
unsigned sharedMemBytesPerBlock = rngsPerBlock * WarpCorrelated_K * 4;
unsigned seedBytes = totalRngs * 4 * WarpCorrelated_STATE_WORDS;
std::vector<uint32_t> seedHost(seedBytes / 4);
cudaEventRecord(prima_init);
if (cudaMalloc(&seedDevice, seedBytes))
{
fprintf(stderr, "Error couldn't allocate state array of size %u\n", seedBytes);
exit(1);
}
int fr = open("/dev/urandom", O_RDONLY);
if (seedBytes != read(fr, &seedHost[0], seedBytes))
{
fprintf(stderr, "Couldn't seed RNGs.\n");
exit(1);
}
//cudaMemcpy(seedDevice, &seedHost[0], seedBytes, cudaMemcpyHostToDevice);
cudaEventRecord(dopo_init);
//printf("Dopo init / prima mc\n");
//eseguo il kernel
cudaEventRecord(prima_mc);
mc_warp<<<gridSize, BLOCK_SIZE, sharedMemBytesPerBlock>>>(d_hits, (unsigned *)seedDevice, Nt);
//printf("sync2\n");
errore = cudaDeviceSynchronize();
check_error(errore, "sync2");
cudaEventRecord(dopo_mc);
} //MWC64X
if (t == 0)
{
curandState_t *states;
errore = cudaMalloc((void **)&states, totThrds * sizeof(curandState_t));
check_error(errore, "alloc states");
//printf("Prima init\n");
cudaEventRecord(prima_init);
init<<<gridSize, BLOCK_SIZE>>>(time(0), states); //genero i semi nella GPU
cudaEventRecord(dopo_init);
//printf("sync1\n");
errore = cudaDeviceSynchronize();
check_error(errore, "sync1");
//printf("Dopo init / prima mc\n");
cudaEventRecord(prima_mc);
mc_curand<<<gridSize, BLOCK_SIZE>>>(d_hits, states, Nt); //eseguo il kernel
//printf("sync2\n");
errore = cudaDeviceSynchronize();
check_error(errore, "sync2");
cudaEventRecord(dopo_mc);
} //curand
//printf("Dopo mc / prima cpy\n");
float totRuntime = 0.0;
errore = cudaEventSynchronize(dopo_mc);
check_error(errore, "fine di tutto");
errore = cudaEventElapsedTime(&runtime, prima_init, dopo_init);
check_error(errore, "time init");
printf("init %gms\n", runtime);
totRuntime += runtime;
errore = cudaEventElapsedTime(&runtime, prima_mc, dopo_mc);
check_error(errore, "time mc");
iobytes = totThrds * Nt * sizeof(int);
printf("mc %gms, %g GB/s\n", runtime, iobytes / (runtime * 1.0e6));
totRuntime += runtime;
//effettuo la riduzione (chiamando l'apposito kernel) dell'array delle hits in modo di ottenerne il totale
int *dvrid = NULL;
int h_hits_sum = 0;
size_t ridsize = (numBlocks + 1) * sizeof(int);
errore = cudaEventCreate(evt);
check_error(errore, "evt create 0");
errore = cudaEventCreate(evt + 1);
check_error(errore, "evt create 1");
errore = cudaMalloc(&dvrid, ridsize);
check_error(errore, "malloc dvsum");
errore = cudaMemset(dvrid, -1, ridsize);
check_error(errore, "memset dvsum");
cudaEventRecord(evt[0]);
do_rid(dvrid + 1, d_hits, totThrds, numBlocks);
do_rid(dvrid, dvrid + 1, numBlocks, 1);
cudaEventRecord(evt[1]);
cudaEventSynchronize(evt[1]);
errore = cudaEventSynchronize(evt[1]);
check_error(errore, "rid");
errore = cudaEventElapsedTime(&runtime, evt[0], evt[1]);
check_error(errore, "time rid");
iobytes = (totThrds + 1) * sizeof(int);
printf("rid %gms, %g GB/s\n", runtime, iobytes / (runtime * 1.0e6));
totRuntime += runtime;
cudaEventRecord(evt[0]);
cudaMemcpy(&h_hits_sum, dvrid, sizeof(h_hits_sum), cudaMemcpyDeviceToHost);
cudaEventRecord(evt[1]);
errore = cudaEventSynchronize(evt[1]);
check_error(errore, "memcpy D2H");
cudaEventDestroy(evt[0]);
cudaEventDestroy(evt[1]);
cudaFree(dvrid);
printf("Total: %gms\n", totRuntime);
//effettuo il passo finale dell'algoritmo e stampo i risultati
/*std::vector<uint32_t>hitsHost(totalThreads, 0);
cudaMemcpy(&hitsHost[0], d_hits, 4*totalThreads, cudaMemcpyDeviceToHost);
int totalHits = std::accumulate(hitsHost.begin(), hitsHost.end(), 0.0);*/
cout << "Hits: " << h_hits_sum << "\n";
//cout << "HitsT: " << totalHits << "\n";
float AS = float(h_hits_sum) / (totThrds * float(Nt));
float estPi = AS * 4;
cout << "Est PI: " << estPi << "\n";
float rError = fabs(1.0 - estPi / M_PI);
cout << "Relative error: " << rError << "\n";
//cudaFree(states);
cudaFree(d_hits);
return 0;
}