分类： ARM

从一个复数点积算法看NEON的汇编优化

摘要：本文通过一个真实案例（4096点双精度浮点复数点积算法），描述了使用 Zynq-7000 NEON进行算法优化的过程以及一些关键技巧，相对于使用编译器对C代码做优化，性能提升了大约4.8倍。 本文介绍的内容对需要用到NEON实现高性能计算的开发者非常有帮助。

开发工具：Xilinx SDK 2013.4
开发板： Xilinx ZC702或者ZC706

一般来说，使用NEON优化算法有以下几种方式：
> 用现成的开源库，例如Ne10，FFTW等等
>  使用编译器的auto vectorization功能，通过配置合适的编译选项让编译器来生成NEON代码
>  使用NEON Intrinsic
>  使用NEON汇编

前三种方式简单易用，但是受限于编译器的优化能力，往往性能无法达到极致。在需要使用NEON进行高性能计算的时候，就需要用手工汇编来进行优化 了。用汇编语言写程序看起来很让人望而生畏，但实际上在一个算法里面需要高度优化的核心部分并不会太大，只要静下心来，掌握一些书写NEON汇编的技巧， 完全可以轻松的完成核心汇编程序。

下面我们来看一个例子：求两个矢量（vector）的点积（dot product），矢量的长度为4096，矢量的数据为复数（complex），复数的实部和虚部都是双精度浮点类型。，每个矢量占用内存64KB。

对Zynq-7000芯片来说，数据可以保存在片内高速SRAM（OCM）中，也可以保存在DDR3 memory中。我们首先以数据保存在OCM中为例，最后会讨论数据保存在DDR3中的情况。

1. 使用编译器的auto vectorization进行优化：

对这个算法来说，C语言的标准实现大致是这样的。

double a[N*2] ;
double b[N*2] ;
double result_golden[2] = {0.0, 0.0};

int i;
result_golden[0]=0.0;
result_golden[1]=0.0;
for(i = 0; i < 4096; i++)
{
    result_golden[0] += a[2 * i] * b[2 * i] - a[2 * i + 1] * b[2 * i + 1];
    result_golden[1] += a[2 * i + 1] * b[2 * i] + a[2 * i] * b[2 * i + 1];
}

我们使用以下编译选项来让编译器针对NEON进行优化。

-O2-mcpu=cortex-a9 -mfpu=neon -ftree-vectorize -mvectorize-with-neon-quad -mfloat-abi=softfp -ffast-math

我们将这个算法运行10次，得到平均的执行时间为289.692us。缺省情况下，CPU的主频为667MHz，换算成CPU cycles，执行时时间为193128 cycles。将这个数字除以矢量长度4096，我们就得到每个双精度复数乘加需要的平均时间为47.14 CPU cycles.

通过上面的源码，我们可以看出每个复数乘加涉及到3次双精度浮点数的乘加（MLA）运算和1次双精度浮点数的乘减（MLS）运算。查询 ARM Cortex™-A9 NEON™ Media Processing Engine Technical Reference Manual文档，我们注意到对双精度浮点数，在理想情况下，MLA/MLS指令需要2个CPU cycles，即理论上执行一次复数乘加需要的时间为8 CPU cycles。

即使我们把数据在memory和NEON registers之间传输的时间也考虑进去，每个复数乘加的实测执行时间和理论执行时间之间的差距也是相当大的。这样，我们可以初步得出结论：在这个算 法上，编译器的优化能力离性能极限还有很大差距，通过适当的汇编优化，还有很大的性能提升空间。

2. 如何实现最优的指令序列

我们先来做个实验，测试一下以下函数的执行时间。在这个函数的循环体里面，有4条vmla.f64指令。测试完成后，我们从后往前每次注释掉一条vmla.f64指令，然后测试实际的执行时间，直到最后只剩1条vmla.f64指令为止。

.align 4
.global neon_instr_seq_0
.arm
neon_instr_seq_0:
LDR r0, =16384
neon_instr_seq_0_loop:
vmla.f64 d8, d0, d4
vmla.f64 d9, d1, d5
vmla.f64 d10, d2, d6
vmla.f64 d11, d3, d7
SUBS r0, r0, #1
BNE neon_instr_seq_0_loop
bx lr

这样我们就得到了以下的表格

函数的循环体里面vmla.f64指令的数目 每次循环需要的CPU cycles 平均每条NEON指令需要的CPU cycles
4 8 2
3 6 2
2 6 3
1 6 6

从这个表格，我们可以看出：因为CPU内部微架构的原因，对一个NEON指令的目的寄存器写入后，要等待一段时间再执行下一次写入操作，否则会导致 pipeline stall导致性能下降。对双精度浮点数，相邻的三个NEON指令的目的寄存器必须是不同的。值得注意的是，对其他数据类型，例如单精度浮点数或者32- bit整数，结论会略有差异。

一般来说，在memory和register之间用VLDM/VSTM交换数据会有比较高的效率。对NEON来说，共有32个D寄存器，这样我们对 每个矢量每次可以加载4个复数（即8个双精度浮点数），其中D0-D7保存矢量1的数据，D8-D15保存矢量2的数据，其他寄存器用于保存中间结果。按 照这个思路，比较直接的汇编写法是：

.align 4
.global neon_instr_seq_20
.arm
neon_instr_seq_20:
LDR r0, =16384
neon_instr_seq_20_loop:
vmla.f64 d16, d0, d8
vmla.f64 d17, d0, d9
vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8

vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10

vmla.f64 d20, d4, d12
vmla.f64 d21, d4, d13
vmls.f64 d20, d5, d13
vmla.f64 d21, d5, d12

vmla.f64 d22, d6, d14
vmla.f64 d23, d6, d15
vmls.f64 d22, d7, d15
vmla.f64 d23, d7, d14
SUBS r0, r0, #1
BNE neon_instr_seq_20_loop
bx lr

每次循环使用到了16条FPU指令，理论上需要32 CPU cycles，实测每次循环需要40 CPU cycles。这也印证了一开始我们通过测试得出的结论，所以我们要通过重新排布指令的顺序来取得更好的性能。

重新排布后指令顺序如下所示。实测每次循环需要32 CPU cycles，这段代码达到了最优化。

.align 4
.global neon_instr_seq_21
.arm
neon_instr_seq_21:
LDR r0, =16384
neon_instr_seq_21_loop:
vmla.f64 d16, d0, d8
vmla.f64 d17, d0, d9
vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11
vmla.f64 d20, d4, d12
vmla.f64 d21, d4, d13
vmla.f64 d22, d6, d14
vmla.f64 d23, d6, d15

vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10
vmls.f64 d20, d5, d13
vmla.f64 d21, d5, d12
vmls.f64 d22, d7, d15
vmla.f64 d23, d7, d14
SUBS r0, r0, #1
BNE neon_instr_seq_21_loop
bx lr

对上面的代码，使用到的目的寄存器为8个，结合最开始的测试结论，我们只需要4个目的寄存器就够了，所以可以把上面的代码改写成下面的形式。实测每次循环需要32 CPU cycles，性能仍然是最优。

.align 4
.global neon_instr_seq_22
.arm
neon_instr_seq_22:
LDR r0, =16384
neon_instr_seq_22_loop:
vmla.f64 d16, d0, d8
vmla.f64 d17, d0, d9
vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11
vmla.f64 d16, d4, d12
vmla.f64 d17, d4, d13
vmla.f64 d18, d6, d14
vmla.f64 d19, d6, d15

vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10
vmls.f64 d16, d5, d13
vmla.f64 d17, d5, d12
vmls.f64 d18, d7, d15
vmla.f64 d19, d7, d14
SUBS r0, r0, #1
BNE neon_instr_seq_22_loop
bx lr

3. 如何优化使用VLDM

比较直观的使用VLDM指令后的汇编如下所示。实测每次循环需要45 CPU cycles.

.align 4
.global neon_instr_seq_23
.arm
neon_instr_seq_23:
LDR r3, =16384
neon_instr_seq_23_loop:
vldm r0, {d0-d7}
vldm r1, {d8-d15}

vmla.f64 d16, d0, d8
vmla.f64 d17, d0, d9
vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11
vmla.f64 d16, d4, d12
vmla.f64 d17, d4, d13
vmla.f64 d18, d6, d14
vmla.f64 d19, d6, d15

vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10
vmls.f64 d16, d5, d13
vmla.f64 d17, d5, d12
vmls.f64 d18, d7, d15
vmla.f64 d19, d7, d14
SUBS r3, r3, #1
BNE neon_instr_seq_23_loop
bx lr

在这里我们可以做一个有趣的实验，我们可以注释掉任意一条VLDM指令，然后测试每次循环的执行时间。实测的结果是33 CPU cycles。这是一件很有意思的事情。在Cortex-A9内核里面，NEON/FPU是一个独立的处理单元，Load/Store Unit也是一个独立处理单元，两个单元可以并行执行。Load/Store Unit可以在加载到第一个数据的时候立即把数据forward给NEON/FPU，这样在只有一条VLDM指令时，只需要引入1个CPU cycles的延迟了。这就提示我们通过合理的打散VLDM指令并和NEON/FPU指令交织（interleave），可以提升指令执行的并行度，并继 续提升软件的性能。

到了这个层面，就需要对CPU的微架构（micro architecture）有比较深入的了解了。不过对于软件工程师，没有必要在这方面花太多的时间，只需要注意一些基本的原则，反复尝试几种可能的组合就好了：
> 在VLDR/VLDM和使用到相关寄存器的数据处理指令之间留出足够的时间，避免pipeline stall
>  VLDM指令能够快速的把数据forward给寄存器，所以VLDM和数据处理指令之间不要有其他load/store指令

按照上面的原则，我们对代码的顺序重新调整，得到了下面的汇编代码。实测每次循环的执行时间从45 CPU cycles缩短为39 CPU cycles。

值得注意的是在这里我们假定数据已经在L1 cache里面了。在实际中，数据要从memory加载到L1 cache中，然后再传到NEON register里面，这个过程非常复杂。通过上述方法找到的最优指令序列在真实环境中有可能未必是最优的，这时就需要测试几个在这种简单情况下最优和次 优的指令序列在真实环境中的性能，从而找出真实环境中性能最优的代码。

.align 4
.global neon_instr_seq_35
.arm
neon_instr_seq_35:
LDR r3, =16384
neon_instr_seq_35_loop:
vldr d0, [r0, #0]
vldm r1, {d8-d15}

vmla.f64 d16, d0, d8
vldr d2, [r0, #16]
vmla.f64 d17, d0, d9

vldr d4, [r0, #32]
vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11

vldr d6, [r0, #48]
vmla.f64 d16, d4, d12
vmla.f64 d17, d4, d13

vldr d1, [r0, #8]
vmla.f64 d18, d6, d14
vmla.f64 d19, d6, d15

vldr d3, [r0, #24]
vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8

vldr d5, [r0, #40]
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10

vldr d7, [r0, #56]
vmls.f64 d16, d5, d13
vmla.f64 d17, d5, d12

vmls.f64 d18, d7, d15
vmla.f64 d19, d7, d14
SUBS r3, r3, #1
BNE neon_instr_seq_35_loop
bx lr

4. 使用PLD指令优化性能

上面的讨论都假设数据已经存在于L1 cache中了。在真实的环境中未必会这样，这时就需要用PLD指令提前把数据加载到L1 cache中。算法计算量越大，PLD指令对性能的提升也越明显。

使用PLD指令需要注意：
PLD指令只能加载一个cache line，对Cortex-A9来说一个cache line为32 bytes。
Cortex-A9硬件限制最多有4条PLD指令并行执行，过多的插入PLD指令不会有额外的性能提升。

增加了PLD指令后的完整的汇编代码如下所示：

.align 4
.global complex_dot_product_neon_var4
.arm
complex_dot_product_neon_var4:
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@ r0: address of source vector a
@ r1: address of source vector b
@ r2: address of destination complex
@ r3: vector length
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

vmov.i32 d16, #0
vmov.i32 d17, #0
vmov.i32 d18, #0
vmov.i32 d19, #0

loop_var4:

vldr d0, [r0, #0]
vldm r1!, {d8-d15}

pld [r0, #64]
pld [r0, #96]
pld [r1, #64]
pld [r1, #96]

vmla.f64 d16, d0, d8
vldr d2, [r0, #16]
vmla.f64 d17, d0, d9

vldr d4, [r0, #32]
vmla.f64 d18, d2, d10
vmla.f64 d19, d2, d11

vldr d6, [r0, #48]
vmla.f64 d16, d4, d12
vmla.f64 d17, d4, d13

vldr d1, [r0, #8]
vmla.f64 d18, d6, d14
vmla.f64 d19, d6, d15

vldr d3, [r0, #24]
vmls.f64 d16, d1, d9
vmla.f64 d17, d1, d8

vldr d5, [r0, #40]
vmls.f64 d18, d3, d11
vmla.f64 d19, d3, d10

vldr d7, [r0, #56]
vmls.f64 d16, d5, d13
vmla.f64 d17, d5, d12

vmls.f64 d18, d7, d15
vmla.f64 d19, d7, d14

add r0, r0, #64
SUBS r3, r3, #4
BNE loop_var4

vadd.f64 d16, d16, d18
vadd.f64 d17, d17, d19

vstm r2!, {d16-d17}
bx lr

5. 测试结果及分析：

数据在OCM和DDR3上的执行结果为：

测试项 在OCM上的执行时间 在DDR3上的执行时间
Compiler Auto-vectorization 289.677 272.526
asm_1 175.563 396.615
asm_1 with PLD 70.203 180.765
asm_2 147.861 264.484
asm_2 with PLD 60.522 170.733

>  Compiler Auto-vectorization为本文最开始的C代码通过GCC编译器优化后的测试
>  asm_1为不考虑第3步VLDM/VLDR和NEON数据处理指令交织的汇编代码
>  asm_2为考虑到第3步VLDM/VLDR和NEON数据处理指令交织后的汇编代码

从测试结果来看，我们可以注意到几点：
>  充分优化后的代码（asm_2 with PLD）在OCM上运行的结果非常接近NEON的最优性能。NEON数据处理指令需要32 CPU cycles，根据Zynq-7000 TRM， OCM的latency为23 cycles。这时PLD指令可以在处理当前加载的数据的同时，提前把下一批要处理的数据加载到L1 cache中。

>  对代码（asm_2 with PLD）来说，数据在DDR3上的性能要差于在OCM上的性能。这个也比较好理解，因为DDR3的latency要远大于NEON数据处理指令执行时间 （32 CPU cycles），在整个运算过程中，PLD指令并不能有效的提前把下一批待处理的数据加载到L1 cache中。或者说，只有部分数据能够被提前加载到L1 cache中。

>  NEON数据处理指令执行时间越长，就越容易抵消掉内存系统延迟带来的影响。

>  数据在OCM上时，系统性能要优于数据在DDR3上。不过这时还要额外考虑数据在OCM上搬进搬出的开销。当然，如果在系统设计时就考虑到这一点，由硬件直接使用 OCM，系统性能就能够得到明显的提升。

6. 小结
通过这个真实案例，我们演示了一些NEON/FPU汇编开发的技巧。通过活学活用这些技巧，我们就能充分发挥Zynq-7000的计算能力，实现一个高度优化的高性能嵌入式系统。

7.5.1. The PLD instruction -- DDI0434C_cortex_a5_mpcore_trm

参考链接

乘积累加运算(Multiply Accumulate, MAC)

乘积累加运算（英语：Multiply Accumulate, MAC）是在数字信号处理器或一些微处理器中的特殊运算。实现此运算操作的硬件电路单元，被称为“乘数累加器”。这种运算的操作，是将乘法的乘积结果和累加器 A 的值相加，再存入累加器：

若没有使用 MAC 指令，上述的程序可能需要二个指令，但 MAC 指令可以使用一个指令完成。而许多运算（例如卷积运算、点积运算、矩阵运算、数字滤波器运算、乃至多项式的求值运算）都可以分解为数个 MAC 指令，因此可以提高上述运算的效率。

MAC指令的输入及输出的数据类型可以是整数、定点数或是浮点数。若处理浮点数时，会有两次的数值修约（Rounding），这在很多典型的DSP上很常见。若一条MAC指令在处理浮点数时只有一次的数值修约，则这种指令称为“融合乘加运算”/“积和熔加运算”（fused multiply-add, FMA）或“熔合乘法累积运算”（fused multiply–accumulate, FMAC）。

积和熔加运算

融合乘加运算的操作和乘积累加的基本一样，对于浮点数的操作也是一条指令完成。但不同的是，非融合乘加的乘积累加运算，处理浮点数时，会先完成b×c的乘积，将其结果数值修约到N个比特，然后才将修约后的结果与寄存器a的数值相加，再把结果修约到N个比特；融合乘加则是先完成a+b×c的操作，获得最终的完整结果后方才修约到N个比特。由于减少了数值修约次数，这种操作可以提高运算结果的精度，以及提高运算效率和速率。

积和融加运算可以显著提升像是这些运算的性能和精度：

积和融加运算通常被依靠用来获取更精确的运算结果。然而，Kahan指出，如果不加思索地使用这种运算操作，在某些情况下可能会带来问题。[1]像是平方差公式x2 − y2，它等价于 ((x×x) − y×y)，若果x与y已知数值，使用积和融加运算来求结果，哪怕x = y时，因为在进行首次乘法操作时无视低位的有效比特，可能会使运算结果出错，如果是多步运算，第一步就出错则会连累后续的运算结果接连出错，比如前述的平方差求值后，再取结果的平方根，那么这个结果也会出错。

参考链接