设置并行块

add<<<N,1>>>( dev_a, dev_b, dev_c );
//N blocks x 1 thread/block = N parallel threads

这句话里面的1，就是the number of threads per block we want the CUDA runtime to create on our behalf，每个块的线程个数。

设置并行线程

add<<<1,N>>>( dev _ a, dev _ b, dev _ c );

我们把1和N反过来，就是一个N线程的单块程序。当我们获取索引的时候需要换一个表达来获取线程的索引：

混合设置

int tid = threadIdx.x + blockIdx.x * blockDim.x;

add<<< (N+127)/128, 128 >>>( dev _ a, dev _ b, dev _ c );

我们可以把thread per block给顾定成128，然后根据N来分配需要多少个block。(N+127)/128实际上就是求ceil的过程。

if (tid < N)
	c[tid] = a[tid] + b[tid];

因此，前面的历史遗留问题：判断小于N就被用作判断最后的block的这个thread是不是超过了我们问题所需要求解的范围N。

长向量相加

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

可以看到，在每一次运行完加法之后，都向后移动blockDim.x * gridDim.x个位置再做加法，以至超出范围，可以得到N以内而又超过硬件边界范围之外的内容以加入到这一次add的循环中进行累加。

处理图片

这个例子里生成了一个二维的block和一个二维的thread，每个block包含256个thread。

共享显存和同步

新增了一个关键词__shared__，这个关键词标注的变量能够给同一个block里的每一个thread都给单独整一个copy，但是其他block对这个shared变量不可读写。这种共享独取是非常高速的。但是，如果需要同一个块之内的线程交流，我们需要对交互的thread进行同步。

#include "../common/book.h"
#define imin(a,b) (a<b?a:b)
const int N = 33 * 1024;
const int threadsPerBlock = 256;
__global__ void dot( float *a, float *b, float *c ) {
	__shared__ float cache[threadsPerBlock];
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	int cacheIndex = threadIdx.x;
	float temp = 0;
	while (tid < N) {
		temp += a[tid] * b[tid];
		tid += blockDim.x * gridDim.x;
	}
	// set the cache values
	cache[cacheIndex] = temp;

这里面可以看到每一个thread都会进行乘法操作并把乘积存进cache[]中的对应位置（每一个thread都由自己在cache[]中对应的位置）。但是大家读写速度可能会不一样，需要同步一波：

// synchronize threads in this block
__syncthreads();

这句话能等大家都执行完才执行接下来的命令。

// for reductions, threadsPerBlock must be a power of 2
// because of the following code
int i = blockDim.x/2;
while (i != 0) {
	if (cacheIndex < i)
		cache[cacheIndex] += cache[cacheIndex + i];
	__syncthreads();
	i /= 2;
}

求和的时候，每一层都进行两个数的求和，因此求N个数字的和的时间复杂度为O(log n)。这个循环会在每一次循环中进行一次同步，每一次都会把待求和向量长度减半。单词迭代如下所示：
最后，第一个元素的值就是我们要的求和结果：

if (cacheIndex == 0)
	c[blockIdx.x] = cache[0];
}

这个cache[0]只是一个block中的product值，而我们需要计算的是所有的block中的求和值，这时由于维度很小，我们直接用CPU来计算。

const int blocksPerGrid =
	imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );
int main( void ) {
	float *a, *b, c, *partial_c;
	float *dev_a, *dev_b, *dev_partial_c;
	// allocate memory on the CPU side
	a = new float[N];
	b = new float[N];
	partial_c = new float[blocksPerGrid];
	// allocate the memory on the GPU
	HANDLE_ERROR( cudaMalloc( (void**)&dev_a,
		N*sizeof(float) ) );
	HANDLE_ERROR( cudaMalloc( (void**)&dev_b,
		N*sizeof(float) ) );
	HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c,
		blocksPerGrid*sizeof(float) ) );
	// fill in the host memory with data
	for (int i=0; i<N; i++) {
		a[i] = i;
		b[i] = i*2;
	}
	// copy the arrays 'a' and 'b' to the GPU
	HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float),
		cudaMemcpyHostToDevice ) );
	HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float),
		cudaMemcpyHostToDevice ) );
	dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a,
	dev_b,
	dev_partial_c );
	// copy the array 'c' back from the GPU to the CPU
	HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
		blocksPerGrid*sizeof(float),
		cudaMemcpyDeviceToHost ) );
	// finish up on the CPU side
	c = 0;
	for (int i=0; i<blocksPerGrid; i++) {
		c += partial_c[i];
	}
	#define sum_squares(x) (x*(x+1)*(2*x+1)/6)
	printf( “Does GPU value %.6g = %.6g?\n”, c,
		2 * sum_squares( (float)(N - 1) ) );
	// free memory on the GPU side
	cudaFree( dev_a );
	cudaFree( dev_b );
	cudaFree( dev_partial_c );
	// free memory on the CPU side
	free( a );
	free( b );
	free( partial_c );
}

同步语句

当__syncthreads()位于条件语句if内时，就会引发错误。同步命令指当所有thread都运行到这个同步时，才会执行下一步命令。因此，if中有的同步命令无法被执行，导致其他thread都不能继续运行了。例如:

if (cacheIndex < i) {
	cache[cacheIndex] += cache[cacheIndex + i];
	__syncthreads();
}

2020年9月7日