内存的存取模式
显卡上的内存是 DRAM，因此最有效率的存取方式，是以连续的方式存取。前面的程序，虽然看起来是连续存取内存位置(每个 thread 对一块连续的数字计算平方和)，但是我们要考虑到实际上 thread 的执行方式。前面提过，当一个 thread 在等待内存的数据时，GPU 会切换到下一个 thread。也就是说，实际上执行的顺序是类似
thread 0 -> thread 1 -> thread 2 -> ...
因此，在同一个 thread 中连续存取内存，在实际执行时反而不是连续了。要让实际执行结果是连续的存取，我们应该要让 thread 0 读取第一个数字，thread 1 读取第二个数字…依此类推。所以，我们可以把 kernel 程序改成如下：
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 BLOCK_NUM 32
#define THREAD_NUM 256
这表示我们会建立 32 个 blocks，每个 blocks 有 256 个 threads，总共有 32*256 = 8192 个 threads。
接着，我们把 kernel 部份改成：
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 函式部份，修改成：
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(&sum, 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 的计算结果进行加总。
