《SIMD instruction considered harmful》SIMD指令被认为是有害的

作者：大卫·帕特森 (David Patterson) 和安德鲁·沃特曼 (Andrew Waterman)，2017 年 9 月 18 日

原文链接：SIMD Instructions Considered Harmful | SIGARCH

在撰写 《RISC-V 手册》的过程中，我们将 RISC-V 向量代码与 SIMD 进行了比较。 我们对 ARM、MIPS 和 x86 的 SIMD 指令扩展的潜在危害（insidiousness)感到震惊。 基于本书的第 8 章，我们决定在此博客中分享这些见解。

就像阿片类药物一样，SIMD 起步时很简单。 架构师将现有的 64 位寄存器和 ALU 划分为许多 8 位、16 位或 32 位的部分，然后对它们进行并行计算。 操作码提供数据宽度和操作。 数据传输只是加载和存储单个 64 位寄存器。 怎么会有人反对这样做呢？

为了加速 SIMD，架构师随后将寄存器的宽度加倍以同时计算更多部分。 因为传统上 ISA 支持向后二进制兼容性，并且操作码指定了数据宽度，所以扩展 SIMD 寄存器也会扩展 SIMD 指令集。 将 SIMD 寄存器的宽度和 SIMD 指令的数量加倍的每个后续步骤都会导致 ISA 走上复杂性不断增加的道路，这些由处理器设计者、编译器编写者和汇编语言程序员承担。 自 1978 年以来，IA-32 指令集已从 80 条指令增长到大约 1400 条指令，这在很大程度上是由 SIMD 推动的。

在我们看来，利用数据级并行性的一种更古老且更优雅的替代方法是向量架构。 向量计算机从主内存中收集（gather)对象并将它们放入长而连续的向量寄存器中。 流水线执行单元在这些向量寄存器上计算非常高效。 然后向量架构将结果从向量寄存器分散(scatter)回主存。

虽然一个简单的向量处理器可能一次执行一个向量元素，但元素操作根据定义是独立的，因此处理器理论上可以同时计算所有元素。 RISC-V 的最宽数据是 64 bit，今天的向量处理器通常每个时钟周期执行两个、四个或八个 64 bit元素。 当向量长度不是每个时钟周期执行的元素数量的整数倍时，硬件会处理边缘/余数情况。 与 SIMD 一样，向量数据宽度是可变的。 每个寄存器具有 N 个 64 位元素的向量处理器还可以计算具有 2N 个 32 bit、4N 个 16 bit和 8N 个 8 bit元素的向量。

即使是简单的 DAXPY 函数（如下）也说明了我们的论点。

  void daxpy(size_t n, double a, const double x[], double y[])
    {
     for (size_t i = 0; i < n; i++) {
       y[i] = a*x[i] + y[i];
     }
    }

图 1. n = 1000 时执行的指令数和总指令数。簿记代码统计包括：复制标量变量 a 以用于 SIMD 执行，当 n 不是 SIMD 寄存器宽度的倍数时要处理的边缘代码，以及 如果 n 为 0，则跳过主循环。图 2-4 和附录中的文本描述了每个 ISA 的已编译 DAXPY 代码。（簿记代码bookkeeping code）

 图 1 总结了 MIPS-32 SIMD 架构 (MSA)、IA-32 AVX2 和 RV32V 程序的 DAXPY 指令数。 （附录中的图 2 至图 4 显示了 RISC-V、MIPS-32 和 IA-32 的编译器输出。）SIMD 计算代码与簿记代码相形见绌。 MIPS-32 MSA 和 IA-32 AVX2 代码的三分之二到四分之三是 SIMD 开销，用于为主 SIMD 循环准备数据或在 n 不是SIMD 寄存器中的浮点数数量的整数倍数时处理边缘元素 。

附录中图2的RV32V代码不需要这样的簿记代码，这减少了静态指令数。 与 SIMD 不同的是，它有一个向量长度寄存器 vl，这使得向量指令可以在 n 为任意值时工作。

然而，SIMD 和向量处理之间最显着的区别不是静态代码大小。 SIMD 指令执行的指令比 RV32V 多 10 到 20 倍，因为每个 SIMD 循环只执行 2 或 4 个元素，而不是向量情况下的 64 个。 额外的指令获取和指令解码意味着更高的能量来执行相同的任务。

相较于DAXPY的标量版本，附录的图 3 和图 4 中的结果显示，尽管主循环的大小保持不变，SIMD 将静态代码的指令和字节大小大致翻了一番， 执行的动态指令数减少了 2 或 4 倍，这具体取决于 SIMD 寄存器的宽度。 然而，RV32V 矢量静态代码大小仅增加了 1.2 倍，而动态指令数则减少了 43 倍。

动态指令数存在很大差异，这在我们看来，是 SIMD 和向量之间第二大差异。 最糟糕的是 ISA和簿记代码 快速膨胀的大小。 像 MIPS-32 和 IA-32 这样遵循向后二进制兼容原则的 ISA ，每次将 SIMD 宽度加倍时，都必须复制为较窄 SIMD 寄存器定义的所有旧 SIMD 指令。 数以百计的 MIPS-32 和 IA-32 指令是在多代 SIMD ISA 中创建的，未来还会有数百条指令，最近的 AVX-512 扩展再次证明了这种情况。 这种数据级并行的强力方法对汇编语言程序员的认知负担一定是惊人的。 我们如何记住 vfmadd213pd 的含义以及何时使用它？

相比之下，RV32V 代码不受向量寄存器大小的影响。 如果向量寄存器大小扩大，RV32V 不仅不会改变，您甚至不必重新编译。 由于处理器提供最大向量长度 mvl 的值，因此如果处理器设计人员将向量寄存器加倍或减半，图 2 中的代码将保持不变。

总之，由于 SIMD ISA 决定于硬件，实现更高的数据级并行性意味着改变指令集和编译器。 有人可能会争辩说 SIMD 违反了将架构与实现隔离的设计原则。 相比之下，向量 ISA 允许处理器设计人员为其应用程序选择数据并行资源，而不会影响程序员或编译器。

我们认为 RV32V 的向量方法与 ARMv7、MIPS-32 和 IA-32 的不断扩展的 SIMD 架构之间在成本-能量-性能、复杂性和易编程性方面的强烈对比是支持 RISC-V的最有说服力的论据之一 。

附录：RV32V、MIPS MSA 和 x86 AVX 的 DAXPY 代码和备注

图 2 显示了 RV32V 中的矢量代码，RV32V是RISC-V 的可选矢量扩展。 RV32V 允许向量寄存器在不使用时被禁用以减少上下文切换的成本，因此第一步是启用两个双精度浮点寄存器。 假设此示例中的最大向量长度 (mvl) 为 64 个元素。

循环中的第一条指令为后面的向量指令设置向量长度。 指令setvl 将mvl 和n 中的较小者写入向量长度寄存器vl 和t0。 如果循环的迭代次数大于 n，则代码处理数据的最快速度是一次处理 64 个值，因此将 vl 设置为 mvl。 如果 n 小于 mvl，那么我们不能读取或写入超出 x 和 y 末尾的内容，因此我们应该在最后一次迭代中只计算最后 n 个元素。

接下来的指令是两个向量加载 (vld)、一个用于计算的乘加运算 (vfmadd)、一个向量存储 (vst)、四个用于地址算术的指令（slli、add、sub、add）、一个循环测试 完成（bnez），并返回（ret）。

向量架构的强大之处在于，这个 10 指令循环的每次迭代都会启动 3×64 = 192 次内存访问和 2×64 = 128 次浮点乘加运算（假设 n 至少为 64）。 这平均每条指令大约有 19 次内存访问和 13 次操作。 SIMD 的这些比率要差一个数量级。

# a0 is n, a1 is pointer to x[0], a2 is pointer to y[0], fa0 is a
  0:  li t0, 2<<25
  4:  vsetdcfg t0             # enable 2 64b Fl.Pt. registers
loop:
  8:  setvl  t0, a0           # vl = t0 = min(mvl, n)
  c:  vld    v0, a1           # load vector x
  10:  slli   t1, t0, 3        # t1 = vl * 8 (in bytes)
  14:  vld    v1, a2           # load vector y
  18:  add    a1, a1, t1       # increment pointer to x by vl*8
  1c:  vfmadd v1, v0, fa0, v1  # v1 += v0 * fa0 (y = a * x + y)
  20:  sub    a0, a0, t0       # n -= vl (t0)
  24:  vst    v1, a2           # store Y
  28:  add    a2, a2, t1       # increment pointer to y by vl*8
  2c:  bnez   a0, loop         # repeat if n != 0
  30:  ret                     # return

图 2. DAXPY 编译成 RV32V。 这个循环适用于 n 的任何值，包括 0； 当n = 0时，setvl将向量操作变成nop。RV32V将数据类型和大小与向量寄存器相关联，从而减少了向量指令的数量。

我们现在展示 SIMD 的代码是如何快速增长的。 图 3 是 MIPS SIMD 架构 (MSA) 版本，图 4 是使用 SSE 和 AVX2 指令的 IA-32。 （ARMv7 NEON 不支持双精度浮点运算，因此无法帮助 DAXPY）。

由于 MSA 寄存器为 128 位宽，因此每个 MSA SIMD 指令可以对两个浮点数进行运算。 与 RV32V 不同，因为没有矢量长度寄存器，MSA 需要额外的簿记指令来检查 n 的问题值。 当 n 为奇数时，有额外的代码来计算单个浮点乘加，因为 MSA 必须对操作数对进行运算。 在不太可能但可能的情况下，当 n 为零时，位置 10 处的分支将跳过主计算循环。 如果它不绕循环分支，位置 18 (splati.d) 的指令将 a 的副本放入 SIMD 寄存器 w2 的两半。 要在 SIMD 中添加标量数据，我们需要将其复制到与 SIMD 寄存器一样宽。

在主循环内部，两条加载指令（ld.d）将 y 和 x 的两个元素读入 SIMD 寄存器 w0 和 w1，三个指令执行地址运算（addiu、addiu、addu），乘加指令（fmadd.d ) 进行计算，然后测试循环终止 (bne)，然后在分支延迟槽中存储 (st.d) 以保存结果。 主循环终止后，代码检查 n 是否为奇数。 如果是，它使用标量指令执行最后的乘加。 最后一条指令 (jr) 返回到调用站点。

Intel 经历了多代 SIMD 扩展，图 4 使用了这些扩展。 SSE 扩展到 128 位 SIMD 导致了 xmm 寄存器和可以使用它们的指令，而扩展到 256 位 SIMD 作为 AVX 的一部分创建了 ymm 寄存器及其指令。

地址 0 到 25 处的第一组指令从内存加载变量，在 256 位 ymm 寄存器中复制 a 的四份，并在进入主循环之前进行测试以确保 n 至少为 4。 （图 4 的标题更详细地解释了如何进行。）

主循环执行 DAXPY 计算的核心。 地址 27 处的 AVX 指令 vmovapd 将 x 的 4 个元素加载到 ymm0 中。 地址 2c 处的 AVX 指令 vfmadd213pd 将 a (ymm2) 的 4 个副本乘以 x (ymm0) 的 4 个元素，加上 y 的 4 个元素（在地址 ecx+edx*8 的内存中），并将 4 个和放入 ymm0。 以下位于地址 32 的 AVX 指令 vmovapd 将 4 个结果存储到 y 中。 接下来的三个指令递增计数器并在需要时重复循环。

与 MIPS MSA 的代码类似，地址 3e 和 57 之间的“边缘”代码处理 n 不是 4 的倍数的情况。

# a0 is n, a2 is pointer to x[0], a3 is pointer to y[0], $w13 is a
  0:  li       a1,-2
  4:  and      a1,a0,a1         # a1 = floor(n/2)*2 (mask bit 0)
  8:  sll      t0,a1,0x3        # t0 = byte address of a1
  c:  addu     v1,a3,t0         # v1 = &y[a1]
 10:  beq      a3,v1,38         # if y==&y[a1] goto Fringe
                                #    (t0==0 so n is 0 | 1)
 14:  move     v0,a2            # (delay slot) v0 = &x[0]
 18:  splati.d $w2,$w13[0]      # w2 = fill SIMD reg. with copies of a
Main Loop:
 1c:  ld.d     $w0,0(a3)        # w0 = 2 elements of y
 20:  addiu    a3,a3,16         # incr. pointer to y by 2 FP numbers
 24:  ld.d     $w1,0(v0)        # w1 = 2 elements of x
 28:  addiu    v0,v0,16         # incr. pointer to x by 2 FP numbers
 2c:  fmadd.d  $w0,$w1,$w2      # w0 = w0 + w1 * w2
 30:  bne      v1,a3,1c         # if (end of y != ptr to y) go to Loop
 34:  st.d     $w0,-16(a3)      # (delay slot) store 2 elts of y
Fringe:
 38:  beq      a1,a0,50         # if (n is even) goto Done
 3c:  addu     a2,a2,t0         # (delay slot) a2 = &x[n-1]
 40:  ldc1     $f1,0(v1)        # f1 = y[n-1]
 44:  ldc1     $f0,0(a2)        # f0 = x[n-1]
 48:  madd.d   $f13,$f1,$f13,$f0# f13 = f1+f0*f13 (muladd if n is odd)
 4c:  sdc1     $f13,0(v1)       # y[n-1] = f13 (store odd result)
Done:
 50:  jr      ra                # return
 54:  nop                       # (delay slot)

图 3. DAXPY 编译成 MIPS32 MSA。 将此代码与图 2 中的 RV32V 代码进行比较时，SIMD 的簿记开销是显而易见的。MIPS MSA 代码的第一部分（地址 0 到 18）复制 SIMD 寄存器中的标量变量 a 并检查以确保 n 位于 在进入主循环之前至少 2。 MIPS MSA 代码的第三部分（地址 38 到 4c）处理 n 不是 2 的倍数时的边缘情况。在 RV32V 中不需要这样的簿记代码，因为向量长度寄存器 vl 和 setvl 指令处理它。 此代码由 gcc 生成，并设置了标志以生成小代码。

# eax is i, n is esi, a is xmm1,
# pointer to x[0] is ebx, pointer to y[0] is ecx
 0: push   esi
 1: push   ebx
 2: mov    esi,[esp+0xc]   # esi = n
 6: mov    ebx,[esp+0x18]  # ebx = x
 a: vmovsd xmm1,[esp+0x10] # xmm1 = a
 10: mov    ecx,[esp+0x1c]  # ecx = y
 14: vmovddup xmm2,xmm1     # xmm2 = {a,a}
 18: mov    eax,esi
 1a: and    eax,0xfffffffc  # eax = floor(n/4)*4
 1d: vinsertf128 ymm2,ymm2,xmm2,0x1 # ymm2 = {a,a,a,a}
 23: je     3e              # if n < 4 goto Fringe
 25: xor    edx,edx         # edx = 0
Main Loop:
 27: vmovapd ymm0,[ebx+edx*8] # load 4 elements of x
 2c: vfmadd213pd ymm0,ymm2,[ecx+edx*8] # 4 mul adds
 32: vmovapd [ecx+edx*8],ymm0 # store into 4 elements of y
 37: add    edx,0x4
 3a: cmp    edx,eax          # compare to n
 3c: jb     27               # repeat loop if < n
Fringe:
 3e: cmp    esi,eax          # any fringe elements?
 40: jbe    59               # if (n mod 4) == 0 go to Done
Fringe Loop:
 42: vmovsd xmm0,[ebx+eax*8] # load element of x
 47: vfmadd213sd xmm0,xmm1,[ecx+eax*8] # 1 mul add
 4d: vmovsd [ecx+eax*8],xmm0 # store into element of y
 52: add    eax,0x1          # increment Fringe count
 55: cmp    esi,eax          # compare Loop and Fringe counts
 57: jne    42 <daxpy+0x42>  # repeat FringeLoop if != 0
Done:
 59: pop    ebx              # function epilogue
 5a: pop    esi
 5b: ret

图 4. DAXPY 编译成 IA-32 SSE 和 AVX2。 地址 a 处的 SSE 指令 vmovsd 将 a 加载到 128 位 xmm1 寄存器的一半。 地址 14 处的 SSE 指令 vmovddup 将 a 复制到 xmm1 的两半中，用于以后的 SIMD 计算。 地址 1d 处的 AVX 指令 vinsertf128 从 xmm1 中的 a 的两个副本开始，在 ymm2 中制作了四个 a 副本。 地址 42 到 4d 处的三个 AVX 指令（vmovsd、vfmadd213sd、vmovsd）在 mod(n,4) ≄ 0 时处理。它们一次执行一个元素的 DAXPY 计算，循环重复直到函数正好执行 n 个乘加操作。 同样，这样的代码对于 RV32V 是不必要的，因为向量长度寄存器 vl 和 setvl 指令使循环适用于任何 n 值。

作者简介：大卫·帕特森 (David Patterson) 是加州大学伯克利分校计算机科学研究生院的教授，也是谷歌的杰出工程师。 Andrew Waterman 是 SiFive 的总工程师，拥有加州大学伯克利分校的计算机科学博士学位。

免责声明：这些帖子是由个人贡献者撰写的，目的是为了社区的利益在今日计算机体系结构博客上分享他们的想法。 本博客中表达的任何观点或意见均为个人观点，仅属于博客作者，不代表 ACM SIGARCH 或其上级组织 ACM 的观点。