📚 课程目录第4章 GPU编程入门

4. GPU编程入门：从架构到Hello World

好，咱们正式进入GPU编程。我知道很多做生物信息的朋友，平时跟CPU打交道多，突然听说要用GPU来跑BLAST，第一反应可能是：“这东西跟CPU有啥区别？我写代码的方式要变多少？”

嗯，今天我就从底层给你拆开看看。GPU不是“更快更强的CPU”，它是一台异构的协处理器。你写的主机代码跑在CPU上，而计算密集的部分——比如序列比对里的得分矩阵填充——会被扔到GPU上并行执行。

4.1 GPU架构：一个视觉上的“暴力美学”

先看架构。我当年第一次接触CUDA时，最震撼的不是性能，而是芯片布局。CPU里，缓存和控制单元占了大部分面积，真正做算术的ALU（算术逻辑单元）就那么几个。GPU呢？反过来的。

一句话总结：GPU靠成百上千个“小核”堆规模，CPU靠几个“大核”拼单核性能。BLAST这种数据并行（每个查询序列独立计算）的任务，是GPU的菜。

4.2 CUDA核心 vs SM vs 显存层次

记东西不能死记。我习惯从调度视角理解这张图：

CUDA Core：就是个简单的浮点/整数运算单元。没有复杂的乱序执行，没有分支预测。纯干活。
SM（Streaming Multiprocessor）：一个SM里有几十到上百个CUDA Core，还有一块共享内存（Shared Memory）和一大把寄存器。SM负责调度执行线程。
显存层次：Global Memory（显存）最大但最慢，Shared Memory很小但跟L1缓存差不多快。你在BLAST里查数据库，通常数据先放Global Memory，每个block再把需要的小块拉到Shared Memory里复用。我最早做GPU版Smith-Waterman时，没用好Shared Memory，结果性能只比CPU快两倍……后来把得分矩阵的常用列缓存到Shared Memory，直接提升了8倍。

存储层次	大小	延迟	作用域
寄存器 (Register)	每SM几千个	~1 cycle	单个线程
共享内存 (Shared Mem)	48~96 KB / SM	~20 cycles	单个Block
L2缓存	几MB	~100 cycles	所有SM
全局内存 (Global Mem)	几~几十GB	~400 cycles	所有线程 + CPU

我个人习惯：凡是需要在Block内线程之间共享的数据，优先放Shared Memory。比如BLAST的得分矩阵的滑动窗口，用Shared Memory能省下全局内存带宽，这是我调优时最常用的手法。

4.3 CUDA编程模型：比你想象中简单

你写一个kernel函数——就是跑在GPU上的函数——前面加上 __global__ 关键字。然后在主机端用 <<<gridDim, blockDim>>> 指定线程网格和线程块大小，调用它。CPU跟GPU之间用 cudaMemcpy 搬数据。

嗯，举个例子你会更清楚。我带你写一个“Hello World”，但让它真正干点活——让1000个线程各自打印自己的线程ID，同时顺便计算一个简单的加法。

4.4 Hello World 示例（带干活版）

#include <stdio.h>
#include <cuda_runtime.h>

// kernel：每个线程执行一次
__global__ void hello_cuda(int *out, int N) {
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    if (idx < N) {
        // 给每个线程分配一个独一无二的ID
        out[idx] = idx * 2;   // 随便算点东西
        // printf 在GPU里也可以，但尽量少用，影响性能
        printf("Thread %4d (block %2d, local %2d)\n",
               idx, blockIdx.x, threadIdx.x);
    }
}

int main() {
    int N = 1000;
    int bytes = N * sizeof(int);

    // 1. 分配GPU内存
    int *d_out;
    cudaMalloc(&d_out, bytes);

    // 2. 启动kernel：10个block，每个block 100个线程
    // 总共 10 * 100 = 1000 个线程，正好覆盖N
    hello_cuda<<<10, 100>>>(d_out, N);

    // 3. 把结果拷回主机
    int *h_out = (int*)malloc(bytes);
    cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost);

    // 4. 验证
    for (int i = 0; i < 10; i++) {
        printf("h_out[%d] = %d\n", i, h_out[i]);
    }

    // 5. 清理
    cudaFree(d_out);
    free(h_out);
    return 0;
}

编译命令（假设你CUDA Toolkit装好了，并且nvcc在PATH里）：

nvcc -o hello_cuda hello_cuda.cu
./hello_cuda

我曾经踩过的坑：<<<10, 100>>> 的数字乘积不能大于你的GPU能同时执行的线程数上限（通常 blocks * threads_per_block 不超过65535*1024？实际要看Compute Capability）。另外，printf 在kernel里非常慢——我调试时常用，但生产代码一定要删掉，不然性能直接腰斩。

4.5 代码里藏的几个关键点

内置变量：threadIdx.x 是线程在block内的索引，blockIdx.x 是block在grid内的索引，blockDim.x 是block的线程数。你用这三个量就能算出全局唯一ID。
内存拷贝：cudaMemcpyDeviceToHost 是方向，别搞反了。搞反了不会报错，但拷回来的是一堆垃圾——我就干过这事，排查了一下午。
cudaMalloc：注意参数是设备指针的指针。你要是写成 cudaMalloc(d_out, bytes)，直接段错误——因为d_out被当作地址值而非指针。

4.6 线程层次：Grid → Block → Thread

你看下面这个关系，说白了就是三层嵌套：

一个Grid = 多个Block
一个Block = 多个Thread
一个Thread 跑你写的kernel函数中的一份代码

实际映射到硬件上，一个Block会调度到一个SM上执行，Block里的线程被分成Warp（32个线程一组）。SM以Warp为单位发射指令。如果Warp里32个线程走了不同的分支（比如if-else），它们就会串行化，效率变低。这是BLAST加速时你尤其要关注的——尽量保证同一个Warp的线程执行相同路径。我在优化序列比对kernel时，就通过重新组织输入数据，让同一个Warp处理相似长度的序列，分支发散减少了70%。

4.7 本章小结（用自己的话说）

GPU编程入门就这么回事：

理解SM + CUDA Core + 显存层次——这是硬件基础。
学会用 __global__ 写kernel。
会用 <<<grid, block>>> 启动线程。
会用 cudaMemcpy 在主机和设备之间搬数据。