CUDA积分法求PI

本文将设计并实现一个并行计算程序,该程序利用积分法计算 PI 的值。

CUDA 入门书籍推荐

《GPU高性能编程CUDA实战》

基本思路

\[ \pi = 4\int_{0}^{1}\frac{1}{1+x^{2}}dx \]

CUDA 程序执行的主要步骤

  1. CPU在GPU上申请空间

    cudaMalloc(起始地址,大小)

  2. CPU将数据从内存拷贝到显存

    cudaMemcpy(源,目标,大小,方向)

  3. CPU启动GPU上的内核进行计算

    kernel_name <<<blocks,threads>>>(函数参数)

  4. CPU将处理结果从显存拷贝到内存

    cudaMemcpy(源,目标,大小,方向)

源码分析

PI.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
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

static void HandleError(cudaError_t err,
const char *file,
int line) {
if (err != cudaSuccess) {
printf("%s in %s at line %d\n", cudaGetErrorString(err),
file, line);
exit(EXIT_FAILURE);
}
}

#define HANDLE_ERROR(err) (HandleError(err, __FILE__, __LINE__))
#define HANDLE_NULL(a) {if (a == NULL) { \
printf("Host memory failed in %s at line %d\n", \
__FILE__, __LINE__); \
exit(EXIT_FAILURE);}}

const int N = 1024 * 1024 * 64; //积分时划分的份数
const int threadsPerBlock = 256; //block中的线程数
const int blocksPerGrid = 64; //grid中的block数

//缺省__host__,表明CPU运行,CPU调用
double function_for_cpu(double x) {
return 4 / (1 + x * x);
}

//__device__修饰,只能被内核函数调用,表明GPU运行,GPU调用
__device__ double function_for_gpu(double x) {
return 4 / (1 + x * x);
}

/**
* __global__修饰,内核函数,表明GPU运行,CPU调用
* 在GPU上使用积分法并行计算PI
* @param a 积分下界
* @param b 积分上界
* @param integral 存储积分结果
* @return
*/
__global__ void trap(double *a, double *b, double *integral) {

/**
* __shared__修饰,声明一个共享内存缓冲区,名字为cache
* 表示数据存放在共享存储器中,每个线程块都有该变量的一个副本;
* 只有在块内的线程可以访问,其它块内的线程不能访问;
*/
__shared__ double cache[threadsPerBlock];

/**
* 该步的目的是计算初始线程索引,其中threadIdx, blockIdx, blockDim都是内置变量
* threadIdx是存储线程信息的结构体,对于线程0来说,threadIdx.x=0;
* blockDim.x表示block在x维度的线程数量,本例中使用的是一维线程块,
* 因此只需用到blockDim.x
* blockIdx.x表示block的索引,对于第一个线程块来说,blockIdx.x=0;
* 对于第二个线程块来说,blockIdx.x=1...
* 在计算tid线程索引时,需要要在threadIdx.x 的基础上加上一个基地址,
* 实际上就是将二维索引空间转换为线性空间
*/
int tid = threadIdx.x + blockIdx.x * blockDim.x;

/**
* 共享内存缓存中的偏移就等于线程索引,线程块索引与该偏移无关,
* 因为每个线程块都拥有该共享内存的私有副本
*/
int cacheIndex = threadIdx.x;

double x, temp = 0;
while (tid < N) {
x = *a + (double)(*b - *a) / N * (tid + 0.5);
temp += function_for_gpu(x);

/**
* 在每个线程计算完当前索引上的任务后,需要对索引进行递增,
* 其中,递增的步长为线程格中正在运行的线程数量,
* 这个数值等于线程块中的线程数量乘以线程格中线程块的数量,
* 即 blockDim.x * gridDim.x
*
* 该方法类似于多CPU或多核CPU的并行,数据迭代的增量不是1,
* 而是CPU的数量;在GPU实现中,一般将并行线程数量看做处理器的数量
*/
tid += blockDim.x * gridDim.x;
}

//设置cache中相应位置上的值
cache[cacheIndex] = temp;

/**
* 对线程块中的线程进行同步,该操作用于确保
* 所有对共享数组cache[]写入操作在读取cache之前完成
*/
__syncthreads();

/**
* 对于归约运算来说,以下代码要求threadsPerBlock必须是2的幂,
* 因为每次合并,要求分成的两部分数组的长度要一致
* 基本思想:
* 每个线程将cache[]中的两个值相加起来,然后将结果保存回cache[]
* 由于每个线程都将两个值合并为一个值,那么在完成这个步骤后,
* 得到的结果数量就是计算开始时数值数量的一半。接着,对这一半
* 进行相同操作,直到cache[]中256个值归约为1个值
*/
int i = blockDim.x / 2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}

//将这个值保存到全局内存后,内核函数结束
//这里使用了索引为0的线程将cache[0]写入全局内存
if (cacheIndex == 0)
integral[blockIdx.x] = cache[0];
}

/**
* 在CPU上使用积分法串行计算PI
* @param a 积分下界
* @param b 积分上界
* @param integral 存储积分结果
*/
void trap_by_cpu(double a, double b, double *integral) {
int i;
double x, temp = 0;
for (i = 0; i < N; i++) {
x = a + (double)(b - a) / N * (i + 0.5);
temp += function_for_cpu(x);
}
temp *= (double)(b - a) / N;
*integral = temp;
}

int main(void) {
double integral;
double *partial_integral;
double a, b;

double *dev_partial_integral;
double *dev_a, *dev_b;
cudaEvent_t start,stop;
float tm;

clock_t clockBegin, clockEnd;
float duration;

a = 0;
b = 1;

//使用CPU计算PI的值
clockBegin = clock();
trap_by_cpu(a, b, &integral);
clockEnd = clock();
duration = (float)1000 * (clockEnd - clockBegin) / CLOCKS_PER_SEC;
printf("CPU Result: %.20lf\n", integral);
printf("CPU Elapsed time: %.6lfms\n", duration);

getchar();

//使用GPU+CPU计算PI的值
//使用event计算时间
cudaEventCreate(&start); //创建event
cudaEventCreate(&stop); //创建event
cudaEventRecord(start, 0); //记录当前时间

// 申请CPU存储空间
partial_integral = (double*)malloc(blocksPerGrid * sizeof(double));

// 申请GPU存储空间
HANDLE_ERROR(cudaMalloc((void**)&dev_a, sizeof(double)));
HANDLE_ERROR(cudaMalloc((void**)&dev_b, sizeof(double)));
HANDLE_ERROR(cudaMalloc((void**)&dev_partial_integral,
blocksPerGrid * sizeof(double)));

//将'a'和'b'复制到GPU
HANDLE_ERROR(cudaMemcpy(dev_a, &a, sizeof(double), cudaMemcpyHostToDevice));
HANDLE_ERROR(cudaMemcpy(dev_b, &b, sizeof(double), cudaMemcpyHostToDevice));

//启动内核函数,启动的线程块的数量是64个,每个线程块的线程数量是256
trap<<<blocksPerGrid, threadsPerBlock>>>(dev_a, dev_b, dev_partial_integral);

//将计算结果数组'dev_partial_integral'从GPU复制到CPU
HANDLE_ERROR(cudaMemcpy(partial_integral, dev_partial_integral,
blocksPerGrid * sizeof(double),
cudaMemcpyDeviceToHost));

//在CPU上进行归约操作,得到最终的计算结果
integral = 0;
for (int i = 0; i < blocksPerGrid; i++) {
integral += partial_integral[i];
}
integral *= (double)(b - a) / N;

cudaEventRecord(stop,0); //记录当前时间
cudaEventSynchronize(stop); //等待stop event完成
cudaEventElapsedTime(&tm, start, stop); //计算时间差(毫秒级)
printf("GPU Result: %.20lf\n", integral);
printf("GPU Elapsed time:%.6f ms.\n", tm);

//释放GPU内存
HANDLE_ERROR(cudaFree(dev_a));
HANDLE_ERROR(cudaFree(dev_b));
HANDLE_ERROR(cudaFree( dev_partial_integral));

//释放CPU内存
free(partial_integral);

getchar();
}

References

《GPU高性能编程CUDA实战》