2. GPU计算基础：GPU架构简介、CUDA编程模型与线程层次

好，让我们开始真正接触GPU的世界。我见过不少搞生信的朋友，一开始觉得GPU就是个“更快的显卡”，装上就能让BLAST飞起来。结果呢？要么跑不起来，要么效率还不如CPU。为什么？

说白了，GPU和CPU是两种完全不同的思维方式。你想想看，CPU是个“全能选手”，什么活都能干，但一次只能干几件。GPU呢，是个“搬砖大队”，每个工人能力一般，但胜在人多，能同时干成千上万件简单的活。BLAST这种序列比对任务，本质上就是海量的独立计算——太适合GPU了。

2.1 GPU架构：它到底长什么样？

我们先看看GPU的核心部件。下面这张图是我个人习惯用来跟团队新人讲架构的简化模型，你一看就明白。

这张图我每次培训都会用。左边是物理硬件，右边是编程模型，中间是映射关系。很多人搞不清楚这两者的区别，总觉得“线程”就是“核心”。不是的，线程是逻辑概念，核心是物理资源。一个核心可以执行成千上万个线程，靠的是极速的上下文切换——这点和CPU很不一样。

2.2 CUDA编程模型：我该怎么写代码？

CUDA编程模型，说白了就是NVIDIA给你的一套“指挥GPU干活”的工具。你写一个函数，叫kernel（内核函数），然后告诉GPU：“嘿，用10000个线程去跑这个函数！”

一个最简单的CUDA程序长啥样？我直接上代码：

// 一个简单的向量加法kernel
__global__ void vecAdd(float *A, float *B, float *C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    // ... 分配内存，拷贝数据到GPU ...
    int N = 1024;
    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    // 启动kernel
    vecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
    // ... 拷贝结果回CPU ...
}

你看，关键是那行vecAdd<<<blocksPerGrid, threadsPerBlock>>>>(...)。尖括号里的两个参数，就是告诉GPU：我要多少个线程块（Block），每个块里有多少个线程（Thread）。这就是接下来要细聊的线程层次。

2.3 线程层次：理解Grid、Block、Thread

线程层次是CUDA的核心。它分三级，从大到小：

• Grid（网格）：一次kernel启动对应的所有线程。可以是一维、二维或三维。
• Block（线程块）：一组可以被SM独立调度的线程。同一个Block内的线程可以共享数据（通过共享内存）。
• Thread（线程）：最小的执行单元。每个线程都有自己独立的寄存器、程序计数器、局部内存。

为什么分这么多级？你想想看，如果只有“线程”这级，你怎么调度1万个线程？硬件并不知道哪些线程可以一起干活。有了Block这个中间层，硬件就能以Block为单位，把任务分配给不同的SM去并行执行。

还有一个概念特别重要——Warp。这是GPU调度的基本单位，一个Warp固定包含32个线程。SM每次发出一条指令，这个Warp里的32个线程要么一起执行，要么有些因为条件分支被“屏蔽”了。我在项目中遇到过因为线程分支导致性能断崖下跌的情况，后来查出来是warp divergence（线程束分化）在作祟。

2.4 内存层次：数据放哪儿，是个学问

GPU的内存层次很丰富，用错了位置，性能就全毁了。

嗯，这里要注意。全局内存是最常用的，但也是最慢的。BLAST这种应用，如果频繁访问全局内存，性能肯定上不去。我以前优化过一个数据库加速项目，一开始就是无脑用全局内存，结果GPU利用率连20%都不到。后来把频繁访问的热数据搬到共享内存里，利用率直接飙到70%+。

2.5 CUDA的hello world：从代码到启动

来，我们看一个更完整的例子，把刚才的线程层次和内存层次串起来：

#include <stdio.h>

// 使用共享内存的向量求和
__global__ void sumWithShared(const float *A, const float *B, float *C, int N) {
    __shared__ float cache[256];

    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    if (tid < N) {
        cache[threadIdx.x] = A[tid] + B[tid];
    }
    __syncthreads();  // 同步，确保所有线程写完共享内存

    if (threadIdx.x == 0) {
        float sum = 0.0f;
        for (int i = 0; i < blockDim.x; i++) {
            sum += cache[i];
        }
        C[blockIdx.x] = sum;
    }
}

int main() {
    int N = 1024 * 1024;
    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    // 分配、初始化、拷贝（略）
    float *d_A, *d_B, *d_C;
    cudaMalloc(&d_A, N * sizeof(float));
    // ... cudaMemcpy ...

    sumWithShared<<<blocksPerGrid, threadsPerBlock>>>>(d_A, d_B, d_C, N);

    cudaDeviceSynchronize();  // 等GPU干完
    // ... cudaMemcpy回CPU ...
}

__syncthreads()这行很关键。它保证一个Block内的所有线程都执行到这一行了，才能继续。不然，线程A刚写完共享内存，线程B就开始读，数据就乱套了。