最大化GPU内存吞吐量

GPU内存的利用方式对GPU的性能影响很大，因此需要最大程度减少对内存的访问，GPU的内存模型如图1所示。

图1 GPU内存模型示意图

其中最重要的就是尽可能合并全局内存访问的时间，因为不同的内存性能差距比较明显，参考如下优化思路：

尽量使用page-locked(pinned)内存。
相当于单独在CPU上开辟了一块固定内存，当GPU需要数据时，直接从锁页内存中复制，可使用cudaHostAlloc、cudaHostRegister等函数。
尽可能的减少host和device之间的内存传输。
- 将多次传输合并为一次大传输。
- 可以将一些并行度不高但是需要H-D传输的任务放到GPU上。
访问全局内存时，要考虑对齐问题。
如果没有正确的对齐，读写将被编译器拆分为多次操作，降低访存性能。

一维数据使用cudaMalloc()开辟GPU全局内存空间，多维数据建议使用cudaMallocPitch()建立内存空间，以保证段对齐。cudaMallocPitch函数分配的内存中，数组的每一行的第一个元素的开始地址都保证是对齐的。因为每行有多少个数据是不确定的，cudaMallocPitch在分配内存时，每行会多分配一些字节，以保证x方向元素个数*sizeof（元素）+多分配的字节是256的倍数（对齐）。
不要把kernel中的变量放到local Memory中。
尽可能减少每个线程的内存操作次数，注意这个方法不一定有用，因为在有的情况下需要将数据快速复制到shared Memory，然后复制回去，这样虽然总内存操作增多了，但是访存时间下降了。
比如求解数组平方和的案例，当一个thread在等待内存数据的时候，GPU就会切换到下一个thread。所以，实际执行的顺序类似于thread0 --> thread1 --> ... ... --> threadn，这就导致了同一个thread在读取内存是连续的，但是对于整体而言，执行的过程中内存读取就不是连续的。
```
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
    const int size = DATA_SIZE / THREAD_NUM;
    const int tid = threadIdx.x;
    int tmp_sum = 0;
    clock_t start;
    if (tid == 0) start = clock();
    for (int i = tid * size; i < (tid + 1) * size; i++) 
    {
        tmp_sum += data[i] * data[i];
    }
    sum[tid] = tmp_sum;
    if (tid == 0) *time = clock() - start;
}
```
如果在读取数组data时，能够做到每隔N个线程读取一次，所有的内存操作都会统一，GPU线程就不会出现切换，大大减少了开销。下面的代码实测比上面的代码性能提升15倍。
```
__global__ static void squaresSum(int *data, int *sum, clock_t *time)
{
 const int size = DATA_SIZE / THREAD_NUM;
    const int tid = threadIdx.x;
    int tmp_sum = 0;
    clock_t start;
    if (tid == 0) start = clock();
    for (int i = tid; i < DATA_SIZE; i += THREAD_NUM)
    {
      tmp_sum += data[i] * data[i];
    }
    sum[tid] = tmp_sum;
    if (tid == 0) *time = clock() - start;
}
```
避免线程发散（thread divergence），GPU相邻线程执行相同的操作效率会更高。
比如常用的reduce操作，假设给定一个长度为N的数组，需要计算该数组的所有元素之和，算法比较简单，分为三个步骤。
1. 将数据load至shared memory中。
2. 在shared memory中对数据进行reduce操作。
3. 将最后的结果写回global memory中。
```
__global__ void reduce(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    //each thread loads one element from global memory to shared mem
    unsigned int i=blockIdx.x*blockDim.x+threadIdx.x;
    unsigned int tid=threadIdx.x;
    sdata[tid]=d_in[i];
    __syncthreads();

    // do reduction in shared mem
    for(unsigned int s=1; s<blockDim.x; s*=2){
        if(tid%(2*s) == 0){
            sdata[tid]+=sdata[tid+s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if(tid==0)d_out[blockIdx.x]=sdata[tid];
}
```
以每个线程块中包含256个线程为例(THREAD_PER_BLOCK=256)，首先256个线程分为8组warp，每组32个线程，第tid号线程将i号数据从global memory取出放到shared memory中，接下来开始迭代操作：
1. 第1轮迭代：如果tid%2 ==0，则第tid号线程将shared memory中第tid号位置的值和第tid+1号的值进行相加，而后放在第tid号位置。
2. 第2轮迭代：如果tid%4==0，则第tid号线程将shared memory中第tid号位置的值和第tid+2号的值进行相加，而后放在第tid号位置。
3. 不断迭代，直到所有元素都被累加到第0号位置。
目前上述代码存在的最大问题就是warp divergent的问题。对于一个block而言，它所有的thread都是执行同一条指令。如果存在if-else这样的分支情况的话，thread会执行所有的分支。只是不满足条件的分支，所产生的结果不会记录下来。可以在上图中看到，在每一轮迭代中都会产生两个分支，分别是红色和橙色的分支。这严重影响了代码执行的效率。如果能尽可能地让所有线程走到同一个分支里面效率应当会提高。

优化后：
```
__global__ void reduce1(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    //each thread loads one element from global memory to shared mem
    unsigned int i=blockIdx.x*blockDim.x+threadIdx.x;
    unsigned int tid=threadIdx.x;
    sdata[tid]=d_in[i];
    __syncthreads();

    // do reduction in shared mem
    for(unsigned int s=1; s<blockDim.x; s*=2){
        int index = 2*s*tid;
        if(index < blockDim.x){
            sdata[index]+=sdata[index+s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if(tid==0)d_out[blockIdx.x]=sdata[tid];
}
```
虽然代码依旧存在着if语句，但是却与之前的代码有所不同。继续假定block中存在256个thread，即拥有256/32=8个warp。

当进行第1次迭代时，0~3号warp的index<blockDim.x，4~7号warp的index>=blockDim.x。对于每个warp而言，都只是进入到一个分支内，所以并不会存在warp divergence的情况。

当进行第2次迭代时，0、1号两个warp进入计算分支。

当进行第3次迭代时，只有0号warp进入计算分支。

当进行第4次迭代时，只有0号warp的前16个线程进入分支。此时开始产生warp divergence。通过这种方式，消除了前3次迭代的warp divergence。
使用shared Memory时，防止bank conflict。
reduce1的最大问题是bank冲突。我们把目光聚焦在这个for循环中。并且只聚焦在0号warp。

在第一次迭代中，0号线程需要去load shared memory的0号地址以及1号地址的数，然后写回到0号地址。而此时，这个warp中的16号线程，需要去load shared memory中的32号地址和33号地址。可以发现，0号地址跟32号地址产生了2路的bank冲突。

在第2次迭代中，0号线程需要去load shared memory中的0号地址和2号地址。这个warp中的8号线程需要load shared memory中的32号地址以及34号地址，16号线程需要load shared memory中的64号地址和68号地址，24号线程需要load shared memory中的96号地址和100号地址。又因为0、32、64、96号地址对应着同一个bank，所以此时产生了4路的bank冲突。

现在，可以继续算下去，8路bank冲突，16路bank冲突。由于bank冲突，所以reduce1性能受限。下图说明了在load第一个数据时所产生的bank冲突。

解决bank冲突的方式就是把for循环逆着来。原来stride从0到256，现在stride从128到0。
```
__global__ void reduce2(float *d_in,float *d_out){
    __shared__ float sdata[THREAD_PER_BLOCK];

    //each thread loads one element from global memory to shared mem
    unsigned int i=blockIdx.x*blockDim.x+threadIdx.x;
    unsigned int tid=threadIdx.x;
    sdata[tid]=d_in[i];
    __syncthreads();

    // do reduction in shared mem
    for(unsigned int s=blockDim.x/2; s>0; s>>=1){
        if(tid < s){
            sdata[tid]+=sdata[tid+s];
        }
        __syncthreads();
    }

    // write result for this block to global mem
    if(tid==0)d_out[blockIdx.x]=sdata[tid];
}
```
继续分析如何通过一个小小的改变就能消除bank冲突。

首先，注意这个for循环，并且只分析0号warp。0号线程需要load shared memory的0号元素以及128号元素。1号线程需要load shared memory中的1号元素和129号元素。这一轮迭代中，在读取第一个数时，warp中的32个线程刚好load一行shared memory数据。

再分析第2轮迭代，0号线程load 0号元素和64号元素，1号线程load 1号元素和65号元素，每个线程load shared memory的一行。

再分析第3轮迭代，0号线程load 0号元素和32号元素，以此类推，一个线程load shared memory的一行，没有bank冲突。

到了4轮迭代，0号线程load 0号元素和16号元素。

由于s=16，16~31号被跳过。迭代示意图如下。

更多reduce函数优化请参考：https://github.com/Liu-xiandong/How_to_optimize_in_GPU

尽可能连续的访问全局内存。

比如启动N*N个线程计算两个N*N矩阵的乘积，结果矩阵的每个值使用下述公式计算：

$\text{[math]}$

在计算C_ij过程中，B_kj涉及非连续数据的访问，影响访存效率，此时可以引入矩阵B^T，使得数据连续访问：

$\text{[math]}$

优化前：每个线程只计算1个元素，B数组非连续访问。

__global__ void matrix_mul_gpu(int n, float *a, float*b, float*c)
{
   const int bid = blockIdx.x;
   const int tid = threadIdx.x;
      float s = 0.0;
        for (int k = 0; k < n; ++k)
           s += a[bid*n + k] * b[k*n + tid];               //A(i,k)*B(k,j)
         c[bid*n+tid] = s;
}

优化后：A，B矩阵均连续访问，相比优化前实现了30%性能提升。

__global__ void transpose(int n, float *b, float *bt)
{
bt[blockIdx.x*n + threadIdx.x] = b[threadIdx.x*n + blockIdx.x];
}
__global__ void matrix_mul_gpu_opt(int n, float *a, float*b, float*c)
{
const int bid = blockIdx.x;
const int tid = threadIdx.x;
float s = 0.0;
for (int k = 0; k < n; ++k)
  s+= a[bid*n + k] * b[tid*n + k];          //A*BT  连续访问
 c[bid*n + tid] = s;
}

父主题： kernel代码优化