CUDA入门（二）：数据并行执行模型
内置变量：内置变量的值一般都由运行时系统预初始化，例如 CUDA 的 kernel 函数中， gradDim，blockDim，blockIdx，threadIdx都是内置变量，它们的值由 CUDA 运行时系统预初始化，可以在 kernel 函数中引用。其他地方应避免使用这些变量。
1. CUDA的线程组织
上一节向量加法中的线程被组织成二级的层次结构：一个网络包含一个或更多的线程块，每块包含一个或更多的线程。一个块中所有线程的 blockId 相同，每个块中又可通过唯一的 threadIdx 访问唯一的线程。
网络一般是由线程块组成的三维数组，线程块又是线程组成的三维数组。不需要的维度可以设为 1 。网络确切组织形式是由 kernel 函数启动语句中的配备参数决定的 (<<< 配备参数 1， 配备参数 2>>>)，参数 1 确定了网格的维度和线程块的数量，参数 2 确定了线程块的维度和每个线程块内线程的数量。每个参数都是 dim3 类型，dim3 类型为含有 x, y 和 z 三个无符号整形的结构体，分别对应三个维度。
下面的主机端代码启动 vecAddKernel() 函数，生成了包含 32 个一维线程块的一维网格，每个线程块包含 128 个线程：
dim3 dimGrid(32, 1, 1); dim3 dimBlock(128, 1, 1); vecAddKernel<<<dimGrid, dimBlock>>>(...);
对于向量加法，可以写成如下形式：
dim3 dimGrid(ceil(n/256.0), 1, 1); dim3 dimBlock(256, 1, 1); vecAddKernel<<<dimGrid, dimBlock>>>(...);
CUDA C提供了用一维网络和线程启动 kernel 的简单方式，就跟上一节的一样：
vecAddKernel<<<ceil(n/256.0), 256>>>(...);
在 CUDA C 中 gridDim.x，gridDim.y 和 gridDim.z 的取值范围是 1 ~ 65536。在所有的线程块中，blockIdx.x 的取值范围从 0 到 gridDim.x - 1，blockIdx.y 的取值范围从 0 到 gridDim.y - 1，blockIdx.z 的取值范围从 0 到 gridDim.z - 1。
在线程块的配置上，将 z 字段设为 1 可得到二维线程块，将 y, z 设为 1 可得到一维线程块。计算能力低于 2.0 的设备一个线程块最多包含 512 个线程，大于等于 2.0 的设备一个线程块最多包含 1024 个线程。（512， 1， 1）、（8， 16， 4）都是允许的 blockDim值。网格可以比线程块有更高的维度，反之亦然。可以使用下面的主机代码生成下图的网络：
dim3 dimGrid(2, 2, 1); dim3 dimBlock(4, 2, 2); KernelFunction<<<dimGrid, dimBlock>>>(...);
这个网络包含 4 个线程块，每个线程块可以通过（blockIdx.y, blockIdx.x）标记。例如，（1，0）位置上线程块的 blockIdx.y=1，blockIdx.x=0。每个线程块由 4x2x2 的线程数组组成。例如线程(1, 0, 2)的 threadIdx.z=1, threadIdx.y=0, threadIdx.z=2。注意标签是按照 z, y, x 的顺序来的。
2. 线程与多维数组的映射
线程组织形式的设置要依据数据的特点，灰度图片为二维数组，那就把网络和线程块组织成二维形式处理灰度图片。
以下图为例，对于一张 72x62 的图片用 16x16 的线程块处理，则在 x，y 方向分别需要 5 个和 4 个线程块。但 x，y 方向上线程块分别有 4 个线程和 2 个线程多余，需要像 vecAddKernel 函数一样用 if 语句防止多余的线程进行运算。
//因为标签是按照 y, x 的顺序来的，所以gird中参数是ceil(m/16.0), ceil(n/16.0), 1。 dim3 dimGrid(ceil(m/16.0), ceil(n/16.0), 1); dim3 dimBlock(16, 16, 1); colorToGreyscaleConversion<<<dimGrid,dimBlock>>>(d_Pin,d_Pout,m,n);
对于一张 2000 × 1500 (3-M个像素) 图片，将会产生 11,750 个线程块，x 方向上 125 个，y 方向上 94 个。在这个 kernel 函数中，gridDim.x, gridDim.y, blockDim.x 和 blockDim.y 的值分别为 125, 94, 16 和 16。
对于动态内存分配的二维数组，列数无法确定，然而 CUDA C 编译二维数组时需要知道列数，所以我们还用等价的一维数组表示二维数组。一维数组表示高维数组有行序为主序和列序为主序两种方式。下图为行序为主序的示意图。
C 编译器是以行序为主序线性化二维数组的，但 FORTRAN 是以列序为主序， 用 CUDA FORTRAN 编程要注意。
原理： L = r * 0.21 + g * 0.72 + b * 0.07
结合下图可以看出 Pout 的列索引 Col 可以通过 threadIdx.x + blockIdx.x * blockDim.x 求得；行索引 Row 可以通过 threadIdx.y + blockIdx.y * blockDim.y 求得。而 Pout 以行序为主序在计算机中存储，那么一个像素块在 Pout 中的索引就为：Rowwidth + Col。但 Pin 为 RGB 图像，在以为数组中的存储形式为相邻三个元素分别为一个像素 r, g, b 三通道的值，那么一个像素块的 r 通道在 Pout 中的索引就为：(Row*width + Col) \ 3，g, b 分别加1，2。
colorToGreyscaleConversion() 函数源代码
// we have 3 channels corresponding to RGB
// The input image is encoded as unsigned characters [0, 255]
__global__
void colorToGreyscaleConversion(unsigned char * Pout, unsigned char * Pin, int width, int height){
int Col = threadIdx.x + blockIdx.x * blockDim.x;
int Row = threadIdx.y + blockIdx.y * blockDim.y;
if(Col < width && Row < height){
// get 1D coordinate for the grayscale image
int greyOffset = Row*width + Col;
// one can think of the RGB image having CHANNEL times columns than the grayscale image
int rgbOffset = greyOffset * CHANNELS;
unsigned char r = Pin[rgbOffset ]; // red value for pixel
unsigned char g = Pin[rgbOffset + 1]; // green value for pixel
unsigned char b = Pin[rgbOffset + 2]; // blue value for pixel
Pout[greyOffset] = 0.21f*r + 0.71f*g + 0.07f*b;
}
} 3. 矩阵乘法 kernel 函数
以下图为例看一下矩阵的乘法，d_M 矩阵的第 Row 行与 d_N 矩阵的第 Col 列内积得到 d_P 矩阵的一个元素 d_P[Row][Col]
使每一个线程计算 d_P 中的一个元素，每个线程在网格中的第 blockIdx.y * blockDim.y + threadIdx.y 行，第 blockIdx.x * blockDim.x + threadId.y 列。以下为代码：
__global__ void MatrixMulKernel(float *d_M, float *d_N, float *d_P, int m, int k, int n){
int Row = blockIdx.y*blockDim.y+threadIdx.y;
int Col = blockIdx.x*blockDim.x+threadIdx.x;
float PValue = 0;
if(Row < m && Col < n){
for (int i = 0; i < k; ++i)
PValue += A[Row * k + i] * B[i * n + Col];
d_P[Row * n + Col] = PValue;
}
} 启动 MatrixMulKernel() 函数的主机代码：
#define BLOCK_SIZE 32
void cuda(float *A, float *B, float *C, int m, int k, int n) {
float *d_A, *d_B, *d_C;
size_t size = m * k * sizeof(float);
cudaMalloc((void**)&d_A, size);
cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
size = k * n * sizeof(float);
cudaMalloc((void**)&d_B, size);
cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);
// Allocate C in device memory
size = m * n * sizeof(float);
cudaMalloc((void**)&d_C, size);
// Invoke kernel
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid(ceil(n/float(blockDim.x), ceil(m/float(blockDim.y));
MatrixMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C, m, k, n);
cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
// Free device memory
} 4. 线程同步和透明的可拓展性
CUDA 中需要协调多个线程的执行过程，CUDA 允许同一个线程块里的线程使用同一个栅栏函数 syncthreads() 协调执行过程。当 kernel 函数执行到 __syncthread() 函数时，线程块中的线程执行到调用位置时需要等待所有线程执行到该点才能继续执行。下图是栅栏函数的执行示意图：
在CUDA中,syncthreads()语句(如果存在)必须由块中的所有线程执行.当a syncthreads()放在一个if语句中时,块中的所有线程都会执行包含syncthreads()或不包含它们的路径.对于if-then-else语句,如果每个路径都有一个syncthreads()语句,则块中的所有线程都在路径syncthreads()上执行,或者所有线程then都执行该else路径.这两个__syncthreads()是不同的屏障同步点.如果块中的线程执行then路径而另一个线程执行else路径,则它们将在不同的屏障同步点处等待.他们最终会永远等待对方.程序员有责任编写代码以满足这些要求。
假设有一个内核,该内核使用一个由32个线程组成的线程块启动.
kernel<<<1,32>>>()
//内核的代码如下:
__global__ void kernel()
{
if (threadIdx.x < 16)
{
// do something
__syncthreads();
}
else
{
// do something
__snycthreads();
}
} 线程块的前16个线程将运行if语句.另外16个else语句.如果前16个线程中的每个线程都到达__syncthreads,则它们将阻塞,直到整个线程块到达该语句.但是这种情况永远不会出现,因为其他16个线程卡在了else分支中,并且会出现死锁.
你应该避免__syncthreads在不同的if和else分支中使用,或者你必须确保整个线程块在同一个分支中运行。
5. 线程块的资源分配
当前硬件中，执行资源被组织成多核流处理器 (Streaming Multiprocessor, SM) ，几个线程块分配到一个 SM 上，每个设备对一个 SM 能分配的线程块是有限制的。下图每个 SM 分配了 3 个线程块。
最新的 CUDA 设备中一个 SM 最多可以分配 1536 个线程。如果一个 SM 最多可以分配 8 个线程块，那么可以分配 6 个包含 256 线程的线程块，3个包含 512 线程的线程块。
6.查询 GPU（设备） 属性
#include<cuda.h>
#include<stdio.h>
int main(){
int dev_count;
cudaGetDeviceCount(&dev_count); //获取GPU的数目
printf("number of GPU: %d\n", dev_count);
cudaDeviceProp dev_prop;
for(int i = 0; i < dev_count;++i){ //遍历每个GPU，并获取部分参数
cudaGetDeviceProperties(&dev_prop, i);
printf("GPU %d 的线程块最大线程数: x = %d, y = %d, z = %d\n", i+1, dev_prop.maxThreadsDim[0], dev_prop.maxThreadsDim[1], dev_prop.maxThreadsDim[2]);
printf("GPU %d 的网格最大线程块数量: x = %d, y = %d, z = %d\n", i+1, dev_prop.maxGridSize[0], dev_prop.maxGridSize[1], dev_prop.maxGridSize[2]);
printf("GPU %d 的warpSize: %d\n", i+1, dev_prop.warpSize);
}
return 0;
} 老师分配服务器的输出结果：
7. 线程调度和容许时延
某个线程块分配给 某个SM 后，线程块将被划分为由 32 个线程组成的 warp 。下午展示了3个线程块的划分情况，每个 warp 包含 32 个连续的线程：第一个 warp 包含 0-31 号线程，第二个包含 32-63 号线程，以此类推。
SM 按照 SIMD 的模式执行一个 warp 中的所有线程，warp 中的所有线程只能取一条指令执行。
上图也说明了一个 SM 中包含多个流处理器 (Streaming Processor, SP), 这些 SP 就是真正执行指令的部件。一般情况下，SM 中 SP 的数目小于分配到的线程的数目，也就是说每个 SM 中只有一部分线程执行指令。早期 GPU 中在任意给定时刻，每个 SM 只能执行一个 warp 的一条指令，但当前 GPU 中每个 SM 能同时执行几个 warp 的指令。不管那种情况， SM 中的硬件只执行所有 warp 的一部分，这是为了提高执行长延迟操作的效率。
当一个 warp 的一条指令需要等待先前启动的长延迟操作的结果时，这个 warp 将不会被选中执行。系统将调度另一个空闲的 warp 执行运算。
参考：大规模并行编程处理器实战（第二版）David B. Kirk, Wen-mei W. Hwu 著，赵开勇 汪朝晖 程亦超 译
Programming Massively Parallel Processors A Hands-on Approach（3-rd）David B. Kirk, Wen-mei W. Hwu
https://qa.1r1g.com/sf/ask/1442942231/
摩根士丹利二面2人在聊
查看12道真题和解析
