[转载]http://www.pcinlife.com/article/graphics/2008-06-04/1212575164d532.html
第一个CUDA程序
第一个CUDA程序
CUDA 目前有两种不同的 API:Runtime API 和 Driver API,两种 API 各有其适用的范围。由于 runtime API 较容易使用,一开始我们会以 runetime API 为主。
CUDA 的初始化
|
要使用 runtime API 的时候,需要 include
cuda_runtime.h
。所以,在程序的最前面,加上
#include <stdio.h>
#include <cuda_runtime.h>
#include <cuda_runtime.h>
接下来是一个
InitCUDA
函式,会呼叫 runtime API 中,有关初始化 CUDA 的功能:
bool InitCUDA()
{
int count;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(∝, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
cudaSetDevice(i);
return true;
}
{
int count;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(∝, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
cudaSetDevice(i);
return true;
}
这个函式会先呼叫
cudaGetDeviceCount
函式,取得支持 CUDA 的装置的数目。如果系统上没有支持 CUDA 的装置,则它会传回 1,而 device 0 会是一个仿真的装置,但不支持 CUDA 1.0 以上的功能。所以,要确定系统上是否有支持 CUDA 的装置,需要对每个 device 呼叫
cudaGetDeviceProperties
函式,取得装置的各项数据,并判断装置支持的 CUDA 版本(prop.major 和 prop.minor 分别代表装置支持的版本号码,例如 1.0 则 prop.major 为 1 而 prop.minor 为 0)。
透过 cudaGetDeviceProperties 函式可以取得许多数
据,除了装置支持的 CUDA 版本之外,还有装置的名称、内存的大小、最大的 thread 数目、执行单元的频率等等。详情可参考 NVIDIA 的 CUDA Programming Guide。
在找到支持 CUDA 1.0 以上的装置之后,就可以呼叫
cudaSetDevice
函式,把它设为目前要使用的装置。
最后是 main 函式。在 main 函式中我们直接呼叫刚才的 InitCUDA 函式,并显示适当的讯息:
int main()
{
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
return 0;
}
{
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
return 0;
}
这样就可以利用 nvcc 来 compile 这个程序了。使用 Visual Studio 的话,若按照先前的设定方式,可以直接 Build Project 并执行。
nvcc 是 CUDA 的 compile 工具,它会将 .cu 檔拆解出在 GPU 上执行的部份,及在 host 上执行的部份,并呼叫适当的程序进行 compile 动作。在 GPU 执行的部份会透过 NVIDIA 提供的 compiler 编译成中介码,而 host 执行的部份则会透过系统上的 C++ compiler 编译(在 Windows 上使用 Visual C++ 而在 Linux 上使用 gcc)。
编译后的程序,执行时如果系统上有支持 CUDA 的装置,应该会显示 CUDA initialized. 的讯息,否则会显示相关的错误讯息。
利用 CUDA 进行运算
|
到目前为止,我们的程序并没有做什么有用的工作。所以,现在我们加入一个简单的动作,就是把一大堆数字,计算出它的平方和。
首先,把程序最前面的 include 部份改成:
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#define DATA_SIZE 1048576
int data[DATA_SIZE];
#include <stdlib.h>
#include <cuda_runtime.h>
#define DATA_SIZE 1048576
int data[DATA_SIZE];
并加入一个新函式
GenerateNumbers
:
void GenerateNumbers(int *number, int size)
{
for(int i = 0; i < size; i++) {
number[i] = rand() % 10;
}
}
{
for(int i = 0; i < size; i++) {
number[i] = rand() % 10;
}
}
这个函式会产生一大堆 0 ~ 9 之间的随机数。
要利用 CUDA 进行计算之前,要先把数据复制到显卡内存中,才能让显示芯片使用。因此,需要取得一块适当大小的显卡内存,再把产生好的数据复制进去。在
main
函式中加入:
GenerateNumbers(data, DATA_SIZE);
int* gpudata, *result;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
int* gpudata, *result;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
上面这段程序会先呼叫
GenerateNumbers
产生随机数,并呼叫
cudaMalloc
取得一块显卡内存(
result
则是用来存取计算结果,在稍后会用到),并透过
cudaMemcpy
将产生的随机数复制到显卡内存中。cudaMalloc 和 cudaMemcpy 的用法和一般的 malloc 及 memcpy 类似,不过 cudaMemcpy 则多出一个参数,指示复制内存的方向。在这里因为是从主内存复制到显卡内存,所以使用
cudaMemcpyHostToDevice
。如果是从显卡内存到主内存,则使用
cudaMemcpyDeviceToHost
。这在之后会用到。
接下来是要写在显示芯片上执行的程序。在 CUDA 中,在函式前面加上
__global__
表示这个函式是要在显示芯片上执行的。因此,加入以下的函式:
__global__ static void sumOfSquares(int *num, int* result)
{
int sum = 0;
int i;
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
}
{
int sum = 0;
int i;
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
}
在显示芯片上执行的程序有一些限制,例如它不能有传回值。其它的限制会在之后提到。
接下来是要让 CUDA 执行这个函式。在 CUDA 中,要执行一个函式,使用以下的语法:
函式名称<<<block 数目, thread 数目, shared memory 大小>>>(参数...);
呼叫完后,还要把结果从显示芯片复制回主内存上。在
main
函式中加入以下的程序:
sumOfSquares<<<1, 1, 0>>>(gpudata, result);
int sum;
cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("sum: %d\n", sum);
int sum;
cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("sum: %d\n", sum);
因为这个程序只使用一个 thread,所以 block 数目、thread 数目都是 1。我们也没有使用到任何 shared memory,所以设为 0。编译后执行,应该可以看到执行的结果。
为了确定执行的结果正确,我们可以加上一段以 CPU 执行的程序代码,来验证结果:
sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", sum);
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", sum);
编译后执行,确认两个结果相同。
计算运行时间
|
CUDA 提供了一个
clock
函式,可以取得目前的 timestamp,很适合用来判断一段程序执行所花费的时间(单位为 GPU 执行单元的频率)。这对程序的优化也相当有用。要在我们的程序中记录时间,把
sumOfSquares
函式改成:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
int sum = 0;
int i;
clock_t start = clock();
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
*time = clock() - start;
}
clock_t* time)
{
int sum = 0;
int i;
clock_t start = clock();
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
*time = clock() - start;
}
把
main
函式中间部份改成:
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, 1, 0>>>(gpudata, result, time);
int sum;
clock_t time_used;
cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
int sum;
clock_t time_used;
cudaMemcpy(∑, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("sum: %d time: %d\n", sum, time_used);
编译后执行,就可以看到执行所花费的时间了。
如果计算实际运行时间的话,可能会注意到它的执行效率并不好。这是因为我们的程序并没有利用到 CUDA 的主要的优势,即并行化执行。在下一段文章中,会讨论如何进行优化的动作。
在
上一篇文章
中, 我们做了一个计算一大堆数字的平方和的程序。不过,我们也提到这个程序的执行效率并不理想。当然,实际上来说,如果只是要做计算平方和的动作,用 CPU 做会比用 GPU 快得多。这是因为平方和的计算并不需要太多运算能力,所以几乎都是被内存带宽所限制。因此,光是把数据复制到显卡内存上的这个动作,所需要的时间,可能已 经和直接在 CPU 上进行计算差不多了。
不过,如果进行平方和的计算,只是一个更复杂的计算过程的一部份的话,那么当然在 GPU 上计算还是有它的好处的。而且,如果数据已经在显卡内存上(例如在 GPU 上透过某种算法产生),那么,使用 GPU 进行这样的运算,还是会比较快的。
刚才也提到了,由于这个计算的主要瓶颈是内存带宽,所以,理论上显卡的内存带宽是相当大的。这里我们就来看看,倒底我们的第一个程序,能利用到多少内存带宽。
程序的并行化
|
我们的第一个程序,并没有利用到任何并行化的功能。整个程序只有一个 thread。在 GeForce 8800GT 上面,在 GPU 上执行的部份(称为 "
kernel
") 大约花费 640M 个频率。GeForce 8800GT 的执行单元的频率是 1.5GHz,因此这表示它花费了约 0.43 秒的时间。1M 个 32 bits 数字的数据量是 4MB,因此,这个程序实际上使用的内存带宽,只有 9.3MB/s 左右!这是非常糟糕的表现。
为什么会有这样差的表现呢?这是因为 GPU 的架构特性所造成的。在 CUDA 中,一般的数据复制到的显卡内存的部份,称为
global memory
。 这些内存是没有 cache 的,而且,存取 global memory 所需要的时间(即 latency)是非常长的,通常是数百个 cycles。由于我们的程序只有一个 thread,所以每次它读取 global memory 的内容,就要等到实际读取到数据、累加到 sum 之后,才能进行下一步。这就是为什么它的表现会这么的差。
由于 global memory 并没有 cache,所以要避开巨大的 latency 的方法,就是要利用大量的 threads。假设现在有大量的 threads 在同时执行,那么当一个 thread 读取内存,开始等待结果的时候,GPU 就可以立刻切换到下一个 thread,并读取下一个内存位置。因此,理想上当 thread 的数目够多的时候,就可以完全把 global memory 的巨大 latency 隐藏起来了。
要怎么把计算平方和的程序并行化呢?最简单的方法,似乎就是把数字分成若干组,把各组数字分别计算平方和后,最后再把每组的和加总起来就可以了。一开始,我们可以把最后加总的动作,由 CPU 来进行。
首先,在 first_cuda.cu 中,在
#define DATA_SIZE
的后面增加一个 #define,设定 thread 的数目:
#define DATA_SIZE 1048576
#define THREAD_NUM 256
#define THREAD_NUM 256
接着,把 kernel 程序改成:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int size = DATA_SIZE / THREAD_NUM;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid * size; i < (tid + 1) * size; i++) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
clock_t* time)
{
const int tid = threadIdx.x;
const int size = DATA_SIZE / THREAD_NUM;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid * size; i < (tid + 1) * size; i++) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
程序里的
threadIdx
是 CUDA 的一个内建的变量,表示目前的 thread 是第几个 thread(由 0 开始计算)。以我们的例子来说,会有 256 个 threads,所以同时会有 256 个
sumOfSquares
函式在执行,但每一个的
threadIdx.x
则分别会是 0 ~ 255。利用这个变量,我们就可以让每一份函式执行时,对整个数据不同的部份计算平方和。另外,我们也让计算时间的动作,只在 thread 0(即
threadIdx.x
= 0 的时候)进行。
同样的,由于会有 256 个计算结果,所以原来存放
result
的内存位置也要扩大。把
main
函式中的中间部份改成:
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
int sum[THREAD_NUM];
clock_t time_used;
cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
int sum[THREAD_NUM];
clock_t time_used;
cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
可以注意到我们在呼叫
sumOfSquares
函式时,指定
THREAD_NUM
为 thread 的数目。最后,在 CPU 端把计算好的各组数据的平方和进行加总:
int final_sum = 0;
for(int i = 0; i < THREAD_NUM; i++) {
final_sum += sum[i];
}
printf("sum: %d time: %d\n", final_sum, time_used);
final_sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", final_sum);
for(int i = 0; i < THREAD_NUM; i++) {
final_sum += sum[i];
}
printf("sum: %d time: %d\n", final_sum, time_used);
final_sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", final_sum);
编译后执行,确认结果和原来相同。
这个版本的程序,在 GeForce 8800GT 上执行,只需要约 8.3M cycles,比前一版程序快了 77 倍!这就是透过大量 thread 来隐藏 latency 所带来的效果。
不过,如果计算一下它使用的内存带宽,就会发现其实仍不是很理想,大约只有 723MB/s 而已。这和 GeForce 8800GT 所具有的内存带宽是很大的差距。为什么会这样呢?
内存的存取模式
|
显卡上的内存是 DRAM,因此最有效率的存取方式,是以连续的方式存取。前面的程序,虽然看起来是连续存取内存位置(每个 thread 对一块连续的数字计算平方和),但是我们要考虑到实际上 thread 的执行方式。前面提过,当一个 thread 在等待内存的数据时,GPU 会切换到下一个 thread。也就是说,实际上执行的顺序是类似
thread 0 -> thread 1 -> thread 2 -> ...
因此,在同一个 thread 中连续存取内存,在实际执行时反而不是连续了。要让实际执行结果是连续的存取,我们应该要让 thread 0 读取第一个数字,thread 1 读取第二个数字…依此类推。所以,我们可以把 kernel 程序改成如下:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
clock_t* time)
{
const int tid = threadIdx.x;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
编译后执行,确认结果相同。
仅仅是这样简单的修改,实际执行的效率就有很大的差别。在 GeForce 8800GT 上,上面的程序执行需要的频率是 2.6M cycles,又比前一版程序快了三倍。不过,这样仍只有 2.3GB/s 的带宽而已。
这是因为我们使用的 thread 数目还是不够多的原因。理论上 256 个 threads 最多只能隐藏 256 cycles 的 latency。但是 GPU 存取 global memory 时的 latency 可能高达 500 cycles 以上。如果增加 thread 数目,就可以看到更好的效率。例如,可以把
THREAD_NUM
改成 512。在 GeForce 8800GT 上,这可以让执行花费的时间减少到 1.95M cycles。有些改进,但是仍不够大。不幸的是,目前 GeForce 8800GT 一个 block 最多只能有 512 个 threads,所以不能再增加了,而且,如果 thread 数目增加太多,那么在 CPU 端要做的最后加总工作也会变多。
更多的并行化
|
前面提到了 block。在之前介绍呼叫 CUDA 函式时,也有提到 "block 数目" 这个参数。到目前为止,我们都只使用一个 block。究竟 block 是什么呢?
在 CUDA 中,thread 是可以分组的,也就是 block。一个 block 中的 thread,具有一个共享的 shared memory,也可以进行同步工作。不同 block 之间的 thread 则不行。在我们的程序中,其实不太需要进行 thread 的同步动作,因此我们可以使用多个 block 来进一步增加 thread 的数目。
首先,在
#define DATA_SIZE
的地方,改成如下:
#define DATA_SIZE 1048576
#define BLOCK_NUM 32
#define THREAD_NUM 256
#define BLOCK_NUM 32
#define THREAD_NUM 256
这表示我们会建立 32 个 blocks,每个 blocks 有 256 个 threads,总共有 32*256 = 8192 个 threads。
接着,我们把 kernel 部份改成:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int sum = 0;
int i;
if(tid == 0) time[bid] = clock();
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
sum += num[i] * num[i];
}
result[bid * THREAD_NUM + tid] = sum;
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
clock_t* time)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int sum = 0;
int i;
if(tid == 0) time[bid] = clock();
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
sum += num[i] * num[i];
}
result[bid * THREAD_NUM + tid] = sum;
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
blockIdx.x
和
threadIdx.x
一样是 CUDA 内建的变量,它表示的是目前的 block 编号。另外,注意到我们把计算时间的方式改成每个 block 都会记录开始时间及结束时间。
main 函式部份,修改成:
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result,
sizeof(int) * THREAD_NUM * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM, 0>>>(gpudata, result,
time);
int sum[THREAD_NUM * BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) {
final_sum += sum[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("sum: %d time: %d\n", final_sum, max_end - min_start);
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result,
sizeof(int) * THREAD_NUM * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM, 0>>>(gpudata, result,
time);
int sum[THREAD_NUM * BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(∑, result, sizeof(int) * THREAD_NUM * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) {
final_sum += sum[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("sum: %d time: %d\n", final_sum, max_end - min_start);
基本上我们只是把 result 的大小变大,并修改计算时间的方式,把每个 block 最早的开始时间,和最晚的结束时间相减,取得总运行时间。
这个版本的程序,执行的时间减少很多,在 GeForce 8800GT 上只需要约 150K cycles,相当于 40GB/s 左右的带宽。不过,它在 CPU 上执行的部份,需要的时间加长了(因为 CPU 现在需要加总 8192 个数字)。为了避免这个问题,我们可以让每个 block 把自己的每个 thread 的计算结果进行加总。
Thread 的同步
|
前面提过,一个 block 内的 thread 可以有共享的内存,也可以进行同步。我们可以利用这一点,让每个 block 内的所有 thread 把自己计算的结果加总起来。把 kernel 改成如下:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
result[bid] = shared[0];
}
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
result[bid] = shared[0];
}
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
利用
__shared__
声明的变量表示这是 shared memory,是一个 block 中每个 thread 都共享的内存。它会使用在 GPU 上的内存,所以存取的速度相当快,不需要担心 latency 的问题。
__syncthreads()
是一个 CUDA 的内部函数,表示 block 中所有的 thread 都要同步到这个点,才能继续执行。在我们的例子中,由于之后要把所有 thread 计算的结果进行加总,所以我们需要确定每个 thread 都已经把结果写到 shared[tid] 里面了。
接下来,把
main
函式的一部份改成:
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM,
THREAD_NUM * sizeof(int)>>>(gpudata, result, time);
int sum[BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(∑, result, sizeof(int) * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < BLOCK_NUM; i++) {
final_sum += sum[i];
}
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<BLOCK_NUM, THREAD_NUM,
THREAD_NUM * sizeof(int)>>>(gpudata, result, time);
int sum[BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(∑, result, sizeof(int) * BLOCK_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < BLOCK_NUM; i++) {
final_sum += sum[i];
}
可以注意到,现在 CPU 只需要加总
BLOCK_NUM
也就是 32 个数字就可以了。
不过,这个程序由于在 GPU 上多做了一些动作,所以它的效率会比较差一些。在 GeForce 8800GT 上,它需要约 164K cycles。
当然,效率会变差的一个原因是,在这一版的程序中,最后加总的工作,只由每个 block 的 thread 0 来进行,但这并不是最有效率的方法。理论上,把 256 个数字加总的动作,是可以并行化的。最常见的方法,是透过树状的加法:
![](https://i-blog.csdnimg.cn/blog_migrate/ba37d2dcdf52c99a75a1424552bea2c8.png)
把 kernel 改成如下:
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
int offset = 1, mask = 1;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
while(offset < THREAD_NUM) {
if((tid & mask) == 0) {
shared[tid] += shared[tid + offset];
}
offset += offset;
mask = offset + mask;
__syncthreads();
}
if(tid == 0) {
result[bid] = shared[0];
time[bid + BLOCK_NUM] = clock();
}
}
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
int offset = 1, mask = 1;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
while(offset < THREAD_NUM) {
if((tid & mask) == 0) {
shared[tid] += shared[tid + offset];
}
offset += offset;
mask = offset + mask;
__syncthreads();
}
if(tid == 0) {
result[bid] = shared[0];
time[bid + BLOCK_NUM] = clock();
}
}
后面的 while 循环就是进行树状加法。
main
函式则不需要修改。
这一版的程序,在 GeForce 8800GT 上执行需要的时间,大约是 140K cycles(相当于约 43GB/s),比完全不在 GPU 上进行加总的版本还快!这是因为,在完全不在 GPU 上进行加总的版本,写入到 global memory 的数据数量很大(8192 个数字),也对效率会有影响。所以,这一版程序不但在 CPU 上的运算需求降低,在 GPU 上也能跑的更快。
进一步改善
|
上一个版本的树状加法是一般的写法,但是它在 GPU 上执行的时候,会有 share memory 的 bank conflict 的问题(详情在后面介绍 GPU 架构时会提到)。采用下面的方法,可以避免这个问题:
offset = THREAD_NUM / 2;
while(offset > 0) {
if(tid < offset) {
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
while(offset > 0) {
if(tid < offset) {
shared[tid] += shared[tid + offset];
}
offset >>= 1;
__syncthreads();
}
这样同时也省去了
mask
变数。因此,这个版本的执行的效率就可以再提高一些。在 GeForce 8800GT 上,这个版本执行的时间是约 137K cycles。当然,这时差别已经很小了。如果还要再提高效率,可以把树状加法整个展开:
if(tid < 128) { shared[tid] += shared[tid + 128]; }
__syncthreads();
if(tid < 64) { shared[tid] += shared[tid + 64]; }
__syncthreads();
if(tid < 32) { shared[tid] += shared[tid + 32]; }
__syncthreads();
if(tid < 16) { shared[tid] += shared[tid + 16]; }
__syncthreads();
if(tid < 8) { shared[tid] += shared[tid + 8]; }
__syncthreads();
if(tid < 4) { shared[tid] += shared[tid + 4]; }
__syncthreads();
if(tid < 2) { shared[tid] += shared[tid + 2]; }
__syncthreads();
if(tid < 1) { shared[tid] += shared[tid + 1]; }
__syncthreads();
__syncthreads();
if(tid < 64) { shared[tid] += shared[tid + 64]; }
__syncthreads();
if(tid < 32) { shared[tid] += shared[tid + 32]; }
__syncthreads();
if(tid < 16) { shared[tid] += shared[tid + 16]; }
__syncthreads();
if(tid < 8) { shared[tid] += shared[tid + 8]; }
__syncthreads();
if(tid < 4) { shared[tid] += shared[tid + 4]; }
__syncthreads();
if(tid < 2) { shared[tid] += shared[tid + 2]; }
__syncthreads();
if(tid < 1) { shared[tid] += shared[tid + 1]; }
__syncthreads();
前面介绍的计算平方和的程序,似乎没有什么实用价值。所以我们的第二个 CUDA 程序,要做一个确实有(某些)实用价值的程序,也就是进行矩阵乘法。而且,这次我们会使用浮点数。
虽然矩阵乘法有点老套,不过因为它相当简单,而且也可以用来介绍一些有关 CUDA 的有趣性质。