共享内存的合理使用

参考书目： CUDA编程基础与实践

作者: 樊哲勇

章节： 第8章 共享内存的合理使用

数组归约

比如有一个数组[1,2,3,4,5,6,7,8],归约的一种算法如下

//1, 2, 3, 4, 5, 6, 7, 8
//6, 8, 10,12
//16,20
//36

我们可以以这样的方法写一个用block_size=4计算N=8的数组和

#include"error.cuh"
#include<stdio.h>

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif

// const int REPEATS = 100;
const int N = 8;
const int M = sizeof(real) * N;
const int BLOCK_SIZE = 4;//这种归约算法必须满足BLOCK_SIZE>=N/2

void timing(real* h, real* d_x, const int method);
__global__ void reduce_sum(real* h);
real reduce(real* d_x, real* h, const int method);

int main(void)
{
    initDevice(0);//打印设备信息
    real* h = (real*)malloc(M);//主机开辟空间
    for(int i = 0; i<N; i++)//主机空间赋值
    {
        h[i] = 1.23;
    }
    real* d_x;
    // real ans = 1.23*N;
    CHECK(cudaMalloc(&d_x,M));//设备开辟空间
    timing(h, d_x, 0);
    free(h);//释放主机空间
    CHECK(cudaFree(d_x));//释放显存空间
    return 0;

}

void timing(real* h, real* d_x, const int method=0)
{
    CHECK(cudaMemcpy(d_x, h, M, cudaMemcpyHostToDevice));//主机数据拷贝到设备
    /*====记录开始时间====*/
    cudaEvent_t start, stop;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop));
    CHECK(cudaEventRecord(start));//将计时开始事件提交到队列
    cudaEventQuery(start);//刷新队列

    /*====开始GPU计算====*/
    real res = reduce(d_x,h,method);

    /*====记录结束时间====*/
    CHECK(cudaEventRecord(stop));//将计时结束事件提交到队列
    CHECK(cudaEventSynchronize(stop));//让主机等待stop结束,类似刷新队列
    float elapsed_time;//初始化运行时间
    CHECK(cudaEventElapsedTime(&elapsed_time, start, stop));//计算运行时间
    printf("Time = %g ms. \n", elapsed_time);

    /*消除计时事件*/
    CHECK(cudaEventDestroy(start));
    CHECK(cudaEventDestroy(stop));
    printf("ans1 is %f, ans2 is %f\n",h[0],h[1]);
}

__global__ void reduce_sum(real* d_x)
{
    const int tid = blockDim.x*blockIdx.x+threadIdx.x;//计算d_x对应的一维索引
    for(int offset=N>>1; offset>0; offset>>=1)
    {
        if(tid<offset)//只计算前一半数字
        {
            d_x[tid] += d_x[tid+offset];
        }
        //这时可能其他线程还没有算完，如果进入下一个循环会导致计算错误,需要同步
        __syncthreads();//同步一个线程块中的线程,必须在一个block中!!!
        //也就是说必须满足block_size>=N/2
    }
}

real reduce(real* d_x, real* h, const int method)
{
    const int grid_size = (N+BLOCK_SIZE-1)/BLOCK_SIZE;//计算grid_size, 向上取整
    switch (method)
    {
        case 0:
            reduce_sum<<<grid_size, BLOCK_SIZE>>>(d_x);
            break;
        default:
            printf("Error: wrong method\n");
            exit(1);
            break;
    }
    CHECK(cudaMemcpy(h, d_x, M, cudaMemcpyDeviceToHost));//设备数据拷贝到主机
    return h[0];
}

全局内存

这种情况仅适合于BLOCK_SIZE>=N/2的情况,且N为2的幂次的情况

所以可以再创建一个数组用来存每个block的和,也就是说把原数组长度缩小BLOCK_SIZE倍,然后再通过其他算法(for循环)计算其和

核函数为

/*由于之前的reduce_sum只能在block_size>=N/2才能正确计算归约和,所以要优化*/
__global__ void reduce_global(real* d_x, real* d_y)
{
    const int tid = threadIdx.x;//由于每次只在一个block中算,所以id为threadIdx.x
    real* x = d_x+blockDim.x*blockIdx.x;//当前计算的是哪个block中的数,d_x首地址加几个block中的线程数(不需要乘sizeof(real)是因为编译器会自动处理,根据数组元素的大小进行指针移动)
    for(int offset=blockDim.x>>1; offset>0; offset>>=1)//只计算一个block里的归约
    {
        if(tid<offset)//只计算前一半数字
        {
            x[tid] += x[tid+offset];
        }
        //这时可能其他线程还没有算完，如果进入下一个循环会导致计算错误,需要同步
        __syncthreads();//同步一个线程块中的线程,必须在一个block中!!!
        //也就是说必须满足block_size>=N/2
    }
    if(tid==0)
    {
        d_y[blockIdx.x] = x[0];//这个block计算完了就赋值给d_y[blockIdx.x]
    }
}

这样得到了长度为grid_size的数组d_y,d_y[i]表示第i个block的归约和,之后对d_y求和即可

共享内存

以上算法对全局内存频繁访问,由于全局内存的访问延迟最高,所以要减小对其的使用,所以我们可以使用共享内存进行改写

共享内存是对整个线程块可见的区域,对于每个线程都有一个副本(也就是说不同线程块之间这个共享内存所存的数可以不一样),生命周期在核函数内,所以在核函数结束之前必须把其保存到全局内存中

/*使用共享内存*/
__global__ void reduce_shared_memory(real* d_x, real* d_y)
{
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    const int n = blockDim.x*bid+tid;
    __shared__ real s_y[128];//初始化共享内存,长度为BLOCK_SIZE
    s_y[tid] = (n<N)? d_x[n]: 0.0;//给共享内存赋值
    __syncthreads();//保证线程块中所有数据已经准备好
    for(int offset = blockDim.x>>1; offset>0; offset >>= 1)
    {
        if(tid<offset)
        {
            s_y[tid] += s_y[tid+offset];
        }
        __syncthreads();
    }
    if(tid==0)
    {
        d_y[bid]=s_y[0];
    }
}

根据计时结果可以发现,加速效果并不明显,可能是我们的例子共享内存的访问次数还不够多,只有在大数量的访问时才能体现出加速,使用共享内存的还有一个作用可以改善全局内存的访问方式(将非合并改为合并访问)

动态共享内存

上面的例子我们在初始化静态内存时指定了它的长度__shared__ real s_y[128];,它的长度为BLOCK_SIZE,因为我们每次计算的归约都是在一个block里的,所以最长就是BLOCK_SIZE,但如果我们写错这里的长度,就会影响核函数性能,所以我们需要使用动态共享内存,以提高程序的维护性

//配置核函数时指定第三个参数,即共享内存的长度,默认为0
<<<grid_size, block_size, sizeof(real)*block_size>>>

同时声明方式改为

//加限定词extern 并且不指定其大小
extern __shared__ real s[];
//不能声明为 extern __shared__ real *s;

性能和共享内存几乎无差别,但维护性更高

矩阵转置

基本流程为将大矩阵分成若片大小为TILE_SIZE*TILE_SIZE的小矩阵,然后利用并行计算的思想对其进行转置,也就是行列互换

全局内存合并与非合并访问

__global__ void transpose_global1(real* d_x, real* d_y, const int N)
{
    //这个不好解释,建议画图理解
    int nx = blockDim.x*blockIdx.x + threadIdx.x;//大矩阵的x轴索引
    int ny = blockDim.y*blockIdx.y + threadIdx.y;//大矩阵的y轴索引
    if(nx<N&&ny<N)
    {
        //转换为一维数组后的索引
        d_y[ny*N+nx]=d_x[nx*N+ny];//非合并读取,合并写入
        //等价于
        d_y[ny*N+nx]=__ldg(&d_x[nx*N+ny]);
    }
}

可以注意到上面例子在从d_x中读数据时,根据最内层threadIdx.x,可以看出其是非合并读取的,也就是相邻线程(相邻threadIdx.x)所访问的全局内存d_x[nx*N+ny]不是连续的,但是写入的时候是合并的,由于编译器的优化,在读取的时候,会自动用函数__ldg()对数据的读取进行缓存,缓解非合并读取带来的影响

__global__ void transpose_global2(real* d_x, real* d_y, const int N)
{
    int nx = blockDim.x*blockIdx.x + threadIdx.x;
    int ny = blockDim.y*blockIdx.y + threadIdx.y;
    if(nx<N&&ny<N)
    {
        d_y[nx*N+ny]=d_x[ny*N+nx];//合并读取,非合并写入
    }
}

这个例子相反,是合并读取,非合并写入

共享内存改为合并访问

通过共享内存来将非合并改为合并,也就是将共享内存作为中介,从而避开了对全局内存的非合并访问.

__global__ void transpose_shared(real* d_x, real* d_y, const int N)
{
    __shared__ real s_x[TILE_DIM][TILE_DIM];
    int bx = TILE_DIM*blockIdx.x;
    int by = TILE_DIM*blockIdx.y;
    int nx1 = threadIdx.x + bx;
    int ny1 = threadIdx.y + by;
    if(nx1<N&&ny1<N)
    {
        s_x[threadIdx.y][threadIdx.x] = d_x[ny1*N+nx1];//合并读取
    }
    __syncthreads();
    // if(nx1<N&&ny1<N)
    // {
    //     d_y[nx1*N+ny1]=s_x[threadIdx.y][threadIdx.x];// 非合并写入
    //     // d_y[ny1*N+nx1]=s_x[threadIdx.x][threadIdx.y];//这样只能在每个block里转置
    // }

    //合并写入
    //d_y[nx1*N+ny1]=s_x[threadIdx.y][threadIdx.x];// 将这个threadIdx.y与threadIdx.x互换即可
    int nx2 = threadIdx.y + bx;
    int ny2 = threadIdx.x + by;
    if(nx2<N&&ny2<N)
    {
        d_y[nx2*N+ny2]=s_x[threadIdx.x][threadIdx.y];
    }
    __syncthreads();
}

注意不能简单地用d_y[ny1*N+nx1]=s_x[threadIdx.x][threadIdx.y];,这样会使每个block分别转置,与答案不符

消除bank冲突

共享内存在物理上被分为32个相同宽度,又能被同时访问的bank, 可以对其进行0~31编号,每个bank又可以对其分为不同层,开普勒架构之后,每个bank宽度为4字节,所以第一层的32个bank占用字节为0~32*4也就是0~127字节, 第二层开始地址128字节处

何为bank冲突?

只要同一个线程束(相邻32个线程)内多个线程不同时访问同一个bank中的不同层的数据,该线程束对共享内存的方位就只需要一次内存事务,换言之,如果一个线程束中多个线程试图访问一个bank的不同层数据,就会发生bank冲突

在我们之前的例子中,声明了一个TILE_SIZE*TILE_SIZE的单精度浮点型共享内存数组(32*32),我们在转置时d_y[nx2*N+ny2]=s_x[threadIdx.x][threadIdx.y]可以看出在一个线程束内的线程(相邻threadIdx.x)会访问s_x[threadIdx.x][threadIdx.y]也就是相同bank的不同层,所以会引发bank冲突

可以通过__shared__ real s_x[TILE_DIM][TILE_DIM+1];消除bank冲突

__global__ void transpose_shared_bank(real* d_x, real* d_y, const int N)
{
    __shared__ real s_x[TILE_DIM][TILE_DIM+1];//在TILE_DIM=32时, 消除bank冲突
    int bx = TILE_DIM*blockIdx.x;
    int by = TILE_DIM*blockIdx.y;
    int nx1 = threadIdx.x + bx;
    int ny1 = threadIdx.y + by;
    if(nx1<N&&ny1<N)
    {
        s_x[threadIdx.y][threadIdx.x] = d_x[ny1*N+nx1];//合并读取
    }
    __syncthreads();
    int nx2 = threadIdx.y + bx;
    int ny2 = threadIdx.x + by;
    if(nx2<N&&ny2<N)
    {
        d_y[nx2*N+ny2]=s_x[threadIdx.x][threadIdx.y];//这里连续的32个threadIdx.x对应了跨度为TILE_DIM+1的数据，所以不会引起bank冲突
    }
    __syncthreads();
}