12.8. Core2上相同的例子
归功于更好的μop融合与更强大的执行单元,相同的循环在Core2运行得更高效。例子12.6c与12.6d都是Core2上运行的好候选。在例子12.6c中的j1 L1指令应该改为jb L1,以启用32位模式中的宏操作融合。这个改变是可能的,因为在例子12.6c中,循环计数器不能是负数。
现在,我们将分析对一个Core2处理器,例子12.6d中的循环的资源使用。循环中每条指令使用的资源列在下面的表12.2中。
指令 | 指令长度 | μop融合域 | 每个执行端口的μop | ||||||||
|
|
| 端口0 | 端口1 | 端口5 | 端口2 | 端口3 | 端口4 | |||
movapd xmm1, [esi+eax] | 5 | 1 |
|
|
| 1 |
|
| |||
mulpd xmm1, xmm2 | 4 | 1 | 1 |
|
|
|
|
| |||
addpd xmm1, [edi+eax] | 5 | 1 |
| 1 |
| 1 |
|
| |||
movapd [edi+eax], xmm1 | 5 | 1 |
|
|
|
| 1 | 1 | |||
add eax, 16 | 3 | 1 | x | x | x |
|
|
| |||
js L1 | 2 | 1 |
|
| 1 |
|
|
| |||
Total | 24 | 6 | 1.33 | 1.33 | 1.33 | 2 | 1 | 1 | |||
表12.2. Core2上DAXPY循环的μop总数
这里预解码不是一个问题,因为循环的大小小于64字节。没有必要对齐循环,因为其大小小于50字节。如果循环超过64字节,那么我们将注意到预解码器在一个时钟周期里仅可以处理前3条指令,因为它不能处理超过16字节。对解码器没有限制,因为所有的指令每条仅产生一个μop。解码器每时钟周期可以处理4条指令,因此解码时间将是每迭代6/5 = 1.5时钟周期。
端口0、1与2每个有一个不去别处的μop。另外,指令ADD EAX, 16可去往这三个端口的其中一个。如果这些指令在端口0、1与2间均匀分布,那么这些端口平均每迭代忙碌1.33时钟周期。
端口3每迭代接收2个μop,用于读X[i]及Y[i]。因此,我们可以断定端口3上的内存读是DAXPY循环的瓶颈,执行时间是每迭代2时钟周期。
例子12.6c比例子12.6d多一条指令,我们已经分析过它了。这个额外的指令是CMP EAX, ECX,它去往端口0、1与2的其中一个。这将增加这些端口的压力到1.67,这仍然小于端口3上的2个μop。在32位模式中,额外的μop可以通过宏操作融合消除。浮点加法与乘法单元都有吞吐率1。因此,这些单元不是DAXPY循环中的瓶颈。
我们需要考虑循环中寄存器读暂停的可能性。例子12.6c循环里有4个读,但不修改的寄存器。它们是寄存器ESI,EDI,XMM2与ECX。这会导致一次寄存器读暂停,如果在4个连续μop内读所有这些寄存器。不过读这些寄存器的μop被分得如此开,在4个类型μop中,不可能读超过其中3个寄存器。例子12.6d少读一个寄存器,因为ECX不再使用。因此,在这些循环里,寄存器读暂停不太可能造成麻烦。
结论是,Core2处理器运行DAXPY循环的时间只需PM的一半。在PM上使优化特别困难的μop冲突问题,在Core2设计中被消除了。
12.9. Sandy Bridge上相同的例子
向量寄存器中元素数量,在具有AVX指令集的处理器上,使用YMM寄存器加倍了, 如下例所示:
; Example 12.6f. Loop of DAXPY for processors with AVX
.data
SignBit DD 0, 80000000H ; qword with sign bit set
n = 100 ; Define constant n (divisible by 4)
.code
vbroadcastsd ymm0, [SignBit] ; preload sign bit x 4
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea rsi, [X + 8 * n] ; Point to end of array X (aligned)
lea rdi, [Y + 8 * n] ; Point to end of array Y (aligned)
vbroadcastsd ymm2, [DA] ; Load DA x 4
vxorpd ymm2, ymm2, ymm0 ; Change sign
align 32
L1: vmulpd ymm1, ymm2, [rsi+rax] ; X[i]*(-DA)
vaddpd ymm1, ymm1, [rdi+rax] ; Y[i]-X[i]*DA
vmovapd [rdi+rax], ymm1 ; Store (vmovupd if not aligned by 32)
add rax, 32 ; Add size of four elements to index
jl L1 ; Loop back
vzeroupper ; If subsequent code is non-VEX
现在我们将分析,对Sandy Bridge处理器,例子12.6f中循环的资源使用。循环中的每条指令的资源列出在下面表12.3中。
指令 | 指令长度 | μop融合域 | 每个执行端口的μop | ||||||||
|
|
| 端口0 | 端口1 | 端口5 | 端口2 | 端口3 | 端口4 | |||
movapd xmm1, [esi+eax] | 5 | 1 | 1 |
|
| 1+ |
|
| |||
mulpd xmm1, xmm2 | 5 | 1 |
| 1 |
|
| 1+ |
| |||
addpd xmm1, [edi+eax] | 5 | 1 |
|
|
| 1 |
| 1+ | |||
movapd [edi+eax], xmm1 | 4 | ½ |
|
| ½ |
|
|
| |||
add eax, 16 | 2 | ½ |
|
| ½ |
|
|
| |||
Total | 21 | 4 | 1 | 1 | 1 | 2 | 1+ | 1+ | |||
表12.3. Sandy Bridge上DAXPY循环的μop总数(例子12.6f)
指令长度不重要,如果我们假定循环在μop缓存中。ADD与JL指令很可能融合在一起,一起仅使用一个μop。有两个256位读操作,每个连续两个时钟周期使用一个读端口,在表中表示为1+。使用两个读端口(端口2与3),在两个时钟周期里,将有两次256位读的吞吐率。其中一个读端口,对在第二个时钟周期里的写,进行地址计算。写端口(端口4)被256位写占据两个时钟周期。限制因素将是读与写操作,将两个读端口以及写端口用到最大限度。端口0、1与5用了一半的容量。没有循环携带依赖链。因此,预期吞吐率是两个时钟周期,一次迭代的四个计算。
不幸的是,由于大量的内存操作,缓存库冲突的可能性很高。很难预测哪些操作将同时读或写,因为同时进行许多迭代。我在Sandy Bridge上的测试显示,每迭代2时钟周期的理论时间,从来没有得到过。最好的结果是每迭代2.2时钟周期,当RSI与RDI间的距离是4k字节倍数时。在其他情形下,时间上每得到2.6时钟周期或更多。
12.10. 使用FMA4相同的例子
AMD Bulldozer是第一个提供融合乘加指令的x86处理器。这些指令可以显著提高吞吐率。
; Example 12.6g. Loop of DAXPY for processors with FMA4
.code
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea rsi, [X + 8 * n] ; Point to end of array X (aligned)
lea rdi, [Y + 8 * n] ; Point to end of array Y (aligned)
vmovddup xmm2, [DA] ; Load DA x 2
align 32
L1: vmovapd xmm1, [rsi+rax] ; X[i]
vfnmaddpd xmm1, xmm1, xmm2, [rdi+rax] ; Y[i]-X[i]*DA
vmovapd [rdi+rax], xmm1 ; Store
add rax, 16 ; Add size of two elements to index
jl L1 ; Loop back
这个循环不受执行单元的限制,而是内存访问。不幸的是,在Bulldozer上缓存库冲突非常频繁。这将消耗时间从理论上的每迭代2时钟周期,增加到几乎4。在Bulldozer上,使用更大的YMM寄存器没有优势。另一方面,YMM双倍加长指令的解码有些低效。
12.11. 使用FMA3相同的例子
Intel Haswell处理器提供3操作数形式、称为FMA3的融合乘加指令。AMD Piledriver同时支持FMA3与FMA4。
; Example 12.6h. Loop of DAXPY with FMA3, using xmm registers
.code
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea rsi, [X + 8 * n] ; Point to end of array X (aligned)
lea rdi, [Y + 8 * n] ; Point to end of array Y (aligned)
vmovddup xmm2, [DA] ; Load DA x 2
align 32
L1: vmovapd xmm1, [rdi+rax] ; Y[i]
vfnmadd231pd xmm1, xmm2, [rsi+rax] ; Y[i]-X[i]*DA
vmovapd [rdi+rax], xmm1 ; Store
add rax, 16 ; Add size of two elements to index
jl L1 ; Loop back
在Intel处理器上,这个循环每迭代需要1 ~ 2时钟周期。使用YMM寄存器,吞吐率几乎加倍:
; Example 12.6i. Loop of DAXPY with FMA3, using ymm registers
.code
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea rsi, [X + 8 * n] ; Point to end of array X
lea rdi, [Y + 8 * n] ; Point to end of array Y
vbroadcastsd xmm2, [DA] ; Load DA x 4
align 32
L1: vmovupd ymm1, [rsi+rax] ; X[i]
vfnmadd231pd ymm1, ymm2, [rdi+rax] ; Y[i]-X[i]*DA
vmovupd [rdi+rax], ymm1 ; Store
add rax, 32 ; Add size of two elements to index
jl L1 ; Loop back
例子12.6i无需假定数组对齐到32。因此,VMOVAPD可以改为VMOVUPD。在Haswell与Broadwell上,这个循环的执行时间测得为2时钟周期,在Skylake上是1.5时钟周期。瓶颈不是这里执行的单元,而是缓存效果以及可能分支指令。
12.12. 循环展开
一个重复n次的循环可以被一个重复n / r次、每次执行r次计算的循环所替代,其中r是展开因子。n最好能被r整除。
循环展开可用于以下目的:
- 减少循环开销。每次计算的循环开销除展开因子r。仅在循环开销对计算时间有显著贡献时,这才有用。如果其他某个瓶颈限制了执行速度,没有理由展开一个循环。例如,上面例子12.6e中的循环就不能从展开获益。
- 向量化。为了使用有r元素的向量寄存器,循环必须展开r或r的倍数。例子12.6e中的循环展开2次,以使用有两个双精度值的向量。如果我们使用了单精度数,那么可以展开循环4次,使用有4个元素的向量。
- 改进分支预测。通过展开循环,使重复计数n / r不超过特定CPU可以预测的最大重复计数,可以改进循环退出分支的预测。
- 改进缓存。如果循环遭遇许多数据缓存不命中或缓存竞争,以对特定处理器最优的方式安排内存读写,可能是有利的。在现代处理器上,这很少需要。细节参考微处理器提供商提供的优化手册。
- 消除整数除法。如果循环包含一个表达式,其中循环计数器i除一个整数r,或者计算i模r,通过展开循环r次,可以避免这个整数除法。
- 消除循环内分支。如果带有周期r重复模式的循环内有一个分支或switch语句,可以通过展开该循环r次消除之。例如,如果一个if-else分支每二次去往同一个方向,这个分支可以通过展开2次消除。
- 打破循环携带依赖链。在某些情形里,可以通过使用多个累加器,打破循环携带依赖链。展开因子r等于累加器数。参考第59页,例子9.3b。
- 减少对归纳变量的依赖。如果从前面的迭代值计算一个归纳变量的时延很高、成为一个瓶颈,那么有可能可以通过展开r次,从序列中后r个位置处的值计算归纳变量的每个值。
- 完全展开。在r = n时,完全展开一个循环。这完全消除了循环开销。每个是循环计数器函数的表达式可被常量替代。可以消除仅依赖循环计数器的每个分支。参考第111页例子。
在重复计数n不能被展开因子r整除时,循环指令存在一个问题。将存在n modulo r的余数个、不能在循环里完成的额外计算。这些额外计算必须在主循环前或后完成。
在重复计数不能被展开因子整除时,正确获取额外计算,可能相当困难;如果使用如例子12.6d与e中的负索引,它变得特别棘手。下面的例子再次展示DAXPY算法,这次使用单精度并展开4次。在这个例子中,n是不一定被4整除的变量。数组x与y必须对齐到16。(为了清晰起见,忽略特定于PM处理器的优化)。
; Example 12.7. Unrolled Loop of DAXPY, single precision.
.data
align 16
SignBitS DD 80000000H ; dword with sign bit set
.code
mov eax, n ; Number of calculations, n
sub eax, 4 ; n - 4
lea esi, [X + eax*4] ; Point to X[n-4]
lea edi, [Y + eax*4] ; Point to Y[n-4]
movss xmm2, DA ; Load DA
xorps xmm2, SignBitS ; Change sign
shufps xmm2, xmm2, 0 ; Get -DA into all four dwords of xmm2
neg eax ; -(n-4)
jg L2 ; Skip main loop if n < 4
L1: ; Main loop rolled out by 4
movaps xmm1, [esi+eax*4] ; Load 4 values from X
mulps xmm1, xmm2 ; Multiply with -DA
addps xmm1, [edi+eax*4] ; Add 4 values from Y
movaps [edi+eax*4],xmm1 ; Store 4 results in Y
add eax, 4 ; i += 4
jle L1 ; Loop as long as <= 0
L2: ; Check for remaining calculations
sub eax, 4 ; = -remainder
jns L4 ; Skip extra loop if remainder = 0
L3: ; Extra loop for up to 3 remaining calculations
movss xmm1, [esi+eax*4+16] ; Load 1 value from X
mulss xmm1, xmm2 ; Multiply with -DA
addss xmm1, [edi+eax*4+16] ; Add 1 value from Y
movss [edi+eax*4+16],xmm1 ; Store 1 result in Y
add eax, 1 ; i += 1
js L3 ; Loop as long as negative
L4:
展开在数组上进行计算的循环的另一个解决方案是,将数组扩展到具有r-1个不使用的空间,将重复计数n取整到展开因子r最接近的倍数。这消除了计算余数(n mod r)以及余下计算的额外循环的需要。不使用的数组元素必须初始化为0或其他有效的浮点值,以避免次正规数、NAN、溢出、下溢,或其他任何会减慢浮点计算的情形。如果数组是整数类型,只需避免除0。
仅在有原因这样做,且可以获得速度的显著提升时,才进行循环展开。应该避免过分的循环展开。循环展开的坏处有:
- 代码变得更大,在代码缓存中占据更多空间。这会导致代码缓存不命中,使展开得不偿失。注意,在孤立测试循环时,不能检测到代码缓存不命中。
- 许多处理器有一个循环缓冲来提升非常小循环的速度。循环缓冲被限制为20 ~ 40条指令或64字节代码,依赖于处理器。如果超出循环缓冲的大小,循环展开很可能降低性能,(处理器特定细节在手册3《Intel,AMD与VIA CPU微架构》中提供)。
- 某些处理器有一个有限大小的μop缓存。这个μop缓存非常有价值,应该经济地使用之。展开的循环在μop缓存中占据更多的空间。注意,在孤立测试循环时,不能检测到μop缓存不命中。
- 在n不能被r整除时,在展开循环外部需要进行额外计算,使得代码更复杂、笨拙且增加了分支的数量。
- 展开的循环可能需要更多的寄存器,如用于多个累加器。
12.13. 使用掩码寄存器的向量循环(AVX512)
AVX512指令集提供了受掩蔽的操作。在循环计数不能被向量大小整除时,这个特性对掩蔽超出的向量元素是有用的。
; Example 12.8. Vectorized DAXPY loop using mask registers
.data
align 64
countdown dd 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0
sixteen dd 16
.code
lea rsi, [X] ; point to beginning of X
lea rdi, [Y] ; point to beginning of Y
mov edx, n ; number of elements, n
vmovd xmm0, edx
vpbroadcastd zmm0, xmm0 ; broadcast n
vpaddd zmm0, zmm0, zmmword [countdown] ; counts = n+countdown
vpbroadcastd zmm1, dword [sixteen] ; broadcast 16
vbroadcastss zmm2, dword [DA] ; broadcast DA
xor eax, eax ; loop counter i = 0
L1: ; Main loop rolled out by 16
vpcmpnltud k1, zmm0, zmm1 ; mask k1 = (counts >= 16)
vpsubd zmm0, zmm0, zmm1 ; counts -= 16
vmovups zmm3 {k1}{z}, zmmword [rdi+rax*4] ; load Y[i]
vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi+rax*4] ; Y[i]-DA*X[I]
vmovups zmmword [rdi+rax*4] {k1}, zmm3 ; store Y[i]
add rax, 16 ; i += 16
cmp eax, rdx ; while i < n
jb L1
这里,在k1寄存器中的掩码,对所有有效的元素置1比特,对超出n的元素是0。保存指令被k1掩蔽,避免保存超出n个Y元素。使用相同的掩码,掩蔽读与计算指令同样也是好的。这里,必须使用0掩蔽读指令,以避免对向量寄存器先前值的一个假依赖。掩蔽计算节省了能耗,避免可能的次正规值等导致的异常及性能损失。向一条512位向量指令添加掩码,没有代价。
产生掩码的指令不需要与计算相同的向量大小与粒度。例如,如果计算有双精度(64位粒度),我们可以从具有32位粒度的256位寄存器,甚至16位粒度的128位寄存器,产生所需的8位掩码,如果n保证小于216且支持AVX512BW。
要性能最好,数组X与Y最好对齐到64。
如果向量加法的执行单元是一个瓶颈,制作掩码可以使用整数指令而不是向量指令。这在下面的例子中展示:
; Example 12.9. Vectorized DAXPY, mask calculated with integer instr.
; Note: This version works only for n < 256!
.code
lea rsi, [X] ; point to beginning of X
lea rdi, [Y] ; point to beginning of Y
vbroadcastss zmm2, dword[DA] ; broadcast DA
mov edx, n ; number of elements, n
mov ecx, -1 ; fill with all 1's
xor eax, eax ; loop counter i = 0
L1: ; Main loop rolled out by 16
bzhi ecx, ecx, edx ; zero bit positions > edx
kmovw k1, ecx ; copy 16 bits to mask register
vmovups zmm3 {k1}{z}, zmmword [rdi+rax*4] ; load Y[i]
vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi+rax*4] ; Y[i]-DA*X[I]
vmovups zmmword [rdi+rax*4] {k1}, zmm3 ; store Y[i]
add rax, 16 ; i += 16
sub edx, 16 ; count down n
ja L1 ; repeat if n positive
例子12.9倒数在edx中余下的元素数,如果edx < 16,在最后一次迭代中,使用bzhi清除ecx里余下的比特。指令bzhi属于BMI2指令集,所有已知支持AVX512的处理器都支持它。注意,这仅对n < 256可工作,因为bzhi仅读edx的低8位。
如果迭代计数很高且计算掩码的指令减慢了循环的执行,主循环不带掩码,在主循环后进行要求掩码的操作,会更好:
; Example 12.10. Vectorized DAXPY, masking only after main loop
.code
lea rsi, [X] ; point to beginning of X
lea rdi, [Y] ; point to beginning of Y
vbroadcastss zmm2, dword[DA] ; broadcast DA
mov edx, n ; number of elements, n
and edx, -16 ; round down n to multiple of 16
lea rsi, [rsi+rdx*4] ; point to the end
lea rdi, [rdi+rdx*4] ; point to the end
neg rdx ; use negative index from the end
L1: ; Main loop rolled out by 16
vmovups zmm3, zmmword [rdi+rdx*4] ; load Y[i]
vfnmadd231ps zmm3, zmm2, zword [rsi+rdx*4] ; Y[i]-DA*X[I]
vmovups zmmword [rdi+rdx*4] {k1}, zmm3 ; store Y[i]
add rdx, 16 ; next 16 elements
jnz L1 ; count up to zero
mov edx, n ; calculate remaining elements
and edx, 15 ; n modulo 16
jz L2 ; optionally skip if zero
mov eax, -1 ; all 1's
bzhi eax, eax, edx ; zero bit positions > edx
kmovw k1, eax ; copy 16 bits to mask register
; rsi and rdi are pointing to where the main loop ended.
; Now do the last 0-15 elements
vmovups zmm3 {k1}{z}, zmmword [rdi] ; load Y[i]
vfnmadd231ps zmm3 {k1}, zmm2, zword [rsi] ; Y[i]-DA*X[I]
vmovups zmmword [rdi] {k1}, zmm3 ; store Y[i]
L2:
例子12.10包含循环体两次。第一个主体在L1循环重复n/16次。第二个主体在n不能被16整除时,计算余下的0 ~ 15个元素。如果这个分支预测不好,没有必要使用jz L2跳转到最后的主体。
例子12.10在循环体中仅有5条指令。如果循环计数很高,这将是最快的解决方案。如果循环主体很大,我们不希望包含它两次,或者循环计数小,或者瓶颈是内存访问或循环携带依赖链,例子12.8与12.9是有用的。例子12.9仅用在n < 256时。
12.14. 优化缓存
在访问非缓存内存的循环中,内存访问很可能需要最多时间。如果可能且连续访问,数据应该连续保存,如第11章77页所述。
在循环中数组访问数不应该超出微处理器中读写缓冲数量。减少数据流的一个方式是将多个数组合并为一个结构体数组,使多个数据流被交织为单个流。
现代微处理器有先进的预取机制。这些机制可以检测数据访问模式中的规律,比如以特定步长访问数据。建议尽量减少不同数据流数,尽可能保持固定步长,以利用这样的预取机制。如果数据访问模式足够规律,自动数据预取通常比显式数据预取工作得更好。
在数据访问模式不够规律,自动预取机制无法预测时,使用prefetch指令显式预取数据可能是必须的。对一个以不规律方式访问数据的程序,找出最优的预取策略,通常需要大量的实验。
如果系统有多个CPU核,可以将数据预取放入一个独立的线程。Intel C++有这样的一个特性。
使用一个2指数步长的数据访问很可能导致缓存行竞争。这可以通过改变步长或循环分块来避免。细节参考手册《优化C++软件》中优化内存访问的章节。
对写不太可能很快访问的非缓存内存,非临时写指令是有用的。可以使用向量指令尽量减少非临时写指令,
12.15. 并行
提升CPU密集代码性能的最重要方式是并行。并行的主要方法有:
- 提升CPU进行乱序执行的可能性。这通过打破长依赖链(参考第59页),以及在不同执行单元或执行端口间均匀分配μop(参考第90页)来完成。
- 使用向量指令。参考第13章,113页。
- 使用多线程。参考第14章,142页。
可以使用多个累加器,打破循环携带依赖链,如第59页所述。如果CPU没有别的事情可做,最优累加器数是依赖链中最关键指令的时延除该指令的吞吐率倒数。例如,在AMD处理器上,浮点加法的时延是4时钟周期,吞吐率倒数是1。这意味着最优累加器数是4。下面的例子12.1展示看使用4个浮点寄存器作为累加器,执行加法的循环。
// Example 12.11a, Loop-carried dependency chain
// (Same as example 9.3a page 65)
double list[100], sum = 0.;
for (int i = 0; i < 100; i++) sum += list[i];
使用4个浮点寄存器作为累加器的实现,看起来像这样:
; Example 12.11b, Four floating point accumulators
lea esi, list ; Pointer to list
fld qword ptr [esi] ; accum1 = list[0]
fld qword ptr [esi+8] ; accum2 = list[1]
fld qword ptr [esi+16] ; accum3 = list[2]
fld qword ptr [esi+24] ; accum4 = list[3]
fxch st(3) ; Get accum1 to top
add esi, 800 ; Point to end of list
mov eax, 32-800 ; Index to list[4] from end of list
L1:
fadd qword ptr [esi+eax] ; Add list[i]
fxch st(1) ; Swap accumulators
fadd qword ptr [esi+eax+8] ; Add list[i+1]
fxch st(2) ; Swap accumulators
fadd qword ptr [esi+eax+16] ; Add list[i+2]
fxch st(3) ; Swap accumulators
add eax, 24 ; i += 3
js L1 ; Loop
faddp st(1), st(0) ; Add two accumulators together
fxch st(1) ; Swap accumulators
faddp st(2), st(0) ; Add the two other accumulators
faddp st(1), st(0) ; Add these sums
fstp qword ptr [sum] ; Store the result
在例子12.11中,从list向4个累加器载入前4个值。
在循环中执行的加法数恰好被展开因子3整除。使用浮点寄存器作为累加器有趣的地方是,累加器数等于展开因子加1。这是交换累加器的fxch指令的使用方式的结果。你必须密切注意浮点寄存器栈上的每个累加器,验证在循环的每次迭代后,4个累加器旋转一个位置,使每个累加器每4个加法得到使用,虽然循环仅展开3次。
例子12.11b中的循环每次加法需要1时钟周期,这是浮点加法器的最大吞吐率。浮点加法4时钟周期的时延,由4个累加器处理。指令fxch有零时延,因为在大多数Intel、AMD与VIA处理器(除了Intel Atom)上,它们被翻译为寄存器重命名。
在具有SSE2指令集的处理器上,通过使用XMM寄存器,而不是浮点栈寄存器,可以避免fxch指令,如例子12.11c所示。在AMD上,浮点向量加法的时延是4,吞吐率倒数是2,因此,最优累加器数是2个向量寄存器。
; Example 12.11c, Two XMM vector accumulators
lea esi, list ; list must be aligned by 16
movapd xmm0, [esi] ; list[0], list[1]
movapd xmm1, [esi+16] ; list[2], list[3]
add esi, 800 ; Point to end of list
mov eax, 32-800 ; Index to list[4] from end of list
L1:
addpd xmm0, [esi+eax] ; Add list[i], list[i+1]
addpd xmm1, [esi+eax+16] ; Add list[i+2], list[i+3]
add eax, 32 ; i += 4
js L1 ; Loop
addpd xmm0, xmm1 ; Add the two accumulators together
movhlps xmm1, xmm0 ; There is no movhlpd instruction
addsd xmm0, xmm1 ; Add the two vector elements
movsd [sum], xmm0 ; Store the result
在AMD处理器上,例子12.11b与12.11c一样快,因为两者都受到浮点加法器吞吐率的限制。
在Intel Core2处理器上,例子12.11c比12.11b更快,因为这个处理器有可以在一次操作中处理整个向量的128位浮点加法器、
12.16. 分析依赖性
循环可能有几个互锁的依赖链。这样复杂的情形要求仔细的分析。
下一个例子是泰勒展开。你可能知道,许多函数可通过这个形式的泰勒多项式逼近
fx≈ i=0ncixi
每个指数xi通过之前指数xi-1乘以x,方便地计算。系数ci保存在一个表中。
; Example 12.12a. Taylor expansion
.data
align 16
one dq 1.0 ; 1.0
x dq ? ; x
coeff dq c0, c1, c2, ... ; Taylor coefficients
coeff_end label qword ; end of coeff. list
.code
movsd xmm2, [x] ; xmm2 = x
movsd xmm1, [one] ; xmm1 = x^i
xorps xmm0, xmm0 ; xmm0 = sum. init. to 0
mov eax, offset coeff ; point to c[i]
L1: movsd xmm3, [eax] ; c[i]
mulsd xmm3, xmm1 ; c[i] * x^i
mulsd xmm1, xmm2 ; x^(i+1)
addsd xmm0, xmm3 ; sum += c[i] * x^i
add eax, 8 ; point to c[i+1]
cmp eax, offset coeff_end ; stop at end of list
jb L1 ; loop
(如果你的汇编器混淆了movsd指令与同名的字符串指令,把它编码为DB 0F2H / movups)。
现在关注分析。系数列表很短,我们可以预期它待在缓存中。
为了检查时延是否重要,我们必须查看代码中的依赖。这些依赖显示在图12.1中。
图12.1. 例子12.12泰勒展开中的依赖
有两个连续的依赖链,一个用于计算xi,一个用于计算和。在Intel处理器上,指令mulsd的时延比addsd长。因此,垂直乘法链比加法链更关键。加法必须等待cixi,这个值比xi晚一个乘法时延,且比前面的加法晚。如果没有别的限制性能,那么可以预期这个代码每迭代需要一个乘法时延,这是4 ~ 7时钟周期,取决于处理器。
吞吐率不是限制因素,因为在具有128位执行单元的处理器上,乘法单元的吞吐率是每时钟周期一次乘法。
测量例子12.12a的时间,我们发现循环比乘法多用了至少一个时钟周期。对此解释如下:两个乘法都要等待前面等待在xmm1中的值xi-1。因此,两个乘法在同一时间就绪。我们希望图12.1中阶梯里的垂直乘法先开始,因为它是最关键依赖链的部分。但是微处理器看不到交换这两个乘法的理由,因此图12.1中的水平乘法先开始。垂直乘法被推迟1时钟周期,这是浮点乘法单元的吞吐率倒数。可以通过把图12.1中的水平乘法放在垂直乘法后,可以解决这个问题:
; Example 12.12b. Taylor expansion
movsd xmm2, [x] ; xmm2 = x
movsd xmm1, [one] ; xmm1 = x^i
xorps xmm0, xmm0 ; xmm0 = sum. initialize to 0
mov eax, offset coeff ; point to c[i]
L1: movapd xmm3, xmm1 ; copy x^i
mulsd xmm1, xmm2 ; x^(i+1) (vertical multipl.)
mulsd xmm3, [eax] ; c[i]*x^i (horizontal multipl.)
add eax, 8 ; point to c[i+1]
cmp eax, offset coeff_end ; stop at end of list
addsd xmm0, xmm3 ; sum += c[i] * x^i
jb L1 ; loop
这里,我们需要额外的指令来拷贝xi(movapd xmm3, xmm1),因此我们可以在水平乘法前进行垂直乘法。这使得循环每迭代的执行快1时钟周期。
我们仍然仅使用了xmm寄存器低半部。同时执行两个乘法,一次进行泰勒展开的两步,增益更大。技巧是从xi-2乘上x2计算xi的每个值:
; Example 12.12c. Taylor expansion, double steps
movsd xmm4, [x] ; x
mulsd xmm4, xmm4 ; x^2
movlhps xmm4, xmm4 ; x^2, x^2
movapd xmm1, xmmword ptr [one] ; xmm1(L)=1.0, xmm1(H)=x
xorps xmm0, xmm0 ; xmm0 = sum. init. to 0
mov eax, offset coeff ; point to c[i]
L1: movapd xmm3, xmm1 ; copy x^i, x^(i+1)
mulpd xmm1, xmm4 ; x^(i+2), x^(i+3)
mulpd xmm3, [eax] ; c[i]*x^i, c[i+1]*x^(i+1)
addpd xmm0, xmm3 ; sum += c[i] * x^i
add eax, 16 ; point to c[i+2]
cmp eax, offset coeff_end ; stop at end of list
jb L1 ; loop
haddpd xmm0, xmm0 ; join the two parts of sum
现在速度加倍了。循环每迭代仍然需要乘法时延,但迭代次数减半。乘法时延是4 ~ 5时钟周期,取决于处理器。每迭代执行两个向量乘法,我们仅使用了最大乘法器吞吐率的40或50%。这意味着提升并行度,有更大的增益。可以一次进行4步泰勒展开,xi-4乘上x4计算xi每个值来做到:
; Example 12.12d. Taylor expansion, quadruple steps
movsd xmm4, [x] ; x
mulsd xmm4, xmm4 ; x^2
movlhps xmm4, xmm4 ; x^2, x^2
movapd xmm2, xmmword ptr [one] ; xmm2(L)=1.0, xmm2(H)=x
movapd xmm1, xmm2 ; xmm1 = 1, x
mulpd xmm2, xmm4 ; xmm2 = x^2, x^3
mulpd xmm4, xmm4 ; xmm4 = x^4, x^4
xorps xmm5, xmm5 ; xmm5 = sum. init. to 0
xorps xmm6, xmm6 ; xmm6 = sum. init. to 0
mov eax, offset coeff ; point to c[i]
L1: movapd xmm3, xmm1 ; copy x^i, x^(i+1)
movapd xmm0, xmm2 ; copy x^(i+2), x^(i+3)
mulpd xmm1, xmm4 ; x^(i+4), x^(i+5)
mulpd xmm2, xmm4 ; x^(i+6), x^(i+7)
mulpd xmm3, xmmword ptr [eax] ; term(i), term(i+1)
mulpd xmm0, xmmword ptr [eax+16] ; term(i+2), term(i+3)
addpd xmm5, xmm3 ; add to sum
addpd xmm6, xmm0 ; add to sum
add eax, 32 ; point to c[i+2]
cmp eax, offset coeff_end ; stop at end of list
jb L1 ; loop
addpd xmm5, xmm6 ; join two accumulators
haddpd xmm5, xmm5 ; final sum in xmm5
在例子12.12d中,我们使用两个寄存器xmm1与xmm2来保存x的连续指数,以将图12.1中的垂直依赖链分解为4个并行的链。每个xi指数从xi-4计算。类似的,我们对和使用两个累加器,xmm5与xmm6,以将加法链分解为4条并行的链。在循环后,四个和被加起来。在Core2上的一个测量显示,对例子12.12d,每迭代5.04时钟周期,或每泰勒项1.26时钟周期。这接近于由乘法时延确定的、每迭代5时钟周期的期望。这个解决方案使用了乘法器最大吞吐率的80%,这是一个非常令人满意的结果。进一步展开几乎没有增益。向量加法与其他指令不会增加执行时间,因为它们不需要使用与乘法相同的执行单元。在这个情形里,指令获取与解码不是瓶颈。
在带有256位YMM寄存器、支持AVX指令集的处理器上,吞吐率可以进一步提高:
; Example 12.12e. Taylor expansion using AVX
vmovddup xmm2, [x] ; x, x
vmulpd xmm5, xmm2, xmm2 ; x^2, x^2
vmulpd xmm4, xmm5, xmm5 ; x^4, x^4
vmovsd xmm2, [one] ; 1
vmovhpd xmm2, xmm2, [x] ; 1, x
vmulpd xmm0, xmm2, xmm5 ; x^2, x^3
vinsertf128 ymm2, ymm2, xmm0, 1 ; 1, x, x^2, x^3
vinsertf128 ymm4, ymm4, xmm4, 1 ; x^4, x^4, x^4, x^4
vmulpd ymm3, ymm2, ymm4 ; x^4, x^5, x^6, x^7
vmulpd ymm4, ymm4, ymm4 ; x^8, x^8, x^8, x^8
vxorps xmm0, xmm0 ; sum0 = 0
vxorps xmm1, xmm1 ; sum1 = 0
lea rax, [coeff_end] ; point to end of coeff[i]
mov rcx, -n ; counter
jmp L2 ; jump into loop
align 32 ; loop
L1: vmulpd ymm2, ymm2, ymm4 ; multiply powers of x by x^8
vmulpd ymm3, ymm3, ymm4 ; multiply powers of x by x^8
L2: vmulpd ymm5, ymm2, [rax+rcx] ; first four terms
vaddpd ymm0, ymm0, ymm5 ; sum0
vmulpd ymm5, ymm3, [rax+rcx+32] ; next four terms
vaddpd ymm1, ymm1, ymm5 ; sum1
add rcx, 64
jl L1
; make total sum
vaddpd ymm0, ymm0, ymm1 ; join two accumulators
vextractf128 xmm5, ymm0, 1 ; get high part
vaddpd xmm0, xmm0, xmm5 ; add high and low part
vhaddpd xmm0, xmm0, xmm0 ; final sum in xmm0
vzeroupper ; if subsequent code is non-VEX
例子12.12e每迭代计算8项。在Sandy Bridge上,测得执行时间为每迭代5.5时钟周期,接近由5时钟周期乘法时延确定的理论最小值。在Skylake上,测得执行时间大致相同,尽管在这个处理器上乘法时延仅是4时钟周期。
当然,在例子12.12d与12.12e中,花在循环外的时间比前面的例子多,但性能增益仍然可观,即使对于不多的迭代。最后的一小点别扭是,在循环前将第一个系数放入累加器,从x1开始,而不是x0。可以通过最后加第一项得到更好的精度,但代价是一个额外的加法时延。
可以在循环中使用融合乘加指令(FMA3或FMA4),进一步改进例子12.12e的代码,如果支持这些指令。在Skylake上测得,在使用FMA指令时,执行时间从5.4减为4.5时钟周期。
总之,垂直乘法应该在水平乘法前完成。如果循环迭代次数很高,最佳累加器寄存器数是使用执行单元最大可能吞吐率的值。如果循环迭代次数很低,循环前的设置时间与循环后的加法都要紧,那么最佳累加器寄存器数可能要小一些。
一个泰勒展开通常在项变成负数时结束。不过,总是包括最坏情形的最大项数,以避免检查停止条件所需的浮点比较,可能更好。重复次数常量也将防止循环退出的误预测。计算比所需更多的项的代价,在每项可能使用不到1时钟周期的优化例子中实际上相当小。
记住将mxcsr寄存器置为“Flush to zero模式”,以避免在执行比所需更多次迭代时,可能下溢导致的性能损失。
-
- 没有乱序执行处理器上的循环
P1与PMMX处理器没有乱序执行能力。当然,这些处理器今天已经过时了,但这里描述的原则可能适用于更小的嵌入式处理器。
我选择了一个简单的例子——一个从数组读整数,改变每个整数的符号,把结果保存到另一个数组的过程。这个过程的C++语言代码应该是:
// Example 12.13a. Loop to change sign
void ChangeSign (int * A, int * B, int N) {
for (int i = 0; i < N; i++) B[i] = -A[i];
}
一个对P1优化的汇编实现看起来像这样:
; Example 12.13b
mov esi, [A]
mov eax, [N]
mov edi, [B]
xor ecx, ecx
lea esi, [esi+4*eax] ; point to end of array a
sub ecx, eax ; -n
lea edi, [edi+4*eax] ; point to end of array b
jz short L3
xor ebx, ebx ; start first calculation
mov eax, [esi+4*ecx]
inc ecx
jz short L2
L1: sub ebx, eax ; u
mov eax, [esi+4*ecx] ; v (pairs)
mov [edi+4*ecx-4], ebx ; u
inc ecx ; v (pairs)
mov ebx, 0 ; u
jnz L1 ; v (pairs)
L2: sub ebx, eax ; end last calculation
mov [edi+4*ecx-4], ebx
L3:
这里,迭代是重叠的,以提升成对的机会。在保存第一个值之前,开始读第二个值。指令mov ebx, 0放在inc ecx与jnz L1之间,不是改进成对,而是第一个指令对使用ecx作为地址索引,避免在递增ecx后导致的AGI暂停。
带有浮点操作的循环有点不同,因为浮点指令是流水线化重叠的,而不是成对。考虑第90页例子12.6的DAXPY循环。对P1优化的一个解决方案如下:
; Example 12.14. DAXPY optimized for P1 and PMMX
mov eax, [n] ; number of elements
mov esi, [X] ; pointer to X
mov edi, [Y] ; pointer to Y
xor ecx, ecx
lea esi, [esi+8*eax] ; point to end of X
sub ecx, eax ; -n
lea edi, [edi+8*eax] ; point to end of Y
jz short L3 ; test for n = 0
fld qword ptr [DA] ; start first calc.
fmul qword ptr [esi+8*ecx] ; DA * X[0]
jmp short L2 ; jump into loop
L1: fld qword ptr [DA]
fmul qword ptr [esi+8*ecx] ; DA * X[i]
fxch ; get old result
fstp qword ptr [edi+8*ecx-8] ; store Y[i]
L2: fsubr qword ptr [edi+8*ecx] ; subtract from Y[i]
inc ecx ; increment index
jnz L1 ; loop
fstp qword ptr [edi+8*ecx-8] ; store last result
L3:
这里我们将循环计数器用作数组索引,从负数数到0。
每个操作在之前一个结束前开始,以提升计算重叠。具有乱序执行的处理器自动这样做,但对没有乱序能力的处理器,我们必须显式进行这个重叠。
这里,浮点操作的交错完美工作:之前结果的FSTP填充到FMUL与FSUBR间2时钟的暂停。循环开销与下一个操作的前两条指令填充到FSUBR与FSTP之间3时钟暂停。在递增索引后的第一个时钟周期中读不依赖于索引计数器参数,避免了地址生成互锁(AGI)暂停。
这个解决方案每迭代需要6时钟周期,比别处发布的展开解决方案要好。
12.18. 宏循环
如果一个循环的重复计数小且为常量,完全展开该循环是可能的。这样的好处是,仅依赖于循环计数器的计算可以在汇编时刻完成,而不是在执行时刻。当然坏处是,如果重复次数很高,它占据更多的代码缓存。
MASM语法提供了一个完成支持元编程的强大宏语言,它相当有用。其他汇编器有类似的能力,但语法不同。例如,如果我们需要一个平方数列表,C++代码看起来像这样:
// Example 12.15a. Loop to make list of squares
int squares[10];
for (int i = 0; i < 10; i++) squares[i] = i*i;
在MASM语言中可以通过一个宏循环产生相同的列表:
; Example 12.15b. Macro loop to produce data
.DATA
squares LABEL DWORD ; label at start of array
I = 0 ; temporary counter
REPT 10 ; repeat 10 times
DD I * I ; define one array element
I = I + 1 ; increment counter
ENDM ; end of REPT loop
这里,I是一个预处理变量。I循环在汇编时刻运行,但不在执行时刻。变量I与语句I = I + 1不会进入最终的代码,因此无需时间执行。实际上,例子12.15b不产生可执行代码,仅数据。宏预处理器将上面的代码翻译为:
; Example 12.15c. Resuls of macro loop expansion
squares LABEL DWORD ; label at start of array
DD 0
DD 1
DD 4
DD 9
DD 16
DD 25
DD 36
DD 49
DD 64
DD 81
宏循环对产生代码也是有用的。下一个例子计算xn,其中x是一个浮点值,n是一个正整数。通过反复取x平方,与n里对应二进制数字的因子相乘,这可以最高效地执行。这个算法的C++代码:
// Example 12.16a. Calculate pow(x,n) where n is a positive integer
double x, xp, power;
unsigned int n, i;
xp = x; power = 1.0;
for (i = n; i != 0; i >>= 1) {
if (i & 1) power *= xp;
xp *= xp;
}
如果n在汇编时刻已知,那么可以使用性能的宏循环实现这个指数函数:
; Example 12.16b.
; This macro will raise a double-precision float in X
; to the power of N, where N is a positive integer constant.
; The result is returned in Y. X and Y must be two different
; XMM registers. X is not preserved.
; (Only for processors with SSE2)
INTPOWER MACRO X, Y, N
LOCAL I, YUSED ; define local identifiers
I = N ; I used for shifting N
YUSED = 0 ; remember if Y contains valid data
REPT 32 ; maximum repeat count is 32
IF I AND 1 ; test bit 0
IF YUSED ; If Y already contains data
mulsd Y, X ; multiply Y with a power of X
ELSE ; If this is first time Y is used:
movsd Y, X ; copy data to Y
YUSED = 1 ; remember that Y now contains data
ENDIF ; end of IF YUSED
ENDIF ; end of IF I AND 1
I = I SHR 1 ; shift right I one place
IF I EQ 0 ; stop when I = 0
EXITM ; exit REPT 32 loop prematurely
ENDIF ; end of IF I EQ 0
mulsd X, X ; square X
ENDM ; end of REPT 32 loop
ENDM ; end of INTPOWER macro definition
这个宏产生了完成这项工作所需的最少指令。在最终代码中没有循环开销、prolog与epilog。最重要的,没有分支。所有的分支都被宏预处理器解决了。要计算xmm0的12次方,写:
; Example 12.16c. Macro invocation
INTPOWER xmm0, xmm1, 12
这将被展开为:
; Example 12.16d. Result of macro expansion
mulsd xmm0, xmm0 ; x^2
mulsd xmm0, xmm0 ; x^4
movsd xmm1, xmm0 ; save x^4
mulsd xmm0, xmm0 ; x^8
mulsd xmm1, xmm0 ; x^4 * x^8 = x^12
这甚至比没有展开的优化汇编循环的指令还要少。当mulpd替换mulsd,movapd替换movsd时,这个宏还可以在向量上工作。