CUDA By Examples 5 - 共享内存 Shared Memory

1. 定义

a. 在变量前加上 __shared__, 此变量存储于共享内存中. 1-D或2-D都可以.
b. CUDA在每一个block中都copy了一份此变量.
c. 同一个block中的threads共享此内存, 但是无法读写其他block的拷贝.
d. 共享内存的访问延迟远低于常见的buffer.
e. 需要同步机制(synchronization)来协调threads的读写行为. 防止计算没有完成就执行后续指令.
f. 所有的thread都应该做一样的事. 如果做的事不一样(有if判断), 有些thread可能永远无法执行到同步阶段.

2. 使用共享内存求内积

#include "../common/book.h"

#define imin(a,b) (a<b?a:b)

const int N = 33 * 1024;
const int threadPerBlock = 256;
const int blockPerGrid = imin( 32, (N+threadPerBlock-1) / threadPerBlock );

__global__ void dot( float *a, float *b, float *c)
{
    //共享内存, 每个block都有一份拷贝
    __shared__ float cache[threadPerBlock];
    // thread的索引
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    // 共享内存的索引,每个block都有cache, 故只用threadIdx.x即可
    int cacheIdx = threadIdx.x;

    float temp = 0;
    while(tid<N)
    {
        //当前tid的thread负责把tid,和tid间隔threadIdx总量整数倍的向量做乘-加操作.
        temp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }
    // 完成求和之后,当前thread把和放在对应的cache中
    cache[cacheIdx] = temp;
    // 在当前block内做同步操作, 等所有thread都完成乘-加运算之后才能做reduction.
    __syncthreads();

    //reduction, 向量缩减.
    //缩减后的结果在cache[0]里.
    int i = blockDim.x/2;
    while (i!=0)
    {
        if (cacheIdx<i)
        {
            cache[cacheIdx] += cache[cacheIdx + i];

        }
        //同步, 等所有thread都完成了当次缩减了才能做下一次的缩减.
        //书上说: 同步不能放在if里面, 否则报错.
        //经过试验没有报错, 结果正确.
        __syncthreads();
        i /= 2;
    }
    // 一个block输出一个值,即cache[0]. 所以c的长度和block数量相同.
    // 限制cacheIdx == 0是为了只做一次赋值操作,节省时间.
    if (cacheIdx == 0)
    {
        c[blockIdx.x] = cache[0];
    }
    // 没有做剩下的累加操作是因为在CPU上做小批量的累加更加有效.
}

int main(void)
{
    float *a, *b, c, *partial_c;
    float *dev_a, *dev_b, *dev_partial_c;

    //分配CPU端的内存
    a = (float *)malloc( N*sizeof(float) );
    b = (float *)malloc( N*sizeof(float) );
    partial_c = (float *)malloc( blockPerGrid*sizeof(float));

    //分配GPU端的内存
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a, N*sizeof(float)));
    HANDLE_ERROR( cudaMalloc( (void**)&dev_b, N*sizeof(float)));
    HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c, blockPerGrid*sizeof(float)));

    //将主机内存填入数据
    for (int i=0; i<N; i++)
    {
        a[i] = i;
        b[i] = i*2;
    }

    //将向量a和b拷入GPU
    HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice));
    HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice));

    //GPU上做点积运算
    dot<<<blockPerGrid, threadPerBlock>>>(dev_a, dev_b, dev_partial_c);

    //将向量拷入主机
    HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c, blockPerGrid*sizeof(float), cudaMemcpyDeviceToHost));

    //剩余CPU运算, 求累加和
    c = 0;
    for (int i=0; i<blockPerGrid; i++)
    {
        c += partial_c[i];
    }

    //验证结果是否正确
#define sum_square(x) (x*(x+1)*(2*x+1)/6)
    printf( "Does GPU value %.6g = %.6g?\n",c,
            2 * sum_square( (float)(N-1) ) );
    //释放内存
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_partial_c);

    free( a );
    free( b );
    free( partial_c);

    return 0;
}

3. 不使用同步指令的问题

不使用同步指令, 可能共享内存的计算还没有完成就去执行下面的步骤了. 会造成错误.

#include "cuda.h"
#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define DIM 1024
#define PI 3.1415926535897932f

__global__ void kernel( unsigned char *ptr)
{
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;
    __shared__ float shared[16][16];
    const float period = 128.0f;

    shared[threadIdx.x][threadIdx.y] = 
        255 * (sinf(x*2.0f*PI/period) + 1.0f)*
              (sinf(y*2.0f*PI/period) + 1.0f)/4.0f;
    // 必须加上同步指令.否则有的计算没有完成.
    __syncthreads();
    ptr[offset*4 + 0] = 0;
    ptr[offset*4 + 1] = shared[15-threadIdx.x][15-threadIdx.y];
    ptr[offset*4 + 2] = 0;
    ptr[offset*4 + 3] =255;
}
int main(void)
{
    CPUBitmap bitmap(DIM, DIM);
    unsigned char *dev_bitmap;

    HANDLE_ERROR( cudaMalloc( (void**)&dev_bitmap,
                                bitmap.image_size() ) );
    dim3 grids(DIM/16, DIM/16);
    dim3 threads(16, 16);

    kernel<<<grids, threads>>>( dev_bitmap );

    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
                                bitmap.image_size(),
                                cudaMemcpyDeviceToHost ) );
    bitmap.display_and_exit();

    cudaFree( dev_bitmap );

}

有同步指令的输出(部分图):

没有同步指令的输出(部分图, 错误):