CUDA向量点乘

两个向量的点乘是重要的数学运算，也将会解释CUDA编程中的一个重要概念:归约运算。它和之前的元素两两相加的向量加法操作很类似。不同的是你需要将元素两两相乘。线程需要将它们的所有单个乘法结果连续累加起来，因为所有的一对对的乘法结果需要被累加起来,才能得到点乘的最终结果。最终的点乘的结果将是一个单一值。这种原始输入是两个数组而输出却缩减为一个(单一值)的运算，在CUDA里叫作归约运算。归约运算在很多应用程序里都有用。

#include "stdio.h"
#include<iostream>
#include <cuda.h>
#include <cuda_runtime.h>

#define N 1024
#define threadsPerBlock 512

__global__ void gpu_dot(float *d_a, float *d_b, float *d_c)
{
	//Declare shared memory
	__shared__ float partial_sum[threadsPerBlock];
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	//Calculate index for shared memory 
	int index = threadIdx.x;
	//Calculate Partial Sum
	float sum = 0;
	while (tid < N) 
	{
		sum += d_a[tid] * d_b[tid];
		tid += blockDim.x * gridDim.x;
	}

	// Store partial sum in shared memory
	partial_sum[index] = sum;

	// synchronize threads 
	__syncthreads();

	// Calculating partial sum for whole block in reduce operation
	int i = blockDim.x / 2;
	while (i != 0) 
	{
		if (index < i)
			partial_sum[index] += partial_sum[index + i];
		__syncthreads();
		i /= 2;
	}
	//Store block partial sum in global memory
	if (index == 0)
		d_c[blockIdx.x] = partial_sum[0];
}

int main() 
{
	//Declare Host Array
	float *h_a, *h_b, h_c, *partial_sum;
	//Declare device Array
	float *d_a, *d_b, *d_partial_sum;
	//Calculate total number of blocks per grid
	int block_calc = (N + threadsPerBlock - 1) / threadsPerBlock;
	int blocksPerGrid = (32 < block_calc ? 32 : block_calc);
	// allocate memory on the host side
	h_a = (float*)malloc(N * sizeof(float));
	h_b = (float*)malloc(N * sizeof(float));
	partial_sum = (float*)malloc(blocksPerGrid * sizeof(float));

	// allocate the memory on the device
	cudaMalloc((void**)&d_a, N * sizeof(float));
	cudaMalloc((void**)&d_b, N * sizeof(float));
	cudaMalloc((void**)&d_partial_sum, blocksPerGrid * sizeof(float));

	// fill the host array with data
	for (int i = 0; i < N; i++) 
	{
		h_a[i] = i;
		h_b[i] = 2;
	}

	cudaMemcpy(d_a, h_a, N * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_b, h_b, N * sizeof(float), cudaMemcpyHostToDevice);
	//Call kernel 
	gpu_dot << <blocksPerGrid, threadsPerBlock >> >(d_a, d_b, d_partial_sum);

	// copy the array back to host memory
	cudaMemcpy(partial_sum, d_partial_sum, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost);

	// Calculate final dot product on host
	h_c = 0;
	for (int i = 0; i<blocksPerGrid; i++)
	{
		h_c += partial_sum[i];
	}
	printf("The computed dot product is: %f\n", h_c);
#define cpu_sum(x) (x*(x+1))
	if (h_c == cpu_sum((float)(N - 1)))
	{
		printf("The dot product computed by GPU is correct\n");
	}
	else
	{
		printf("Error in dot product computation");
	}
	
	// free memory on host and device
	cudaFree(d_a);
	cudaFree(d_b);
	cudaFree(d_partial_sum);
	free(h_a);
	free(h_b);
	free(partial_sum);
}

该内核函数使用两个数组作为输入（参数中的a，b)，并将最终得到的部分和放入第三
个数组（参数中的c)，然后在内核内部用共享内存来存储中间的部分和计算结果。具体的共享内存中的元素个数等于每个块的线程数，因为每个块都将有单独的一份这个共享内存的副本。(内核中)定义共享内存后，我们计算出来两个索引值:第一个索引值tid用来计算唯一的线程ID，第二个索引值index变量计算为线程在块内部中的局部ID，后面的归约运算中会用到。再次强调:每个块都有单独的一份共享内存副本，所以每个线程ID索引到的共享内存只能是当前块自己的那个副本。
再往下的 while循环将会对通过线程ID ( tid变量)索引读取每对元素，将它们相乘并累加到sum变量上。当线程总数小于元素数量的时候，它也会循环将tid索引累加偏移到当前线程总数，继续索引下一对元素，并进行运算。每个线程得到的部分和结果被写入到共享内存。我们将继续使用共享内存上的这些线程的部分和计算出当前块的总体部分和。所以，在下面归约运算我们对这些共享内存中的数据进行读取之前，必须确保每个线程都已经完成了对共享内存的写入。可以通过__syncthreads()同步函数做到这一点。
现在，一种计算当前块的部分和的方法，就是让一个线程串行循环将这些所有线程的部分和进行累加，这样就可以得到最终当前块的部分和。只要一个线程就能完成这个归约运算，只是这将需要N次运算，N等于这些部分和的个数(又等于块中的线程数量)。
问题来了，我们就不能并行完成这个归约运算吗?答案是肯定的!并行化的点子就是，每个线程累加2个数的操作，并将每个线程得到的1个数结果覆盖写入这两个数中第1个数的位置。因为每个线程都累加了2个数，因此可以在第一个数中完成操作。现在，对剩余的部分我们重复这个过程，直到最终得到当前块的整体部分和。算法的复杂度是log(N),这远比我们之前的单个线程串行归约累加所需要的N复杂度好得多。
这种并行归约运算是通过条件为( i!=0 ）的while循环进行的。当前块中每个满足一定条件的线程每次循环累加当前它的ID索引的部分和，和ID + blockDim/2作为索引的部分和，并将结果回写共享内存。重复这个过程直到得到最终的单一结果，即整个块最终的一个部分和。最终每个块的单一部分和结果被写入到共享内存中(的特定位置)。通过块的唯一ID进行索引，每个块能确定这个单独的写入特定位置。然而，我们还是没有得到最终结果。这个最后的结果通过设备上的函数，或者(主机上）的 main 函数执行都可以。
我们定义了3种数组来容纳输入数据和输出数据，并分别为它们在主机上和设备上分配空间。通过循环，对主机上的a， b两个数组进行初始化。其中的一个元素被初始化为О到N的值，另一个里面的元素值都恒定为2。下一步则是计算Grid中的块数量和每个块中的线程数，就像我们在本章开头计算过的那样。记住，你也可以简单地将这两个值都设定为常数，就像本章的第一个程序以避免繁杂。
将这两个输入数组复制到显存，并将它们作为内核函数的参数。内核函数会将每个块计算得到的部分和写人到对应每个块ID位置的结果数组的对应元素中。然后该结果数组被复制回主机上的 partial_sum数组中。最终的点乘结果将通过对该数组的从О到块数量的遍历进行累加计算得到，将在h _c中得到最终的点乘结果。
GPU的计算结果和CPU上通过数学公式计算出来的结果进行校验。在2个输人数组中，如果一个数组中的元素值是从О到N-1，而另外数组中的元素值恒定为2，则它们的点乘结果为N*(N-1)。我们先打印出来GPU的计算结果，然后打印出来它和CPU上的数学公式的比对结果。