3、向量点乘和矩阵乘法的例子
3.1、向量点乘
两个向量的点乘是重要的数学运算，也将会解释CUDA编程中的一个重要概念：归约运算。两个向量的点乘运算定义如下：
其实显示的应用中真正的向量肯定会很长很长，两个向量里面有多个元素，而不仅仅只有三个。最终也会将多个乘法结果累加（归约运算）起来，而不仅仅是3个。现在，你看下这个运算，它和之前的元素两两相加的向量加法操作很类似。不同的是你需要将元素两两相乘。线程需要将它们的所有单个乘法结果连续累加起来，因为所有的一对对的乘法结果需要被累加起来，才能得到点乘的最终结果。最终的点乘的结果将是一个单一值。这种原始输入是两个数组而输出却缩减为一个（单一值）的运算，在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) { //声明共享内存 __shared__ float partial_sum[threadsPerBlock]; int tid = threadIdx.x + blockIdx.x * blockDim.x; //计算共享内存的index int index = threadIdx.x; //计算部分和 float sum = 0; while (tid < N) { sum += d_a[tid] * d_b[tid]; tid += blockDim.x * gridDim.x; } // 存储部分和结果到共享内存中 partial_sum[index] = sum; // 同步线程 __syncthreads(); //为整个Block计算部分和 int i = blockDim.x / 2; while (i != 0) { if (index < i) partial_sum[index] += partial_sum[index + i]; __syncthreads(); i /= 2; } //存储部分和到全局内存 if (index == 0) d_c[blockIdx.x] = partial_sum[0]; } int main(void) { //声明Host数组 float *h_a, *h_b, h_c, *partial_sum; //声明device数组 float *d_a, *d_b, *d_partial_sum; //计算每个grid全部block的数目 int block_calc = (N + threadsPerBlock - 1) / threadsPerBlock; int blocksPerGrid = (32 < block_calc ? 32 : block_calc); // 申请主机内存 h_a = (float*)malloc(N * sizeof(float)); h_b = (float*)malloc(N * sizeof(float)); partial_sum = (float*)malloc(blocksPerGrid * sizeof(float)); // 申请显存 cudaMalloc((void**)&d_a, N * sizeof(float)); cudaMalloc((void**)&d_b, N * sizeof(float)); cudaMalloc((void**)&d_partial_sum, blocksPerGrid * sizeof(float)); // 喂主机数组数据 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); //调用Kernel函数 gpu_dot << <blocksPerGrid, threadsPerBlock >> > (d_a, d_b, d_partial_sum); cudaMemcpy(partial_sum, d_partial_sum, blocksPerGrid * sizeof(float), cudaMemcpyDeviceToHost); 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"); } 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等于这些部分和的个数（又等于块中的线程数量）。另一种就是使用并行化的方法，而且复杂度是log2(N)比我们第一中方法的单线程N的复杂度好很多。
这种并行归约运算是通过条件为（i!=0）的while循环进行的。当前块中每个满足一定条件的线程每次循环累加当前它的ID索引的部分和，和ID+blockDim/2作为索引的部分和，并将结果回写共享内存。重复这个过程直到得到最终的单一结果，即整个块最终的一个部分和。最终每个块的单一部分和结果被写入到共享内存中（的特定位置）。通过块的唯一ID进行索引，每个块能确定这个单独的写入特定位置。然而，我们还是没有得到最终结果。这个最后的结果通过设备上的函数，或者（主机上）的main函数执行都可以。
一般情况下，最后几次归约累加只需要很少的资源。如果我们在GPU上进行的话，大部分GPU的资源都将空闲，无法有效利用GPU。但是根据CUDA手册的建议，即使在并行度不佳的情况下，也应当尽量开展在GPU上计算，本例子最后可以使用atomicAdd来完成最终计算，根据实际经验，这种效果往往较好。因为随着显存读写速度和PCI-E传输带宽的差距进一步加大，一次传输的成本，可能远比在GPU上，哪怕是并行度不佳的时候的计算成本要高。复制回来并不一定合算。
3.2、矩阵乘法
除了向量点乘之外，GPU上通过CUDA进行的最重要的数学运算是矩阵乘法。当矩阵非常大的时候，数学运算将非常复杂。需要记住的是，当进行矩阵乘法的时候，乘号前矩阵的列数（即行宽）需要等于乘号后矩阵的行数（即高度）。
矩阵乘法不满足交换律。示例矩阵运算如下：
矩阵乘法中同样的数据能被使用多次，所以这是共享内存的一个理想用例。
这里的blockDim.x和blockDim.y等于TILE_SIZE的大小。现在，结果矩阵中的每个元素都将是第一个矩阵的（对应）行和第二个矩阵的（对应）列进行向量点乘的结果。因为两个矩阵都是变量size*size这么大的方阵，所以需要是size个元素这么大的两个向量进行点乘。所以（每个线程）通过内核函数代码里的for循环从0循环到size，进行了多次乘加。
为了能够在两个矩阵中一个一个的索引元素（计算元素的索引），可以将矩阵看成是在存储器中按照行主序的方式线性存储的。行主序的意思是，第一行的元素一个接一个地连续在存储器中排列，然后下一行再继续（这样排列）在前一行的后面。如下所示：
每个元素的线性索引可以这样计算：用它的行号乘以矩阵的宽度，再加上它的列号即可。所以，对于M1,0来说，它的线性索引值是2。因为它的行号是1，矩阵的宽度是2，它的列号是0。这种方法用来计算两个源矩阵中元素的索引。
为了计算结果矩阵中在[row,col]位置的元素的值，需要前一个矩阵中的一行和后一个矩阵的一列进行点乘。内核里分别用了row*size+k的下标和k*size+col的下标来索引这一行和一列中的元素。
要计算最终的目标矩阵中的[row,col]的元素，需要将第一个矩阵中的第row行和第二个矩阵中的第col列进行向量点乘。而这向量点乘则需要循环多组对应的点，分别对来自第一个矩阵从2D索引线性化到1D索引后的坐标为row*size+k的点，和来自第二个矩阵的同样1D化坐标为k*size+col处的点进行乘法累加操作，其中k是每次循环时候的变量。这样就完成了向量点乘，以及最后的矩阵乘法。
用共享内存中的大小为TILE_SIZE*TILE_SIZE的两个方块来保存需要重用的数据。就像你之前看到的那样的方式来计算行和列索引。首次for循环的时候，先填充共享内存，然后调用__syncthreads()，从而确保后续从共享内存的读取只会发生在块中的所有线程都完成填充共享内存之后。然后循环的后面部分则是计算（部分）向量点乘。因为（大量的重复读取）都只发生在共享内存上，从而显著地降低对全局内存的访存流量，进一步地能提升大矩阵维度下的乘法程序的性能。
当定义了主机和设备指针变量并为它们分配空间后，主机上的数组被填充了一些任意的数字（注意不是随机，因为是写死的），然后这些数组被复制到显存上，这样它们才能作为参数传输给内核函数。接着使用dim3结构体定义Grid中的块和块中的线程形状，具体的形状信息是之前几行算好的，你可以调用这两个内核中的任何一个（使用共享内存的和不使用共享内存的），然后将结果复制回主机来。
代码如下：
//Matrix multiplication using shared and non shared kernal #include "stdio.h" #include<iostream> #include <cuda.h> #include <cuda_runtime.h> #include <math.h> #define TILE_SIZE 2 //不使用共享内存运算矩阵相乘的Kernel函数 __global__ void gpu_Matrix_Mul_nonshared(float *d_a, float *d_b, float *d_c, const int size) { int row, col; col = TILE_SIZE * blockIdx.x + threadIdx.x; row = TILE_SIZE * blockIdx.y + threadIdx.y; for (int k = 0; k < size; k++) { d_c[row*size + col] += d_a[row * size + k] * d_b[k * size + col]; } } // 使用共享内存运算矩阵相乘的Kernel函数 __global__ void gpu_Matrix_Mul_shared(float *d_a, float *d_b, float *d_c, const int size) { int row, col; //Defining Shared Memory __shared__ float shared_a[TILE_SIZE][TILE_SIZE]; __shared__ float shared_b[TILE_SIZE][TILE_SIZE]; col = TILE_SIZE * blockIdx.x + threadIdx.x; row = TILE_SIZE * blockIdx.y + threadIdx.y; for (int i = 0; i < size / TILE_SIZE; i++) { shared_a[threadIdx.y][threadIdx.x] = d_a[row* size + (i*TILE_SIZE + threadIdx.x)]; shared_b[threadIdx.y][threadIdx.x] = d_b[(i*TILE_SIZE + threadIdx.y) * size + col]; __syncthreads(); for (int j = 0; j < TILE_SIZE; j++) d_c[row*size + col] += shared_a[threadIdx.y][j] * shared_b[j][threadIdx.x]; __syncthreads(); } } int main() { const int size = 4; float h_a[size][size], h_b[size][size], h_result[size][size]; float *d_a, *d_b, *d_result; for (int i = 0; i < size; i++) { for (int j = 0; j < size; j++) { h_a[i][j] = i; h_b[i][j] = j; } } cudaMalloc((void **)&d_a, size*size * sizeof(int)); cudaMalloc((void **)&d_b, size*size * sizeof(int)); cudaMalloc((void **)&d_result, size*size * sizeof(int)); cudaMemcpy(d_a, h_a, size*size * sizeof(int), cudaMemcpyHostToDevice); cudaMemcpy(d_b, h_b, size*size * sizeof(int), cudaMemcpyHostToDevice); //定义grid和block的维度 dim3 dimGrid(size / TILE_SIZE, size / TILE_SIZE, 1); dim3 dimBlock(TILE_SIZE, TILE_SIZE, 1); //gpu_Matrix_Mul_shared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size); gpu_Matrix_Mul_nonshared << <dimGrid, dimBlock >> > (d_a, d_b, d_result, size); cudaMemcpy(h_result, d_result, size*size * sizeof(int), cudaMemcpyDeviceToHost); printf("The result of Matrix multiplication is: \n"); for (int i = 0; i < size; i++) { for (int j = 0; j < size; j++) { printf("%f ", h_result[i][j]); } printf("\n"); } cudaFree(d_a); cudaFree(d_b); cudaFree(d_result); return 0; }
