本文只是一个CUDA编程小白记录学习过程的。更全面的理解和学习请跳转链接

一些前置知识的参考链接：
线程束分化
bank conflict
本文主要学习以下手段对cuda代码进行优化：

• 解决warp divergence
• 解决bank conflict
• 解决idle线程

reduce实现向量元素累加

下面的程序实现一个kernel函数，通过reduce方法实现向量的累加和。代码采用多block实现的，每次调用kernel函数只计算blocksize个元素的累加和。

#include <stdio.h>
#include <cuda_runtime.h>
#include "nvToolsExt.h"
#define THREAD_PER_BLOCK 256
#define N 2048


__host__ void initialData(float *h_a, int n)
{
    for(int i = 0; i < n; i++)
    {
        h_a[i] = 1.0;
	    //printf("h_a=%f\n", h_a[i]);
    }
}

__global__ void reduce0(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];
}

int main()
{
    cudaError_t res;
    int BLOCKS =  (N + THREAD_PER_BLOCK -1) / THREAD_PER_BLOCK;

    float *h_a = (float *)malloc(N*sizeof(float));
    float *h_out = (float*)malloc(BLOCKS*sizeof(float));
    initialData(h_a, N);


    float *d_a=NULL, *d_out = NULL;
  

     res =  cudaMalloc((float**)&d_a, N*sizeof(float));
     if(res != cudaSuccess ) printf("test 1\n");
     res = cudaMalloc((float**)&d_out, BLOCKS*sizeof(float));
     if(res != cudaSuccess ) printf("test 3\n");
     
     res = cudaMemcpy(d_a, h_a, N*sizeof(float), cudaMemcpyHostToDevice);   
     if(res != cudaSuccess ) printf("test 2\n");

     //调用kernel函数
     dim3 blocks(BLOCKS, 1, 1);
     dim3 threadPerBlock(THREAD_PER_BLOCK, 1, 1);
     printf("BLOCKS=%d\n", BLOCKS);
     printf("THREAD_PER_BLOCK=%d\n", THREAD_PER_BLOCK);
     nvtxRangePushA("reduce0 cuda kernel func");
     reduce0<<<blocks, threadPerBlock>>>(d_a, d_out);
     nvtxRangePop();
     cudaError_t err = cudaGetLastError();
     if(err != cudaSuccess ) printf("Error: %s\n", cudaGetErrorString(err));
     
     res = cudaMemcpyAsync(h_out, d_out, BLOCKS * sizeof(float), cudaMemcpyDeviceToHost);
     if(res != cudaSuccess ) printf("test 5\n");

     for(int i = 0; i < BLOCKS; ++i)
     {
	    printf("block sum = %f\n", h_out[i]);
     }

     //分配的显存需要手动释放
     free(h_a);
     free(h_out);
     cudaFree(d_a);
     cudaFree(d_out);
     
    return 0;
}

编译

nvcc -o main -G main.cu -lnvToolsExt

使用nsys测试性能：

nsys profile --stats=true ./main
调用nvtx接口，测试关键代码的运行时间。

	nvtxRangePushA("reduce0 cuda kernel func");
     reduce0<<<blocks, threadPerBlock>>>(d_a, d_out);
     nvtxRangePop();

使用nsys-ui查看

解决warp divergence

通过与最基础的代码做对比。一个block中有256个线程，对应256/32=4个warp。在reduce0中通过判断线程id的奇偶做reduce，这势必会带来线程束分化。
因此在代码中引入index，index = 2* s* tid; 通过判断index< blockDim.x ， 在第一轮reduce时将每次进入if语句的线程id控制在0 ~ 127号。第二轮reduce进入if的线程id为0 ~ 63号。因此避免了线程束内部线程分化。

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

解决bank conflict

THREAD_PER_BLOCK=256，在当前的一个block中，shared memory存储的元素个数256。因此在同时访问下标为0， 32， 64 的数据时会出现bank conflict。因此解决bank conflict的方法是for循环逆着来。这样在第一轮reduce分别计算累加的元素下标为 0+128、 1+129、 2+130 …32 + 160（已经在第2个warp中了） 。第二轮reduce计算元素和的下标为：0 + 64、 1 + 65 、 … 、 63 + 127。

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