首页 > 代码库 > 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 个部分和求和得到最终结果。


此处对数据的遍历方式请注意一下,我们是根据顺序给每一个线程的,也就是如下表格所示:

数据顺序遍历方式
线程编号数据下标
00 ~ 2047
... ...... ...
5111046528 ~ 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 在读取内存是连续的, 但是对于整体而言,执行的过程中读取就不是连续的了(这里自己仔细想想,就明白了)。所以,正确的做法如下表格所示:

数据连续读取方式
线程编号数据下标
00,512,... ...,
... ...... ...
511511,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;
}