CUDA_编程：基础与实践的阅读笔记和program理解 Chapter_07-usage_of_the_global_memory

CUDA 编程：基础与实践的阅读笔记和program理解 Chapter 07-usage of the global memory

《CUDA 编程：基础与实践》(樊哲勇)【简介_书评_在线阅读】 - 当当图书 (dangdang.com)

ZouJiu1/CUDA-Programming: Sample codes for my CUDA programming book (github.com)

全局内存的使用

顺序的合并访问 ，x, y 和 z 是由cudaMalloc()分配全局内存的指针，首地址至少是256的整数倍，核函数对这几个指针所指内存区域的访问都是合并的。grid size = gridDim.x = 64x2，block size = blockDim.x = 32。所以第一个线程块block内的线程束(32个线程)会访问数组x的第0-31个数值，sizeof(float) = 4，对应32x4个字节的连续内存，而且首地址至少是256的整数倍，也就是首地址是最小粒度的整数倍，访问只需要4次数据传输即可，是合并访问，合并度=100%。

void  __global__  add(float *x, float *y, float *z)
{
    int n = threadIdx.x + blockIdx.x * blockDim.x;
    z[n] = x[n] + y[n];
}
add<<<128, 32>>>(x, y, z);

随机的合并访问 ，交换两个相邻的数字，第一个线程块内的线程束仍然是访问的x内的第0-31个数据，不过线程号和数组的index不完全一致。访问是交叉的合并访问，合并度=100%，此时线程束对内存的访问还是合并的

void __global__ add_permuted(float *x, float *y, float *z)
{
    int tid_permuted = threadIdx.x ^ 0x1;
    int n = tid_permuted + blockIdx.x * blockDim.x;
    z[n] = x[n] + y[n];
}
add_permuted<<<128, 32>>>(x, y, z);

不对齐的非合并访问 ，第一个线程块内的线程束会访问数组x的第1-32个数值，假若数组x的首地址是256字节，该线程束会访问设备内存的260-(260+127)字节，会触发5次数据传输，数据传输对数据地址的要求， 在一次数据传输内，从全局内存转移到L2缓存的一片内存的首地址一定是一个最小粒度(上述是32字节)的整数倍。 所以，对应的内存地址分别是256_{(256+31)，(256+32)}(256+32x2-1)，…，属于不对齐的非合并访问，合并度则是：16 / 20

上述最后还是使用的设备内存的260~(260+127)字节内的数值，但是传输要求从最小粒度的整数倍开始传输，所以要从256开始传输，也就是多传输了很多个字节，这些多的字节是没用的。浪费了显存带宽。

int n = threadIdx.x + blockIdx.x * blockDim.x +1; 会尝试读取index所在的后一个数值

void __global__ add_offset(float *x, float *y, float *z)
{
    int n = threadIdx.x + blockIdx.x * blockDim.x + 1;
    z[n] = x[n] + y[n];
}
add_offset<<<128, 32>>>(x, y, z);

跨越式的非合并访问， 第一个线程块内的线程束会访问数组x内指标index是0、64x2、64x4、64x6的数。第一个线程块的blockIdx.x = 0，第一个线程块内的线程束的threadIdx.x = 0，1，2，…，对应的网格gridDim.x = 64x2，所以需要访问跨越的内存地址才能拿到数值，像0, 64x2, 64x4, 64x6等。线程束共有32个线程，所以需要传输数据32次，每次传输32个字节但是只需要用到4个字节，所以属于跨越式的非合并访问，合并度则是4/32 = 12.5%

void __global__ add_stride(float *x, float *y, float *z)
{
    int n = blockIdx.x + threadIdx.x * gridDim.x;
    z[n] = x[n] + y[n];
}
add_stride<<<128, 32>>>(x, y, z);

广播式的非合并访问， 第一个线程块内的线程束会一致的访问数组x内的第0个数值，只需要一次数据传输（处理32字节的数据），但由于整个线程束只使用了4字节的数据，也就是传输了32字节但是只需要用到4字节，故合并度是4/32 = 12.5%。属于广播式的非合并访问，适合使用之前chapter 6提到的常量内存，也就是 constant int a = x[0]，或者直接用传递参数的方式传递过来。会放到寄存器内部。

void __global__ add_broadcast(float *x, float *y, float *z)
{
    int n = threadIdx.x + blockIdx.x * blockDim.x;
    z[n] = x[0] + y[n];
}
add_broadcast<<<128, 32>>>(x, y, z);

矩阵转置

转置也就是 $矩阵A，C=A^T，C_{ij}=(A^T{)}_{ij}=A_{ji}$

${\left[\begin{array}{ccc}1& 2& 3\\ 10& 5& 6\\ 7& 8& 9\end{array}\right]}^T=\left[\begin{array}{ccc}1& 10& 7\\ 2& 5& 8\\ 3& 6& 9\end{array}\right]$

复制

对矩阵做了两次分片，第一个是粗粒度的分片，像一个900 x 900的矩阵，可以先分片，每个分片大小是30 x 30，也就是行分到30份，列也分到30份，也就是(30x30) x (30x30)；然后每个小网格再次分片，也就是30x30的网格再次分片，此时就可以直接行分30，列分30，就是(30x1)x(30x1)最后的分片结果是(30x(30x1))x(30x(30x1))；还可以行分3份，列分3份就是(10x3)x(10x3)，最后的分片结果是(30x(10x3))x(30x(10x3))。

上面的图片就是对矩阵做了两次分片，第一次粗粒度分片是黑色实线，第二次分片是虚线，所以上面的网格分片的结果是黑色线3x3，虚线3x3，总体是(3x(3x1))x(3x(3x1))。

{
    const int TILE_DIM = 32;
    const int grid_size_x = (N + TILE_DIM - 1) / TILE_DIM;
    const int grid_size_y = grid_size_x;
    const dim3 block_size(TILE_DIM, TILE_DIM);
    const dim3 grid_size(grid_size_x, grid_size_y);
    copy<<<grid_size, block_size>>>(d_A, d_B, N);
}
__global__   void  copy(const real *A, real *B, const int N)
{
    const int nx = blockIdx.x * TILE_DIM + threadIdx.x;
    const int ny = blockIdx.y * TILE_DIM + threadIdx.y;
    const int index = ny * N + nx;
    if (nx < N && ny < N)
    {
        B[index] = A[index];
    }
}

  上面的code就是复制矩阵到另一个数列，grid_size和block_size都是2dim的，也就是grid size是粗粒度分片，然后block size细粒度分片，和上面的图片表示的内容是类似的，宽和高的方向都有相应的grid size，每个小矩形内都有相应的block size。grid size是每个grid所在方向包含的block块个数，block size是每个block所在方向包含的thread线程的个数。

  上面的index运算方式是： **复合线程索引** ，首先算出横轴和纵轴方向的坐标，TILE_DIM也就是每个方向的block size也就是每个方向的线程数量，blockIdx.x是x方向的block的Index，threadIdx.x是x方向线程在block内部的Index，所以nx拿到的就是x方向的线程坐标，同样ny拿到的就是y方向的线程坐标。N是方阵的长或者宽。此时index = ny * N + nx 拿到的就是线程实际的标号index。因分配的是1dim的内存，所以拿到的就是 1dim 内存的标号index。if(nx < N && ny < N) 判断条件是防止行指标和列指标越界的。

  而且还需要注意的就是：核函数内部可以直接使用在函数外部由 #define 或者 const 定义的常量，包括整型常量和浮点型常量，在MSVC使用时有个限制，不能在核函数内使用在函数外部由const定义的浮点型常量。可以直接使用常量的值，但不能在核函数内使用这种常量的引用或者地址。

    const int nx = blockIdx.x * TILE_DIM + threadIdx.x;
    const int ny = blockIdx.y * TILE_DIM + threadIdx.y;
    const int index = ny * N + nx;

  核函数内对全局内存的访问模式，x方向的线程指标threadIdx.x是最内层的呢，所以相邻的threadIdx.x对应相邻的线程，从而相邻的nx对应相邻的线程，也对应相邻的数组element，所以核函数内，相邻的线程访问了相邻的数组element，在内存对齐的情况下属于上个 chapter 介绍的顺序的合并访问。合并度=100%。作者令N = 10000并在GeForce RTX2080ti上进行实测。使用单精度浮点数，copy复制核函数执行时间是1.6ms，根据该执行时间，有效的显存带宽是  $2 * 10000 * 10000 * (2 * 2)  \text{B} / 1.6 \text{ms} \\ = 2 * 10000 * 10000 * (2 * 2) \text{B}  / (1.6*0.001) \text{s} \\ =  500000000000 \text{(b/s)}$ 

$500000000000/1024/1024/1024\text{(gb/s)}=465.66\text{(gb/s)}$

465.66 gb / s ~ 500 gb / s ，(2 x 2) B = sizeof(float) 单精度占用的字节数，10000 * 10000 是矩阵数组的长度，2是需要访问2个数组A和B，略小于该GPU的理论显存带宽616 gb / s 。

使用全局内存进行矩阵转置

上述核函数内的主要赋值code是，相当是做了 $B_{ij}=A_{ji}$ 操作的

const int index = ny * N + nx;
if (nx < N && ny < N)  B[index] = A{index];
还可以改写到：
if (nx < N && ny < N) B[ny * N + nx] = A[ny * N + nx];

若要实现转置，只需要修改code就可以，也就是 $B_{ij}=A_{ji}$

const int index = ny * N + nx;
if (nx < N && ny < N) B[nx * N + ny] = A[ny * N + nx];     (1)
或者
if (nx < N && ny < N) B[ny * N + nx] = A[nx * N + ny];      (2)

也就是ny行nx列的值，赋值给nx行，ny列，刚好就是转置了，但是上述两类转置的方式，性能是不相同的，第(1)类方式，对矩阵A的数据访问（读取）是顺序的，但对矩阵B内数据的访问（写入）不是顺序的。在第(2)类方式，则正好相反，对矩阵A的数据访问（读取）不是顺序的，但对矩阵B内数据的访问（写入）是顺序的。若是不考虑数据是否对齐的情况下，就说核函数(1)对矩阵A和B的访问分别是合并的和非合并的，核函数(2)对矩阵A和B的访问分别是非合并的和合并的。

__global__ void transpose1(const real *A, real *B, const int N) //5.3ms
{
    const int nx = blockIdx.x * blockDim.x + threadIdx.x;
    const int ny = blockIdx.y * blockDim.y + threadIdx.y;
    if (nx < N && ny < N)
    {
        B[nx * N + ny] = A[ny * N + nx];
    }
}
__global__ void transpose2(const real *A, real *B, const int N)  //2.8ms
{
    const int nx = blockIdx.x * blockDim.x + threadIdx.x;
    const int ny = blockIdx.y * blockDim.y + threadIdx.y;
    if (nx < N && ny < N)
    {
        B[ny * N + nx] = A[nx * N + ny];
    }
}

用单精度浮点数运算，GeForce RTX 2080ti实测核函数的执行时间。两者的耗时是不相同的，transpose2读取操作是非合并的，但使用了只读数据缓存的加载函数__ldg()。从Pascal架构开始的话，若是编译器能够判断一个全局内存变量在整个核函数内都只可读 像矩阵A，会自动调用函数__ldg()读取全局内存，从而对数据的读取进行缓存，缓解非合并访问带来的影响。对于全局内存的写入，则没有类似的函数可用。这就是以上两个核函数性能差别的reason。所以在不能同时满足读取和写入都是合并的情况下，一般来说应当尽量做到合并的写入。让函数__ldg()来缓解非合并访问的延迟。

Kepler和Maxwell架构，默认不会使用__ldg()函数，需要人工指定才可以的也就是transpose3核函数。

https://zhuanlan.zhihu.com/p/654562406