引入Thread概念
第4章使用GPU做并行运算的例子，归结起来就是定义一个内核函数，将数组dev_a和dev_b某一对元素相加。GPU发起N个block运行内核函数。每个block有自己的索引，这样kernel就可以凭借这个索引区分自身，来计算数组对应的元素。在这些例子中，每个block只创建了一个线程。
本章开始引入Thread概念，实际上一个block可以创建多个Thread。为了方便说明，前面将两个数组对应元素相加的运算，之后开始用矢量相加来指称。数组a为矢量a的各个分量，数组b为矢量b的各个分量。之前<<<N, 1>>> 的第二个参数一直没有做介绍，现在我们可以猜到，它就是用来指定一个block里创建多少个线程。
既然这样，我们试试创建一个block，block里创建N个线程来完成上一章的矢量相加的例子。main函数中，调用add的方式从
add<<<N, 1>>>(dev_a, dev_b, dev_c);
改为
add<<<1, N>>>(dev_a, dev_b, dev_c);
内核函数add也需要做个很小的修改。之前各个block通过blockIdx.x作为索引，对矢量a，b的分量求和，现在使用同一个block的各个线程来计算各个分量了，索引肯定不能再用blockIdx了。CUDA贴心的送上threadIdx。
int tid = threadIdx.x;
这个代码就只改了这两个地方。
#include "../common/book.h"
#define N 10
__global__ void add( int *a, int *b, int *c ) {
// int tid = blockIdx.x;
int tid = threadIdx.x;
if (tid < N)
c[tid] = a[tid] + b[tid];
}
int main( void ) {
int a[N], b[N], c[N];
int *dev_a, *dev_b, *dev_c;
// allocate the memory on the GPU
HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_c, N * sizeof(int) ) );
// fill the arrays 'a' and 'b' on the CPU
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = i * i;
}
// copy the arrays 'a' and 'b' to the GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(int),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(int),
cudaMemcpyHostToDevice ) );
// add<<<N,1>>>( dev_a, dev_b, dev_c );
add<<<1,N>>>( dev_a, dev_b, dev_c );
// copy the array 'c' back from the GPU to the CPU
HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(int),
cudaMemcpyDeviceToHost ) );
// display the results
for (int i=0; i<N; i++) {
printf( "%d + %d = %d\n", a[i], b[i], c[i] );
}
// free the memory allocated on the GPU
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_c ) );
return 0;
}Block和Thread的限制
受限于硬件资源，GPU能够支持的block数目和每个block里能包含的thread数目都是有限的。书中的说法是Block的最大数目是65535，每个block最多thread数可以通过设备的属性maxThreadsPerBlock查到。2010年最好GPU，这个数字是512。2024年的今天，这个数字也只翻了一倍1024. 无论如何，我们肯定会遇到元素个数大于最大Block数或者maxThreadsPerBlock，甚至两者的乘积65,535 * 512 = 33,553,920。
如何突破这一限制？
首先，我们不能再使用<<<N,1>>>和<<<1,N>>>，而要为block数和每Block里thread数选择合适的值。这里不光有上面所说的原因，还牵涉到CUDA的一个机制，属于同一个block的thread可以共享内存。这个内存跟cudaMalloc分配的内存是有区别的，cudaMalloc分配的内存是从GPU里DDR，甚至是供GPU使用的主板DDR。而给block的共享内存要快得多。后面章节会详细使用方法，这里我们先学习block，thread组合在一起用法。
其次，当元素个数大于限制时，分多次处理。
为了更好的说明这个问题，我们先认识一下block和thread的逻辑结构。
以二维为例（比较好画图，实际上最高支持3维的），通过<<<gridDim, blockDim>>>可以在GPU端生成指定数量的线程，线程以上面的方式进行组合管理。
gridDim指定x，y方向线程块的数量，blockDim指定x，y方向thread的数量。当我们只需要一维来处理，那么y方向设置为1，在图1中对应就是蓝色方框显示的部分。我们用一维来处理矢量相加，二维来处理图像。
那么对于图中Block(1,0)里的Thread(1,0)，它来计算矢量的哪个分量呢？在本block中，它前面有1个线程，同时前面还有一个block即block(0,0)里的2个线程，所以它应该处理分量的索引应该是 1 + 1*2 = 3。由此对于任意thread，它要处理的分量的所有计算方法是：
int tid = threadIdx.x + blockIdx.x * blockDim.x
这样，我们解决组合使用block和thread时，索引的计算问题。接下来我们来看，当元素个数超出block/thread组合一次性处理能力时该怎么办。
图1中，一次性只能处理6个元素。超过6个元素，比如说当N=15时，thread(1,0)的任务是计算dev_a[3] + dev_c[3]。最右边的线程thread(2,1)计算dev_a[5] + dev_b[5]。数组再往后，就没有剩余线程可以用了。这时，我们可以让最左边的线程计算完dev_a[0] + dev_b[0]后，来计算dev_a[6]+dev_b[6], 然后再计算dev_a[12] + dev_b[12]。也就是说完成最初计算所得索引的元素的相加之后，跳过 gridDim.x * blockDim.x个元素后， 继续把对应的两个元素相加，直到超出数组大小。
我们的内核函数将修改为：
__global__ void add( int *a, int *b, int *c ) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += blockDim.x * gridDim.x;
}
}长矢量相加
下面是一个长矢量相加的例子，在随书代码的chapter05\add_loop_long_blocks.cu。分量的个数N 为 33*1024。有了前面的经验，看起来应该没有什么难度。要注意是N比较大，所以CPU这边的a，b，c也改成动态分配了。
#include "../common/book.h"
#define N (33 * 1024)
__global__ void add( int *a, int *b, int *c ) {
int tid = threadIdx.x + blockIdx.x * blockDim.x;
while (tid < N) {
c[tid] = a[tid] + b[tid];
tid += blockDim.x * gridDim.x;
}
}
int main( void ) {
int *a, *b, *c;
int *dev_a, *dev_b, *dev_c;
// allocate the memory on the CPU
a = (int*)malloc( N * sizeof(int) );
b = (int*)malloc( N * sizeof(int) );
c = (int*)malloc( N * sizeof(int) );
// allocate the memory on the GPU
HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N * sizeof(int) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_c, N * sizeof(int) ) );
// fill the arrays 'a' and 'b' on the CPU
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = 2 * i;
}
// copy the arrays 'a' and 'b' to the GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N * sizeof(int),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N * sizeof(int),
cudaMemcpyHostToDevice ) );
add<<<128,128>>>( dev_a, dev_b, dev_c );
// copy the array 'c' back from the GPU to the CPU
HANDLE_ERROR( cudaMemcpy( c, dev_c, N * sizeof(int),
cudaMemcpyDeviceToHost ) );
// verify that the GPU did the work we requested
bool success = true;
for (int i=0; i<N; i++) {
if ((a[i] + b[i]) != c[i]) {
printf( "Error: %d + %d != %d\n", a[i], b[i], c[i] );
success = false;
}
}
if (success) printf( "We did it!\n" );
// free the memory we allocated on the GPU
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_c ) );
// free the memory we allocated on the CPU
free( a );
free( b );
free( c );
return 0;
}
