CUDA初探：GPU的并行计算 - dflower

到目前为止，我们的程序并没有做什么有用的工作。所以，现在我们加入一个简单的动作，就是把一大堆数字，计算出它的平方和。

首先，把程序最前面的 include 部份改成：

#include <stdio.h>
#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;
    }
}

这个函式会产生一大堆 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);

上面这段程序会先呼叫 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;
}

在显示芯片上执行的程序有一些限制，例如它不能有传回值。其它的限制会在之后提到。

接下来是要让 CUDA 执行这个函式。在 CUDA 中，要执行一个函式，使用以下的语法：

    函式名称<<<block 数目, thread 数目, shared memory 大小>>>(参数

);

呼叫完后，还要把结果从显示芯片复制回主内存上。在 main 函式中加入以下的程序：

    sumOfSquares<<<1, 1, 0>>>(gpudata, result);

    int sum;
    cudaMemcpy(&sum, 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);

编译后执行，确认两个结果相同。

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;
}

把 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);

    sumOfSquares<<<1, 1, 0>>>(gpudata, result, time);

    int sum;
    clock_t time_used;
    cudaMemcpy(&sum, 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);

编译后执行，就可以看到执行所花费的时间了。

我们的第一个程序，并没有利用到任何并行化的功能。整个程序只有一个 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

接着，把 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;
}

程序里的 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(&sum, 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);

编译后执行，确认结果和原来相同。

这个版本的程序，在 GeForce 8800GT 上执行，只需要约 8.3M cycles，比前一版程序快了 77 倍！这就是透过大量 thread 来隐藏 latency 所带来的效果。

不过，如果计算一下它使用的内存带宽，就会发现其实仍不是很理想，大约只有 723MB/s 而已。这和 GeForce 8800GT 所具有的内存带宽是很大的差距。为什么会这样呢？

10.内存的存取模式

显卡上的内存是 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;
}

编译后执行，确认结果相同。

仅仅是这样简单的修改，实际执行的效率就有很大的差别。在 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 端要做的最后加总工作也会变多。

#include <stdio.h>

#include <stdlib.h>

#include <cuda_runtime.h>

#define DATA_SIZE 32768

#define THREAD_NUM 1

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(&prop,i) == cudaSuccess)

{

if(prop.major>=1)

{

//枚举详细信息

printf("Identify: %s\n",prop.name);

printf("Host Memory: %d\n",prop.canMapHostMemory);

printf("Clock Rate: %d khz\n",prop.clockRate);

printf("Compute Mode: %d\n",prop.computeMode);

printf("Device Overlap: %d\n",prop.deviceOverlap);

printf("Integrated: %d\n",prop.integrated);

printf("Kernel Exec Timeout Enabled: %d\n",prop.kernelExecTimeoutEnabled);

printf("Max Grid Size: %d * %d * %d\n",prop.maxGridSize[0],prop.maxGridSize[1],prop.maxGridSize[2]);

printf("Max Threads Dim: %d * %d * %d\n",prop.maxThreadsDim[0],prop.maxThreadsDim[1],prop.maxThreadsDim[2]);

printf("Max Threads per Block: %d\n",prop.maxThreadsPerBlock);

printf("Maximum Pitch: %d bytes\n",prop.memPitch);

printf("Minor Compute Capability: %d\n",prop.minor);

printf("Number of Multiprocessors: %d\n",prop.multiProcessorCount);

printf("32bit Registers Availble per Block: %d\n",prop.regsPerBlock);

printf("Shared Memory Available per Block: %d bytes\n",prop.sharedMemPerBlock);

printf("Alignment Requirement for Textures: %d\n",prop.textureAlignment);

printf("Constant Memory Available: %d bytes\n",prop.totalConstMem);

printf("Global Memory Available: %d bytes\n",prop.totalGlobalMem);

printf("Warp Size: %d threads\n",prop.warpSize);

printf("===================================\n");

break;

}

if(i==count)

{

fprintf(stderr,"There is no device supporting CUDA.\n");

return false;

}

cudaSetDevice(i);

return true;

}

void GenerateNums(int *numbers, int size)

{

for(int i=0;i<size;i++)

{

numbers[i]=rand()%10; //产生0-9的随机数

}

__global__ static void SumSquares(int *num, int *result, clock_t *time) //计算平方和线程

{

int sum=0;

int i;

clock_t start;

const int tid = threadIdx.x;

const int size = DATA_SIZE/THREAD_NUM;

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;

}

static void TestGPU()

{

int data[DATA_SIZE];

GenerateNums(data,DATA_SIZE);

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);

clock_t gpustart = clock();

SumSquares<<<1,THREAD_NUM,0>>>(gpudata,result,time);

int sum[THREAD_NUM];

clock_t time_used;

cudaMemcpy(&sum,result,sizeof(int)*THREAD_NUM,cudaMemcpyDeviceToHost);

cudaMemcpy(&time_used,time,sizeof(clock_t),cudaMemcpyDeviceToHost);

cudaFree(gpudata);

cudaFree(result);

cudaFree(time);

int final_sum=0;

for(int i=0;i<THREAD_NUM;i++)

{

final_sum+=sum[i];

}

clock_t gputime = clock()-gpustart;

printf("GPU sum: %d, cycle used %d, time used: %d\n",final_sum, time_used, gputime);

int sum1 = 0;

for(int i=0;i<DATA_SIZE;i++)

{

sum1 += data[i]*data[i];

}

printf("CPU sum: %d\n",sum1);

}

void main()

{

if(!InitCUDA())

{

getchar();

return;

}

printf("CUDA initialized.\n");

/////////////////////////////////

TestGPU();

getchar();

}

__global__ static void SumSquares(int *num, int *result, clock_t *time)  //计算平方和线程
{
    int sum=0;
    int i;
    clock_t start;
    const int tid = threadIdx.x;
    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;
}