首页 > 代码库 > Ubuntu12.04 之 CUDA 编程 (二) ～～～ GPU 程序加速
Ubuntu12.04 之 CUDA 编程 (二) ～～～ GPU 程序加速
关于 Ubuntu12.04 下 CUDA5.5 的安装请参看如下链接
Ubuntu-12.04 安装 CUDA-5.5
关于 Ubuntu12.04 下 CUDA5.5 程序的运行请参看如下链接
Ubuntu12.04 之 CUDA 编程 (一) ～～～ GPU 运行程序
1、程序的并行化
前一篇文章讲到了如何利用 CUDA5.5 在 GPU 中运行一个程序。通过程序的运行，我们看到了 GPU 确实可以作为一个运算器，但是，我们在前面的例子中并没有正真的发挥 GPU 并行处理程序的能力，也就是说之前的例子只利用了 GPU 的一个线程，没有发挥程序的并行性。
先来说说 CUDA5.5 中 GPU 的架构。它是由 grid 组成， 每个 grid 又可以由 block 组成， 而每个 block 又可以细分为 thread。所以，线程是我们处理的最小的单元了。
接下来的例子通过修改前一个例子，把数组分割成若干个组（每个组由一个线程实现），每个组计算出一个和，然后在 CPU 中将分组的这几个和加在一起，得到最终的结果。这种思想叫做 归约 。其实和分治思想差不多，就是先将大规模问题分解为小规模的问题，最后这些小规模问题整合得到最终解。
由于我的 GPU 支持的块内最大的线程数是 512 个，即 cudaGetDeviceProperties 中的 maxThreadsPerBlock 属性。如何获取这个属性，请参看如下链接：Runtime API 函数解析。 我们使用 512 个线程来实现并行加速。
好了，接下来就是写程序的时候了。
首先，在程序头部增加一个关于 线程数量 的宏定义：
// ======== define area ======== #define DATA_SIZE 1048576 // 1M #define THREAD_NUM 512 // thread num
接着，修改内核函数：
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
const int size = DATA_SIZE / THREAD_NUM;
const int tid = threadIdx.x;
int tmp_sum = 0;
clock_t start;
if (tid == 0) start = clock();
for (int i = tid * size; i < (tid + 1) * size; i++) {
tmp_sum += data[i] * data[i];
}
sum[tid] = tmp_sum;
if (tid == 0) *time = clock() - start;
}此内核程序的目的是把输入的数据分摊到 512 个线程上去计算 部分和，并且 512 个部分和存放到 sum 数组中，最后在 CPU 中对 512 个部分和求和得到最终结果。
此处对数据的遍历方式请注意一下，我们是根据顺序给每一个线程的，也就是如下表格所示：
| 线程编号 | 数据下标 |
| 0 | 0 ～ 2047 |
| ... ... | ... ... |
| 511 | 1046528 ～ 1048575 |
其次，修改主函数部分，也就是修改 部分和 部分。
// init CUDA device
if (!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
// generate rand datas
generateData(data, DATA_SIZE);
// malloc space for datas in GPU
int *gpuData, *sum;
clock_t *time;
cudaMalloc((void**) &gpuData, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &sum, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpuData, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
// calculate the squares‘s sum
squaresSum<<<1, THREAD_NUM, 0>>>(gpuData, sum, time);
// copy the result from GPU to HOST
int result[THREAD_NUM];
clock_t time_used;
cudaMemcpy(&result, sum, sizeof(int) * THREAD_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost);
// free GPU spaces
cudaFree(gpuData);
cudaFree(sum);
cudaFree(time);此处只是把 sum 改成数组，其他都没有变化。接下来就是添加 CPU 求取 部分和 的代码了：
// print result
int tmp_result = 0;
for (int i = 0; i < THREAD_NUM; ++i) {
tmp_result += result[i];
}
printf("(GPU) sum:%d time:%ld\n", tmp_result, time_used);
// CPU calculate
tmp_result = 0;
for (int i = 0; i < DATA_SIZE; ++i){
tmp_result += data[i] * data[i];
}
printf("(CPU) sum:%d\n", tmp_result);至此，程序已经修改完成，编译运行：
经过修改以后的程序，比之前的快了将近 36 倍，可见并行化处理还是存在优势的。不过仔细想一下，我们使用了 512 个线程， 可是性能怎么才提升了 36 倍，不应该是 512 倍吗？？？
这里就涉及到内存的存取模式了，显卡上面的内存是 DRAM，是效率最高的存取方式，它是一种连续的存取方式。前面我们的程序确实的连续读取的呀，都挨个读取了，怎么还是没有达到预期的效果呢？？？
这里还需要考虑 thread 的执行方式，前一节说到，当一个 thread 在等待内存数据的时候， GPU 就会切换到下一个 thread。所以，实际执行的顺序类似于 thread0 --> thread1 --> ... ... --> thread511
这就导致了同一个 thread 在读取内存是连续的， 但是对于整体而言，执行的过程中读取就不是连续的了（这里自己仔细想想，就明白了）。所以，正确的做法如下表格所示：
| 线程编号 | 数据下标 |
| 0 | 0,512,... ..., |
| ... ... | ... ... |
| 511 | 511,1023,... ... |
根据这个原理，修改内核函数如下：
for (int i = tid; i < DATA_SIZE; i += THREAD_NUM) {
tmp_sum += data[i] * data[i];
}编译运行如下：
修改后程序，比之前的又快了 13 倍左右，可见，对内存的读取方式对于性能的影响很大。
至此，并行化后的程序较未并行化之前的程序，速度上快了 493 倍左右，可见，基本上发挥了 512 个线程的优势。
让我们再来分析一下性能：
此 GPU 消耗的时钟周期： 1595788 cycles
GeForce G 103M 的 clockRate： 1.6 GHz
所以可以计算出 GPU 上运行时间是： 时钟周期 / clockRate = 997.3675 us
1 M 个 int 型数据有 4M Byte 的数据量，实际使用的 GPU 内存带宽是：数据量 / 运行时间 = 4.01 GB/s
再来看看我的 GPU GeForce 103m 的内存带宽：运行 SDK 目录下面 /samples/1_Utilities/bandwidthTest
[CUDA Bandwidth Test] - Starting... Running on... Device 0: GeForce G 103M Quick Mode Host to Device Bandwidth, 1 Device(s) PINNED Memory Transfers Transfer Size (Bytes) Bandwidth(MB/s) 33554432 2513.8 Device to Host Bandwidth, 1 Device(s) PINNED Memory Transfers Transfer Size (Bytes) Bandwidth(MB/s) 33554432 1631.5 Device to Device Bandwidth, 1 Device(s) PINNED Memory Transfers Transfer Size (Bytes) Bandwidth(MB/s) 33554432 4580.1 Result = PASS
通过与系统参数的对比，可以知道，基本上达到了系统的极限性能。
2、并行化——进阶
前面提到了 grid、block、thread 三者之间的关系，知道了他们之间是逐渐包含的关系。我们在上面的程序中通过使用 512 个线程达到了 493 倍左右的性能提升，那么是不是可以继续得到提升呢？？？
答案是肯定的，这就要进一步考虑 GPU 的并行化处理了。前面的程序只是使用了单个 block 下的 512 个线程，那么，我们可不可以使用多个 block 来实现？？？
对，就是利用这个思想，达到进一步的并行化。这里使用 8 个 block * 64 threads = 512 threads 实现。
首先，修改主函数宏定义，定义 块数量：
// ======== define area ======== #define DATA_SIZE 1048576 // 1M #define BLOCK_NUM 8 // block num #define THREAD_NUM 64 // thread num
接下来，修改内核函数：
_global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int tmp_sum = 0;
if (tid == 0) time[bid] = clock();
for (int i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
tmp_sum += data[i] * data[i];
}
sum[bid * THREAD_NUM + tid] = tmp_sum;
if (tid == 0) time[bid + BLOCK_NUM] = clock();
}注意：这里的内存遍历方式和前面讲的是一致的，理解一下。
同时记录的时间是一个块的开始和结束时间，因为这里我们最后需要计算的是最早开始和最晚结束的两个时间差，即求出最糟糕的时间。
然后，就是主函数里面的具体实现了：
// init CUDA device
if (!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
// generate rand datas
generateData(data, DATA_SIZE);
// malloc space for datas in GPU
int *gpuData, *sum;
clock_t *time;
cudaMalloc((void**) &gpuData, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &sum, sizeof(int) * THREAD_NUM * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpuData, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
// calculate the squares‘s sum
squaresSum<<<BLOCK_NUM, THREAD_NUM, 0>>>(gpuData, sum, time);
// copy the result from GPU to HOST
int result[THREAD_NUM * BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&result, sum, sizeof(int) * THREAD_NUM * BLOCK_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost);
// free GPU spaces
cudaFree(gpuData);
cudaFree(sum);
cudaFree(time);存取的是开始和结束两个时间，所以时间数组要加倍。
然后就是 CPU 计算求和结果以及 消耗 时钟 cycles：
// print result
int tmp_result = 0;
for (int i = 0; i < THREAD_NUM * BLOCK_NUM; ++i) {
tmp_result += result[i];
}
clock_t min_start, max_end;
min_start = time_used[0];
max_end = time_used[BLOCK_NUM];
for (int i = 1; i < BLOCK_NUM; ++i) {
if (min_start > time_used[i]) min_start = time_used[i];
if (max_end < time_used[i + BLOCK_NUM]) max_end = time_used[i + BLOCK_NUM];
}
printf("(GPU) sum:%d time:%ld\n", tmp_result, max_end - min_start);
// CPU calculate
tmp_result = 0;
for (int i = 0; i < DATA_SIZE; ++i) {
tmp_result += data[i] * data[i];
}
printf("(CPU) sum:%d\n", tmp_result);编译运行程序如下所示：
jping@CJP-PC:~/cuda/2.0_SquareSum$ ./squareSum CUDA initialized. (GPU) sum:29909398 time:1601876 (CPU) sum:29909398 jping@CJP-PC:~/cuda/2.0_SquareSum$
性能与直接使用 512 个线程基本一致。
3、线程同步和共享内存
前面的程序，计算求和的工作在 CPU 中完成，总共需要在 CPU 中做 512 次加法运算，那么有没有办法减少 CPU 中执行加法的次数呢？？？
可以通过同步和共享内存技术，实现在 GPU 上的 block 块内求取部分和，这样最后只需要在 CPU 计算 16 个和就可以了。具体实现方法如下：
首先，在修改内核函数，定义一块共享内存，用 __shared__ 指示：
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
__shared__ int shared[BLOCK_NUM];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
if (tid == 0) time[bid] = clock();
shared[tid] = 0;
for (int i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += data[i] * data[i];
}
__syncthreads();
if (tid == 0) {
for (int i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
sum[bid] = shared[0];
}
if (tid == 0) time[bid + BLOCK_NUM] = clock();
}利用 __shared__ 声明的变量是 shared memory，每个 block 中，各个 thread 之间对于共享内存是共享的，利用的是 GPU 上的内存，所以速度很快，不必担心 latency 的问题。
__syncthreads() 函数是 CUDA 的内部函数，表示所有 threads 都必须同步到这个点，才会执行接下来的代码。我们要做的就是等待每个 thread 计算结束以后，再来计算部分和，所以同步是必不可少的环节。把每个 block 的部分和计算到 shared[0] 里面。
接着，修改主函数部分：
int main(int argc, char const *argv[])
{
// init CUDA device
if (!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
// generate rand datas
generateData(data, DATA_SIZE);
// malloc space for datas in GPU
int *gpuData, *sum;
clock_t *time;
cudaMalloc((void**) &gpuData, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &sum, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpuData, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
// calculate the squares‘s sum
squaresSum<<<BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int)>>>(gpuData, sum, time);
// copy the result from GPU to HOST
int result[BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&result, sum, sizeof(int) * BLOCK_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost);
// free GPU spaces
cudaFree(gpuData);
cudaFree(sum);
cudaFree(time);
// print result
int tmp_result = 0;
for (int i = 0; i < BLOCK_NUM; ++i) {
tmp_result += result[i];
}
clock_t min_start, max_end;
min_start = time_used[0];
max_end = time_used[BLOCK_NUM];
for (int i = 1; i < BLOCK_NUM; ++i) {
if (min_start > time_used[i]) min_start = time_used[i];
if (max_end < time_used[i + BLOCK_NUM]) max_end = time_used[i + BLOCK_NUM];
}
printf("(GPU) sum:%d time:%ld\n", tmp_result, max_end - min_start);
// CPU calculate
tmp_result = 0;
for (int i = 0; i < DATA_SIZE; ++i) {
tmp_result += data[i] * data[i];
}
printf("(CPU) sum:%d\n", tmp_result);
return 0;
}相对与前一个版本，只需要修改 sum 数组的大小，以及调用核函数的代码。
编译运行代码，结果如下：
jping@CJP-PC:~/cuda/2.0_SquareSum$ ./squareSum CUDA initialized. (GPU) sum:29909398 time:1812756 (CPU) sum:29909398 jping@CJP-PC:~/cuda/2.0_SquareSum$
其实和前一版程序相比，时间上没有什么优势，原因在于，我们需要在 GPU 中额外运行求和的这部分代码，导致了运行周期的变长，不过相应的，在 CPU 中的运行时间会减少。
我们在这个程序中，只当每个 block 的 thread0 的时候，计算求和的工作，这样做影响了执行的效率，其实求和可以并行化处理的，也就是通过加法树来实现并行化。举个例子，要计算 8 个数的和，我们没必要用一个 for 循环，逐个相加，而是可以通过第一级流水线实现两两相加，变成 4 个数，第二级流水实现两两相加，变成 2 个数，第三级流水实现两两相加，求得最后的和。
下面通过加法树的方法，实现最后的求和，修改内核函数如下：
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
__shared__ int shared[BLOCK_NUM];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int offset = THREAD_NUM / 2;
if (tid == 0) time[bid] = clock();
shared[tid] = 0;
for (int i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += data[i] * data[i];
}
__syncthreads();
while (offset > 0) {
if (tid < offset) {
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
if (tid == 0) {
sum[bid] = shared[0];
time[bid + BLOCK_NUM] = clock();
}
}此程序实现的就是上诉描述的加法树的结构，注意这里第二个 __syncthreads() 的使用，也就是说，要进行下一级流水线的计算，必须建立在前一级必须已经计算完毕的情况下。
主函数部分不许要修改，最后编译运行结果如下：
jping@CJP-PC:~/cuda/2.0_SquareSum$ ./squareSum CUDA initialized. (GPU) sum:29909398 time:1799680 (CPU) sum:29909398 jping@CJP-PC:~/cuda/2.0_SquareSum$
性能有一部分的改善。
通过使用 GPU 的并行化编程，确实对性能会有很大程度上的提升。由于受限于 Geforce 103m 的内存带宽，程序只能优化到这一步，关于是否还有其他的方式优化，有待进一步学习。
最后给出最终的代码：
/* *******************************************************************
##### File Name: squareSum.cu
##### File Func: calculate the sum of inputs‘s square
##### Author: Caijinping
##### E-mail: caijinping220@gmail.com
##### Create Time: 2014-5-7
* ********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
// ======== define area ========
#define DATA_SIZE 1048576 // 1M
#define BLOCK_NUM 8 // block num
#define THREAD_NUM 64 // thread num
// ======== global area ========
int data[DATA_SIZE];
void printDeviceProp(const cudaDeviceProp &prop);
bool InitCUDA();
void generateData(int *data, int size);
__global__ static void squaresSum(int *data, int *sum, clock_t *time);
int main(int argc, char const *argv[])
{
// init CUDA device
if (!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
// generate rand datas
generateData(data, DATA_SIZE);
// malloc space for datas in GPU
int *gpuData, *sum;
clock_t *time;
cudaMalloc((void**) &gpuData, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &sum, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpuData, data, sizeof(int) * DATA_SIZE, cudaMemcpyHostToDevice);
// calculate the squares‘s sum
squaresSum<<<BLOCK_NUM, THREAD_NUM, THREAD_NUM * sizeof(int)>>>(gpuData, sum, time);
// copy the result from GPU to HOST
int result[BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&result, sum, sizeof(int) * BLOCK_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost);
// free GPU spaces
cudaFree(gpuData);
cudaFree(sum);
cudaFree(time);
// print result
int tmp_result = 0;
for (int i = 0; i < BLOCK_NUM; ++i) {
tmp_result += result[i];
}
clock_t min_start, max_end;
min_start = time_used[0];
max_end = time_used[BLOCK_NUM];
for (int i = 1; i < BLOCK_NUM; ++i) {
if (min_start > time_used[i]) min_start = time_used[i];
if (max_end < time_used[i + BLOCK_NUM]) max_end = time_used[i + BLOCK_NUM];
}
printf("(GPU) sum:%d time:%ld\n", tmp_result, max_end - min_start);
// CPU calculate
tmp_result = 0;
for (int i = 0; i < DATA_SIZE; ++i) {
tmp_result += data[i] * data[i];
}
printf("(CPU) sum:%d\n", tmp_result);
return 0;
}
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
__shared__ int shared[BLOCK_NUM];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int offset = THREAD_NUM / 2;
if (tid == 0) time[bid] = clock();
shared[tid] = 0;
for (int i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += data[i] * data[i];
}
__syncthreads();
while (offset > 0) {
if (tid < offset) {
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
if (tid == 0) {
sum[bid] = shared[0];
time[bid + BLOCK_NUM] = clock();
}
}
// ======== used to generate rand datas ========
void generateData(int *data, int size)
{
for (int i = 0; i < size; ++i) {
data[i] = rand() % 10;
}
}
void printDeviceProp(const cudaDeviceProp &prop)
{
printf("Device Name : %s.\n", prop.name);
printf("totalGlobalMem : %d.\n", prop.totalGlobalMem);
printf("sharedMemPerBlock : %d.\n", prop.sharedMemPerBlock);
printf("regsPerBlock : %d.\n", prop.regsPerBlock);
printf("warpSize : %d.\n", prop.warpSize);
printf("memPitch : %d.\n", prop.memPitch);
printf("maxThreadsPerBlock : %d.\n", prop.maxThreadsPerBlock);
printf("maxThreadsDim[0 - 2] : %d %d %d.\n", prop.maxThreadsDim[0], prop.maxThreadsDim[1], prop.maxThreadsDim[2]);
printf("maxGridSize[0 - 2] : %d %d %d.\n", prop.maxGridSize[0], prop.maxGridSize[1], prop.maxGridSize[2]);
printf("totalConstMem : %d.\n", prop.totalConstMem);
printf("major.minor : %d.%d.\n", prop.major, prop.minor);
printf("clockRate : %d.\n", prop.clockRate);
printf("textureAlignment : %d.\n", prop.textureAlignment);
printf("deviceOverlap : %d.\n", prop.deviceOverlap);
printf("multiProcessorCount : %d.\n", prop.multiProcessorCount);
}
bool InitCUDA()
{
//used to count the device numbers
int count;
// get the cuda device count
cudaGetDeviceCount(&count);
if (count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
// find the device >= 1.X
int i;
for (i = 0; i < count; ++i) {
cudaDeviceProp prop;
if (cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
if (prop.major >= 1) {
//printDeviceProp(prop);
break;
}
}
}
// if can‘t find the device
if (i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
// set cuda device
cudaSetDevice(i);
return true;
}
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。