《大规模并行处理器程序设计》学习笔记

发表于 2026-09-23 分类于笔记阅读次数：

Heterogeneous data parallel computing

CUDA C编程模型：Single-program mutiple-data (SPMD).

cudaMalloc:

函数原型：

1	cudaError_t cudaMalloc(void **devPtr, size_t size);

使用：

1 2	float A_d; cudaMalloc((void )&A_d, n sizeof(float));

cudaMemcpy：

cudaMemcpy(destination, source, size, direction);
cudaMemcpy(A_d, A_h,
           n * sizeof(float),
           cudaMemcpyHostToDevice);

错误处理：

cudaError_t err;

err = cudaMalloc(&A_d, n * sizeof(float));

if (err != cudaSuccess) {
    printf("cudaMalloc failed: %s\n",
           cudaGetErrorString(err));
}

CUDA 的线程组织层级是：

Grid（网格）
 ├── Block 0（线程块）
 │    ├── Thread 0
 │    ├── Thread 1
 │    └── ...
 ├── Block 1
 └── ...

All blocks are of the same size. Up to 1024 threads.

一次 kernel launch 会产生一个 Grid。

blockDim: 表示一个线程块中各个维度的线程数量。

全局唯一索引的写法

1	int i = threadIdx.x + blockDim.x * blockIdx.x;

Cuda C 关键字

声明	在哪里执行	从哪里调用
普通函数 / `__host__`	CPU	CPU
`__device__`	GPU	GPU
`__global__`	GPU	CPU（通常）
`__host__ __device__`	CPU 或 GPU	CPU 或 GPU

内核调用和网格启动：kernel<<<gridDim, blockDim>>>(参数);

Multidimensional grids and data

指定grid和block大小：

1
2
3

dim3 dimGrid(32,1,1);
dim3 dimBlock(128,1,1);
vecAddKernel<<<dimGrid, dimBlock>>>(...);

多维数组要转化成一维数组实现

// 分配一个 2D 数组，行宽 1024 个 float，高 768
size_t pitch;
float *d_data;
cudaMallocPitch(&d_data, &pitch, 1024 * sizeof(float), 768);

// 核函数内访问 (row, col)
__global__ void kernel(float *data, size_t pitch, int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (col < width && row < height) {
        // 注意：这里的 pitch 是字节数，需要先转为元素个数
        int index = row * (pitch / sizeof(float)) + col;
        data[index] = row * 0.5f + col;
    }
}

Compute architecture and scheduling

GPU由SM构成。

不同Block中的线程无法相互同步。

线程以块为单位分配给 SM 执行。一旦某个块被分配给 SM，它会被进一步划分为warps。

warps中的模型遵循SIMP，因此分支会拖慢性能。

SM 的占用率越高，其隐藏长延迟操作的能力就越强。

Heterogeneous data parallel computing

算数强度：从全局内存中访问的浮点运算（FLOP）与字节（B）的比率。

分析程序潜在性能：屋顶线模型。

寄存器速度快的原因：访问延迟短，所需指令数少，消耗能量低。

Variable declaration	Memory	Scope	Lifetime
Automatic variables other than arrays	Register	Thread	Grid
Automatic array variables	Local	Thread	Grid
`__device__ __shared__ int SharedVar;`	Shared	Block	Grid
`__device__ int GlobalVar;`	Global	Grid	Application
`__device__ __constant__ int ConstVar;`	Constant	Grid	Application

Local Memory物理上是放在device memory上的。寄存器不够用发生register spilling时，放在Local Memory。按照经验，很少需要使用自动数组变量。

使用大量寄存器可能会对每个 SM 的占用率产生负面影响。

策略：将数据划分为称为块（tile）的子集，使每个块都能放入共享内存。

#define TILE_WIDTH 16
// blockDim.x = TILE_WIDTH
// blockDim.y = TILE_WIDTH
__global__ void matrixMulKernel(float* M, float* N, float* P, int Width) {

    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int bx = blockIdx.x;  int by = blockIdx.y;
    int tx = threadIdx.x; int ty = threadIdx.y;

    // Identify the row and column of the P element to work on
    int Row = by * TILE_WIDTH + ty;
    int Col = bx * TILE_WIDTH + tx;

    // Loop over the M and N tiles required to compute P element
    float Pvalue = 0;
    for (int ph = 0; ph < Width / TILE_WIDTH; ++ph) {

        // Collaborative loading of M and N tiles into shared memory
        Mds[ty][tx] = M[Row * Width + ph * TILE_WIDTH + tx];
        Nds[ty][tx] = N[(ph * TILE_WIDTH + ty) * Width + Col];

        __syncthreads();

        for (int k = 0; k < TILE_WIDTH; ++k) {
            Pvalue += Mds[ty][k] * Nds[k][tx];
        }

        __syncthreads();
    }

    P[Row * Width + Col] = Pvalue;
}

分phase的思想称为strip-mining

分块之后，全局内存访问次数减少了 TILE_WIDTH 倍

边界检查

每次内存访问都需要有相应的检查，以确保访问中使用的索引在正在访问的数组的边界之内。

// Loop over the M and N tiles required to compute P element
float Pvalue = 0;
for (int ph = 0; ph < ceil(Width/(float)TILE_WIDTH); ++ph) {

    // Collaborative loading of M and N tiles into shared memory
    if ((Row < Width) && (ph*TILE_WIDTH+tx) < Width)
        Mds[ty][tx] = M[Row*Width + ph*TILE_WIDTH + tx];
    else
        Mds[ty][tx] = 0.0f;

    if ((ph*TILE_WIDTH+ty) < Width && Col < Width)
        Nds[ty][tx] = N[(ph*TILE_WIDTH + ty)*Width + Col];
    else
        Nds[ty][tx] = 0.0f;

    __syncthreads();

    for (int k = 0; k < TILE_WIDTH; ++k) {
        Pvalue += Mds[ty][k] * Nds[k][tx];
    }

    __syncthreads();
}

if ((Row < Width) && (Col < Width))
    P[Row*Width + Col] = Pvalue;

动态调整参数

1	extern __shared__ Mds_Nds[];

可以根据设备查询结果动态配置每个块要使用的共享内存量。

核心代码：

__global__ void MatrixMulKernel(float *M, float *N, float *P, int Width)
{
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    int Row = blockIdx.y * TILE_WIDTH + ty;
    int Col = blockIdx.x * TILE_WIDTH + tx;

    // 动态共享内存
    extern __shared__ float Mds_Nds[];

    // 前半部分给 Mds
    float *Mds = Mds_Nds;

    // 后半部分给 Nds
    float *Nds = Mds_Nds + TILE_WIDTH * TILE_WIDTH;

    float Pvalue = 0.0f;

    for (int ph = 0;
         ph < ceil(Width / (float)TILE_WIDTH);
         ++ph) {

        // 加载 M tile
        if ((Row < Width) &&
            (ph * TILE_WIDTH + tx < Width))
            Mds[ty * TILE_WIDTH + tx]
                = M[Row * Width
                    + ph * TILE_WIDTH + tx];
        else
            Mds[ty * TILE_WIDTH + tx] = 0.0f;

        // 加载 N tile
        if ((ph * TILE_WIDTH + ty < Width) &&
            (Col < Width))
            Nds[ty * TILE_WIDTH + tx]
                = N[(ph * TILE_WIDTH + ty)
                    * Width + Col];
        else
            Nds[ty * TILE_WIDTH + tx] = 0.0f;

        __syncthreads();

        for (int k = 0; k < TILE_WIDTH; ++k) {
            Pvalue +=
                Mds[ty * TILE_WIDTH + k]
                *
                Nds[k * TILE_WIDTH + tx];
        }

        __syncthreads();
    }

    if (Row < Width && Col < Width)
        P[Row * Width + Col] = Pvalue;
}

CPU:

dim3 dimBlock(TILE_WIDTH, TILE_WIDTH);
dim3 dimGrid(
    ceil(Width / (float)TILE_WIDTH),
    ceil(Width / (float)TILE_WIDTH)
);

size_t sharedMemSize =
    2 * TILE_WIDTH * TILE_WIDTH * sizeof(float);

MatrixMulKernel<<<dimGrid, dimBlock, sharedMemSize>>>(
    M, N, P, Width
);

Performance considerations

最有利的访问模式是当线程束中的所有线程访问连续的全局内存位置时实现的。

Memory coalescing: cuda自动将访问连续内存的指令合并。

例：

#define TILE_WIDTH 32
#define COARSE_FACTOR 4

__global__ void matrixMulKernel(float* M, float* N, float* P, int width)
{
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    // Identify the row and column of the P element to work on
    int row = by * TILE_WIDTH + ty;
    int colStart = bx * TILE_WIDTH * COARSE_FACTOR + tx;

    // Initialize Pvalue for all output elements
    float Pvalue[COARSE_FACTOR];
    for (int c = 0; c < COARSE_FACTOR; ++c) {
        Pvalue[c] = 0.0f;
    }

    // Loop over the M and N tiles required to compute P element
    for (int ph = 0; ph < width / TILE_WIDTH; ++ph) {

        // Collaborative loading of M tile into shared memory
        Mds[ty][tx] = M[row * width + ph * TILE_WIDTH + tx];

        for (int c = 0; c < COARSE_FACTOR; ++c) {

            int col = colStart + c * TILE_WIDTH;

            // Collaborative loading of N tile into shared memory
            Nds[ty][tx] = N[(ph * TILE_WIDTH + ty) * width + col];
            __syncthreads();

            for (int k = 0; k < TILE_WIDTH; ++k) {
                Pvalue[c] += Mds[ty][k] * Nds[k][tx];
            }

            __syncthreads();
        }
    }

    for (int c = 0; c < COARSE_FACTOR; ++c) {
        int col = colStart + c * TILE_WIDTH;
        P[row * width + col] = Pvalue[c];
    }
}

每个warp中的线程k相同，col相邻，因此对M的访问相邻，可以合并。

Corner Turning: 先按照有利于 Global Memory 合并访存的方向把数据读进 Shared Memory，再在 Shared Memory 中换一个方向访问。

Hiding memory latency： 当一个 warp 因为等待内存数据而停住时，让 SM 去执行别的 warp，这样计算单元尽量不要闲着。故高Occupancy 可以带来更多latency hiding。

Thread coarsening: 让一个线程不再只计算一个输出，而是一次计算多个输出。

收益来源：数据复用

如计算矩阵乘法时，一个线程算4个值：

P[row][colStart]
P[row][colStart + TILE_WIDTH]
P[row][colStart + 2*TILE_WIDTH]
P[row][colStart + 3*TILE_WIDTH]

要权衡复用收益和并行度损失。

例：

#define TILE_WIDTH 32
#define COARSE_FACTOR 4

__global__ void matrixMulKernel(float* M, float* N, float* P, int width)
{
    __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
    __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    int bx = blockIdx.x;
    int by = blockIdx.y;
    int tx = threadIdx.x;
    int ty = threadIdx.y;

    // Identify the row and column of the P element to work on
    int row = by * TILE_WIDTH + ty;
    int colStart = bx * TILE_WIDTH * COARSE_FACTOR + tx;

    // Initialize Pvalue for all output elements
    float Pvalue[COARSE_FACTOR];
    for (int c = 0; c < COARSE_FACTOR; ++c) {
        Pvalue[c] = 0.0f;
    }

    // Loop over the M and N tiles required to compute P element
    for (int ph = 0; ph < width / TILE_WIDTH; ++ph) {

        // Collaborative loading of M tile into shared memory
        Mds[ty][tx] = M[row * width + ph * TILE_WIDTH + tx];

        for (int c = 0; c < COARSE_FACTOR; ++c) {

            int col = colStart + c * TILE_WIDTH;

            // Collaborative loading of N tile into shared memory
            Nds[ty][tx] = N[(ph * TILE_WIDTH + ty) * width + col];
            __syncthreads();

            for (int k = 0; k < TILE_WIDTH; ++k) {
                Pvalue[c] += Mds[ty][k] * Nds[k][tx];
            }

            __syncthreads();
        }
    }

    for (int c = 0; c < COARSE_FACTOR; ++c) {
        int col = colStart + c * TILE_WIDTH;
        P[row * width + col] = Pvalue[c];
    }
}

优化清单

优化	对compute cores的好处	对memory的好处	策略
Maximizing occupancy	更多的计算任务可以隐藏pipline latency	更多的并行内存访问可以隐藏DRAM延迟	调整 SM资源的使用，如每块的线程数、每块的共享内存量和每线程的寄存器数。使得线程数远多于核心数。
Coalesced global memory accesses	减少因等待全局内存访问而产生的流水线停顿	减少全局内存流量，并更好地利用突发传输/缓存行	corner turning或类似操作；重新排列线程到数据的映射方式；重新排列数据的布局。
Minimizing control divergence	高SIMD效率	-	重新安排线程与工作、数据的关系，和数据的布局
Tiling of reused data	减少因等待全局内存访问而产生的pipeline stalls	降低全局内存流量	将块内重复使用的数据放置在共享内存或寄存器中，仅在全局内存与 SM之间传输一次
Privatization	减少因等待atomic updates而产生的流水线停顿	降低原子更新的竞争和串行化	对数据的私有副本进行部分更新，完成后才更新全局副本
Thread coarsening	减少冗余工作、分支发散或同步	减少冗余的全局内存流量	为每个线程分配多个并行单元

卷积

特点：输出元素独立，容易并行。但不同输出会反复读取重叠区域。性能受到带宽限制。

朴素实现：一个线程计算一个输出元素。

优化1：将fliter放入constant cache。

优化2：一个block一起加载input tile到shared memory，片上复用，减少重复加载。

优化3：tile内部使用shared memory，halo直接全局访问，期望halo已经存在L2缓存中。

Stencil

Stencil sweep的特点是使用特殊卷积核的卷积，工作在3D网络上，分块容易碰到线程块大小和共享内存容量限制。

很多时候stencil只使用中心点和坐标轴上的邻居，因此可以将部分输入放到寄存器中，不用放在shared memory。

数值精度：用来模拟偏微分方程，通常要FP64

朴素写法：一个线程负责一个输出。

优化1：将tile搬到shared memory。边界片段，只有映射到tile内部线程参与计算，其余负责加载halo。halo: tile为计算内部输出而额外加载的邻域。

优化2：Thread coarsening。沿着z方向计算一列点，每个线程加一个z方向循环。维护过去、现在、将来三个平面在shared memory中。

优化3：Register tiling，prev和next使用率不高，没必要放在shared memory，可以直接依赖寄存器。

原子操作

当多个线程更新一个输出元素时，会发生竞争。

优化1：使用atomicAdd保证正确性，但竞争严重。

优化2：使用块级私有化，将高竞争的输出结构复制成多个私有副本，让不同线程子集更新不同副本，最后再将副本合并到公共结果中。

收益：

将大范围竞争缩小到了线程子集内部。
副本可以放在shared memory实现低延迟。
大幅提升更新吞吐量

代价：

额外的空间
额外的副本初始化开销
额外的合并开销

优化3：使用Coarsening，减少副本数、合并和初始化次数。使用交错分区实现合并访存。

优化4：Aggregation：合并连续的同-bin更新，先积累在寄存器里，最后再提交，进一步减少原子操作次数。

Reduction

任务：将一串输入值的输入通过一棵二叉归约树压缩成一个值。

朴素实现：每轮将元素数减半。问题：线程映射导致非合并访存和Control Divergence。

收敛式线程映射：让活跃线程保持连续，可以显著降低全局显存请求数，但流量仍偏大。

共享内存规约：只在开始读全局内存，中间结果放在shared memory中。

Hierarchical Reduction：将大输入拆成互补依赖的segment，每个block分别运算，最后将结果累加到全局输出。

Thread Coarsening: 每个线程先独立串行累加，再跑规约树。

Prefix Sum

Inclusive scan：包含当前位置的前缀和。Exclusive scan：排除当前位置的前缀和。

Kogge–Stone并行Scan：

一个block负责一个输入section
每个线程将一个元素从全局搬到shared memory中
在block内执行轮 Kogge–Stone
结果写回global memory

Brent–Kung：

Phase 1：2, 4, 8…的间隔更新，得到部分完整前缀
Phase2：使用完整前缀补齐其余位置

Coarsening：

phase 1：线程内部顺序scan
phase2：对子段和做block-wide scan
phase3：将偏移量加回各子段

处理任意长度的输入：

前面的算法只能在一个block中扫表，大规模输入的话使用三步层次化方法：

Kernel 1：各block做局部scan，section最后一个输出就是该section的总和，将其写入全局数组S[b]。
Kernel 2：扫描section总和数组S，若S仍大于单block可以处理的范围，递归使用同样的分层结构。
Kernel 3：Uniform Add，对于block b的每个局部结果加上S[b-1]

Domino-style Scan: 全局累计偏移量不一定需要一次grid-wide barrier，可以让它从左向右传递。

Merge

并行化merge算法：

按照输出的范围划分输出，对两个输出边界分别调用co-rank
输出分界位置 k称为 rank。其co-rank是唯一的一对 (i, j)，满足i+j=k，即输出的前k个元素，恰好由A的前i个元素与B的前j个元素稳定归并得到。
co-rank可以通过二分计算
问题：相邻线程处理相隔较远的输入子数组，同一条指令下的读取通常不连续；相邻线程写入不同的C子段，同步执行时地址也可能跨较大步长；每个线程的co-rank都在global memory上做不规则二分访问；因此难以合并访存

Tiled Merge Kernel：利用shared memory执行co-rank。

Circular-buffer Merge Kernel：不丢弃 tile 中未消费的元素，新一轮只补充上一轮消耗掉的数量，写指针到达数组末尾后回绕到开头。