GPU Gems 2 · Chapter 43. GPU Computing for Protein Structure Prediction

蛋白质结构预测里，实验设备通常只给出一部分原子之间的距离范围。剩下的范围很宽，像一张只画出少量道路的地图：单条道路的信息不够，但三点之间的绕行关系可以不断排除不可能的距离。本章把这个“收紧范围”的过程变成 GPU 可并行的矩阵更新。

GPU 适合这里不是因为它懂生物学，而是因为每一轮都能让很多矩阵单元同时读取相同的中间节点并写回自己的位置。难点在于更新有时间顺序、上下界有依赖、纹理不能原地读写；处理好这三件事，旧式 fragment unit 也能承担密集的数值计算。

1. 从原子距离到带边界的完全图

distance-bound smoothing 是本章的任务。将含有 n 个原子的分子表示为完全图 Kₙ：节点是原子，边 (i, j) 记录距离的上界 uᵢⱼ 和下界 lᵢⱼ。NMR 只提供部分边的约束，未知边可以先用分子直径和接近零的下界占位。

矩阵存储让 GPU 可以把每条边当作一个独立的输出单元。距离矩阵的对角线是零，通常还保持对称性 D[i, j] = D[j, i]。不过，GPU 版本将上下界分别放在纹理里，即使 (i, j) 与 (j, i) 重复，也能减少 shader 的条件分支和取样次数。

2. 三角不等式如何收紧范围

triangle inequality 给每条边带来一个中间节点 k。对上界，绕 i → k → j 的路径提供更紧的候选值：uᵢⱼ ← min(uᵢⱼ, uᵢₖ + uₖⱼ)。对下界，i 与 j 不能比其中一段减去另一段更近：lᵢⱼ ← max(lᵢⱼ, lᵢₖ − uₖⱼ, lₖⱼ − uᵢₖ)。

上界必须先收紧。lower bound 的两个候选差值使用 upper bound；如果先用宽松的 u 去抬高 l，后面再改变 u，就可能错过本轮本应排除的距离。这个先后关系是算法正确性的一部分，不是 GPU 优化的可选项。

3. Floyd–Warshall：把三重循环拆成 GPU 轮次

Floyd-Warshall algorithm 原本是全对最短路径算法。把它用于 upper smoothing 时，外层固定 k，内层对所有 (i, j) 执行一次 min 更新；lower smoothing 只需替换更新表达式。串行伪代码的关键结构如下：

for k = 0 .. n - 1:
    for every pair (i, j) in parallel:
        U[i, j] = min(U[i, j], U[i, k] + U[k, j])

这段结构适合 GPU，因为一个输出单元只写 U[i, j]，不需要写入其他边。当前 k 的所有 fragment 都可以并行；只有外层的 k 轮必须按顺序推进，下一轮要看到上一轮完整的矩阵。

在 fragment shader 中，(i, j) 由当前像素位置确定，(i, k) 与 (k, j) 的纹理坐标由 rasterizer 或 shader 提供。每个 fragment 读取三个值，计算一个 min 或 max，然后只输出自己的颜色值。这是“读多个、写一个”的数据并行形状。

4. GPU 实现：纹理、动态更新与 ping-pong

ping-pong buffering 解决外层 k 轮的动态依赖。GPU 不能把同一块纹理同时当作 sampler 输入和 render target；因此第 k 轮从 front 纹理读取，把结果渲染到 back 纹理，结束后交换两个句柄。

for (int k = 0; k < n; ++k) {
    bindTexture(frontBounds);
    bindRenderTarget(backBounds);
    setUniform("uMiddle", k);
    drawFullscreenPrimitive();
    swap(frontBounds, backBounds);
}

这里的 drawFullscreenPrimitive 不是在画蛋白质几何，而是在启动一张覆盖距离矩阵的计算网格。每个 fragment 的颜色就是一个新的距离边界；将纹理坐标生成交给 rasterizer，还能减少 shader 指令并让访问模式更规则。

vectorization 把 GPU 四分量运算单元利用起来。连续四条边可以装进一个 RGBA texel：T[i, j] = (D[i, 4j], D[i, 4j+1], D[i, 4j+2], D[i, 4j+3])。这样渲染矩形的宽度约变为原来的四分之一，但 shader 必须把 D[i, k] 广播到正确通道，并处理 n 不是 4 的倍数时的尾部。

5. 动手实验：观察 k 轮怎样收紧矩阵

先预测：把“完成的 k 轮”从 0 拖到 3，哪一类矩阵单元会先变小？打开 RGBA 向量化后，数值会改变，还是只会改变存储和 fetch 估算？

实验中的上界矩阵使用与官方算法相同的 min 更新，下界矩阵使用依赖 upper 的 max 更新；每次推进都按矩阵副本计算一整轮，避免把同一轮的新值提前泄漏给其他单元。绿色单元表示相对初始边界已经收紧，蓝色行列表示当前中间节点的读路径。

6. 三角存储、精度与规模边界

上下界矩阵是对称的，因此只保存上三角和下三角可以节省空间；原章选择把两个三角分别存进纹理，并在几何上只绘制对应区域。这样做会浪费 (i, j) 与 (j, i) 的重复值，却能减少分支，让 texture fetch 和 fragment 区域更规则。对于 GPU，少一次条件判断有时比少一半 texel 更值得，必须在目标设备上测量。

距离平滑对精度敏感。原章的实验使用 32 位 IEEE 浮点纹理，并在不同矩阵规模上比较 CPU 与 GPU 时间；GPU 的优势来自并行性和带宽，而不是改变了 O(n³) 的算法复杂度。规模继续扩大时，纹理边长、GPU 显存和单次 render target 限制会先成为瓶颈，不能把“更大的 n”当作免费收益。

实现检查可以按下面的顺序收敛：

1. 先用 CPU 标量版本验证 upper 和 lower 的更新式与顺序。
2. 再用双缓冲 fragment 版本逐 k 对比每个矩阵单元。
3. 打开三角区域裁剪，确认只读取合法的上三角或下三角。
4. 最后加入 RGBA packing，并测试 n mod 4 的尾部通道和对称性。

小结

• 距离上下界平滑把稀疏 NMR 约束变成可迭代的矩阵更新。
• 三角不等式用中间节点收紧 upper，再用 upper 抬高 lower。
• Floyd–Warshall 的内层 (i, j) 更新具有 GPU 需要的独立写回形状。
• ping-pong 隔离动态轮次的读写，RGBA packing 减少存储宽度与向量 fetch。
• 精度、纹理尺寸、尾部通道和显存限制决定了 GPU 实现的实际边界。