CUDAExample: 点积运算 Dot Product

引入的概念：共享内存

关键字： __share__ 

定义的变量将被存储在共享内存中。对于GPU上启动的每个线程块，CUDA C编译器都将创建该变量的一个副本。线程块中每个线程都将共享这块内存，但是线程却无法看到也不能修改其他线程块中的变量副本。这样就使得一个线程块中的多个线程能够在计算上进行通信和协作。共享内存缓冲区驻留在物理GPU上。因此，在访问共享你日常时的延迟要远远低于访问普通缓冲区的延迟，使得共享内存像每个线程块的攻速缓存或者中间结果暂存器那样高效。

如果想要在线程之间进行通信，还需要一种机制实现线程之间的同步。

#include "cpu_anim.h"
#include "cuda_runtime.h"
#include <device_launch_parameters.h>
#include <device_functions.h>

#define sum_squares(x) (x*(x+1)*(2*x+1)/6)

// 点乘的数据维度
const int N = 33 * 1024;
// 申请的网格、块内线程都是一维的
// 默认每个块中的线程申请512个
const int threadsPerBlock = 256;
// 默认申请的内存块个数不超过32，当前值为
const int blocksPerGrid = min(32, (N + threadsPerBlock - 1)/threadsPerBlock);

__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)
	{
		// 每个线程都计算数组a\b中相应元素乘积的总和 
		temp += a[tid] * b[tid];
		// printf("a = %f, b = %f\n", a[tid], b[tid]);

		// 如果申请的线程数少于序列的个数
		tid += blockDim.x * gridDim.x;
	}

	// 赋值，块中一个系列的线程计算出来的和
	cache[cacheIndex] = temp;
	// printf("cache[cacheIndex] = %f, cacheIndex = %d\n", cache[cacheIndex], cacheIndex);

	// 对线程块中的线程进行同步	
	__syncthreads();

	// 对于归约运算来说，以下代码要求 threadPerBlock必须是2的指数
	// 实现两两元素乘积的加和
	int offset = blockDim.x / 2;
	while (offset != 0)
	{
		if (cacheIndex < offset)
		{
			cache[cacheIndex] += cache[cacheIndex + offset];
		}
		__syncthreads();
		offset /= 2;
	}

	__syncthreads();
	if (cacheIndex == 0)
	{
		// 存储每个线程块的结果
		c[blockIdx.x] = cache[0];
		printf("sum  = %f, block index = %d\n", cache[0], blockIdx.x);
	}

}

int main()
{
	float* a,* b;
	a = new float[N];
	b = new float[N];
	// 每个线程块的结果
	float* partial_c;
	partial_c = new float[blocksPerGrid]; 
	float* dev_a, * dev_b, * dev_partial_c;
	
	cudaMalloc((void**)&dev_a, N * sizeof(float));
	cudaMalloc((void**)&dev_b, N * sizeof(float));
	cudaMalloc((void**)&dev_partial_c, blocksPerGrid * sizeof(float));

	// 初始化a、b序列
	for (size_t i = 0; i < N; i++)
	{
		a[i] = i;
		b[i] = i * 2;
	}

	// 将数组 a、b复制到GPU上
	cudaMemcpy(dev_a, a, N * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_b, b, N * sizeof(float), cudaMemcpyHostToDevice);
	// 核函数
	dot <<< blocksPerGrid, threadsPerBlock>>> (dev_a, dev_b, dev_partial_c);

	cudaMemcpy(partial_c, dev_partial_c, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);

	// 在CPU上完成最终的求和运算
	float sum = 0;
	for (size_t i = 0; i < blocksPerGrid; i++)
	{
		sum += partial_c[i];
	}

	std::cout << "GPU valuen = " << sum << " , real value = " << 2 * sum_squares((float)(N - 1)) << std::endl;
	getchar();
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_partial_c);

	free(a);
	free(b);
	free(partial_c);

	return 0;
}

输出的结果

关键字：

__syncthreads()

笔记：

 错误优化的情况是将第二个__syncthreads()放到了if()线程块中，那么任何cacheIndx大于或等于i的线程将永远都不能执行__syncthreads()。这将使处理器挂起，因为GPU在等在某个有永远都不会发生的事情。

• CUDA同步机制
  • cudaDeviceSychronize(): 用于主机代码
    • CPU等待正在运行的GPU任务结束
  • __syncthreads(): 用于设备代码
    • GPU中线程块内线程同步