Skip to main content

2.1.1 GSSW SIMD

要理解 SIMD Smith-Waterman 如何打包单元格并并行计算，最好的办法是用一个具体的小例子，一步一步走一遍。我们用一个字长 W=4W=4 的 SIMD 寄存器，参考序列长度 88，来看同一行的 88 个单元格怎么塞进 44 个 SIMD 通道里并行计算。

1. 先明确问题：为什么不能直接并行​

Smith-Waterman 计算单元格 (i,j)(i,j) 时依赖三个方向：

H(i,j)=max⁡(0, H(i−1,j−1)+S(i,j), H(i−1,j)+gap, H(i,j−1)+gap)H(i,j) = \max\big(0,\ H(i-1,j-1)+S(i,j),\ H(i-1,j)+gap,\ H(i,j-1)+gap\big)

其中：

  • H(i−1,j−1)H(i-1,j-1)：左上，来自上一行
  • H(i−1,j)H(i-1,j)：上边，来自上一行
  • H(i,j−1)H(i,j-1)：左边，来自同一行

最后这个“左边”依赖，就是同一行内不能直接并行的原因。因为要算 j=2j=2 必须先算 j=1j=1，要算 j=3j=3 必须先算 j=2j=2，串行下去。

2. Farrar 的破解方法：先猜后修正​

Farrar 的洞察是：在真实序列比对中，绝大多数情况下，左边单元格 H(i,j−1)H(i,j-1) 加上空位罚分后，不会成为最大值。也就是说，最大值通常来自上一行的两个方向。

于是可以：

  1. 先忽略左边依赖，只根据上一行计算一个“临时分数”。
  2. 算完后检查：左边依赖是否真的影响结果？如果影响了，再修正。

这样，同一行的多个单元格就可以先并行算临时值，再统一修正。

3. 具体例子​

假设：

  • 查询序列只有一个字符：Q=AQ = A
  • 参考序列：R=A C G T A C G TR = A\ C\ G\ T\ A\ C\ G\ T，长度 88
  • 匹配得分 +2+2，错配得分 −1-1，空位罚分 −2-2
  • SIMD 字长 W=4W = 4

我们现在要计算第 i=1i=1 行（查询字符 AA 对应的一行），共 88 个单元格。

3.1 条带化布局：把一行拆进 SIMD 通道​

SIMD 寄存器有 44 个通道，记为通道 0,1,2,30,1,2,3。我们把参考序列的列按“模 44”分到不同通道：

通道负责的列
0列 1,51, 5
1列 2,62, 6
2列 3,73, 7
3列 4,84, 8

也就是说，一个 SIMD 寄存器一次处理 44 个单元格：

  • 第一组 SIMD 操作：列 1,2,3,41,2,3,4 分别放在通道 0,1,2,30,1,2,3
  • 第二组 SIMD 操作：列 5,6,7,85,6,7,8 分别放在通道 0,1,2,30,1,2,3

论文图 4a 里说的“相同颜色的单元格打包进一个 SIMD 字”，就是这个意思：同一次 SIMD 操作里的单元格颜色相同，它们列号相差 WW，彼此之间没有直接依赖。

3.2 第一组：并行计算列 1 到 4​

上一行 i−1=0i-1=0 全是 00，因为没有比对上任何字符。

匹配得分向量（查询 AA 对参考 A,C,G,TA,C,G,T）：

S=[+2, −1, −1, −1]S = [+2,\ -1,\ -1,\ -1]

忽略左边依赖，只根据上一行计算临时值：

Htemp=max⁡(0, Hprev+S)H_{\text{temp}} = \max(0,\ H_{\text{prev}} + S)

HprevH_{\text{prev}} 是上一行对应位置的向量，这里全是 00。所以：

Htemp=max⁡(0, [0,0,0,0]+[+2,−1,−1,−1])=[2, 0, 0, 0]H_{\text{temp}} = \max(0,\ [0,0,0,0] + [+2,-1,-1,-1]) = [2,\ 0,\ 0,\ 0]

一个 SIMD 寄存器就同时算出了列 11 到列 44 的临时分数。

3.3 检查并修正左边依赖​

现在检查左边依赖是否会改变结果。左边依赖来自同一行的前一个单元格，加上空位罚分 −2-2。

  • 列 11：没有左边，保持 22。
  • 列 22：左边是列 11 的分数 22，2+(−2)=02 + (-2) = 0。临时值是 00，0≥00 \ge 0，不需要修正。
  • 列 33：左边是列 22 的分数 00，0+(−2)=−20 + (-2) = -2。临时值是 00，不需要修正。
  • 列 44：左边是列 33 的分数 00，0+(−2)=−20 + (-2) = -2。临时值是 00，不需要修正。

所以第一组结果就是：

H(1,1..4)=[2, 0, 0, 0]H(1,1..4) = [2,\ 0,\ 0,\ 0]

3.4 第二组：并行计算列 5 到 8​

同样方法。查询 AA 对参考 A,C,G,TA,C,G,T 的匹配得分向量还是：

S=[+2, −1, −1, −1]S = [+2,\ -1,\ -1,\ -1]

上一行仍然是 00，所以：

Htemp=max⁡(0, [0,0,0,0]+[+2,−1,−1,−1])=[2, 0, 0, 0]H_{\text{temp}} = \max(0,\ [0,0,0,0] + [+2,-1,-1,-1]) = [2,\ 0,\ 0,\ 0]

再检查左边依赖，同样不需要修正。所以：

H(1,5..8)=[2, 0, 0, 0]H(1,5..8) = [2,\ 0,\ 0,\ 0]

4. 整体效果​

没有 SIMD 时，这一行 88 个单元格要串行算 88 次。

用 SIMD 后：

  • 第一组 SIMD 操作：一次算列 11 到 44
  • 第二组 SIMD 操作：一次算列 55 到 88

总共只需要 22 次 SIMD 计算，加上少量的修正检查。

这就是“同一行多个单元格并行计算”的核心：

把列号相差 WW 的单元格放进同一个 SIMD 寄存器的不同通道，让它们同时算临时值；左边依赖通过事后检查和修正来处理。

5. 图版本 GSSW 的额外变化​

上面的例子是线性 Smith-Waterman 的 SIMD 并行。GSSW 在此基础上还要处理图：

  • 节点内部：像上面一样，一行内多个单元格打包进 SIMD 并行算。
  • 节点第一行：不能只看上一行，而要看所有父节点。需要先把所有父节点对应位置的分数合并，取最大值，再作为当前节点第一行的初始值。这一步是标量或间接访问，不能直接 SIMD 并行。
  • 所以 GSSW 的计算模式是：密集 SIMD 区域（节点内部）和间接图访问（节点边界）交替出现。

论文里用红色箭头表示“父节点到当前节点第一行”的依赖，说的就是这个节点初始化步骤。

6. 一句话总结​

SIMD Smith-Waterman 把同一行中列号相差 WW 的单元格打包进一个 SIMD 寄存器的不同通道，先忽略左边依赖并行算临时值，再检查并修正。这样一行 88 个单元格只需要 22 次 SIMD 操作，而不是 88 次串行计算。GSSW 在这个基础上，把参考从线性序列换成图：节点内部继续 SIMD 并行，节点边界则要合并多个父节点的信息。