12. 循环
一个CPU密集程序的关键热点几乎总是一个循环。现代计算机的时钟频率非常高,即使最耗时的指令、缓存不命中及低效的异常,都可以不到1毫秒完成。由低效代码导致的时延仅在重复数百万次时,才是显著的。如此高的重复次数很可能仅在一系列嵌套循环的最内层出现。可以做的就是,在本章中讨论的,改进循环的性能。
12.1. 使循环开销最小
循环开销是跳回循环开始以及确定何时退出循环所需的指令。优化这些指令是适用于许多情形的一个相当通用的技术。不过,如果其他某个瓶颈限制了速度,优化循环开销是不需要的。循环中可能出现瓶颈的讨论,参考第89页。
C++中一个典型的循环看起来像这样:
// Example 12.1a. Typical for-loop in C++
for (int i = 0; i < n; i++) {
// (loop body)
}
没有优化,其汇编实现看起来像这样:
; Example 12.1b. For-loop, not optimized
mov ecx, n ; Load n
xor eax, eax ; i = 0
LoopTop:
cmp eax, ecx ; i < n
jge LoopEnd ; Exit when i >= n
; (loop body) ; Loop body goes here
add eax, 1 ; i++
jmp LoopTop ; Jump back
LoopEnd:
使用inc指令对循环计数器加1,可能是不明智的。Inc指令有部分写flags寄存器的问题,这使得它在Intel P4处理器上,比add指令要低效,在其他处理器可能导致假依赖。
在例子12.1b中循环附带的最重要的问题是,有两条跳转指令。通过把分支指令放在末尾,可以消除一个来自循环的跳转:
; Example 12.1c. For-loop with branch in the end
mov ecx, n ; Load n
test ecx, ecx ; Test n
jng LoopEnd ; Skip if n <= 0
xor eax, eax ; i = 0
LoopTop:
; (loop body) ; Loop body goes here
add eax, 1 ; i++
cmp eax, ecx ; i < n
jl LoopTop ; Loop back if i < n
LoopEnd:
现在通过把循环退出分支放在最后,我们去掉了循环里的无条件跳转。在循环前,我们必须放入一个额外的检查来对循环运行0次提供保护。没有这个检查,在n = 0时,循环将运行一次。
把循环退出分支放在末尾的方法甚至可用于退出条件在中间的复杂循环结构、考虑一个退出条件在中间的C++循环:
// Example 12.2a. C++ loop with exit in the middle
int i = 0;
while (true) {
FuncA(); // Upper loop body
if (++i >= n) break; // Exit condition here
FuncB(); // Lower loop body
}
通过重新组织循环,使退出出现在末尾,入口出现在中间,在汇编中可以这样实现:
; Example 12.2b. Assembly loop with entry in the middle
xor eax, eax ; i = 0
jmp LoopEntry ; Jump into middle of loop
LoopTop:
call FuncB ; Lower loop body comes first
LoopEntry:
call FuncA ; Upper loop body comes last
add eax, 1
cmp eax, n
jge LoopTop ; Exit condition in the end
如果计数器最终为0,在例子12.1c与12.2b中的cmp指令可以消除,因为我们依赖add指令来设置零标记。这可以通过从n递减到0,而不是从0递增到n,来完成:
; Example 12.3. Loop with counting down
mov ecx, n ; Load n
test ecx, ecx ; Test n
jng LoopEnd ; Skip if n <= 0
LoopTop:
; (loop body) ; Loop body goes here
sub ecx, 1 ; n--
jnz LoopTop ; Loop back if not zero
LoopEnd:
现在,循环开销被减少为仅两条指令,这是最好的可能。Jecxz与loop指令都应该避免,因为它们较低效。
例子12.3中的解决方案不好,如果在循环里需要i,例如用作一个数组索引。下面的例子对一个整数数组中的所有元素加1:
; Example 12.4a. For-loop with array
mov ecx, n ; Load n
test ecx, ecx ; Test n
jng LoopEnd ; Skip if n <= 0
xor eax, eax ; i = 0
lea esi, Array ; Pointer to an array
LoopTop:
; Loop body: Add 1 to all elements in Array:
add dword ptr [esi+4*eax], 1
add eax, 1 ; i++
cmp eax, ecx ; i < n
jl LoopTop ; Loop back if i < n
LoopEnd:
数组的开始地址在esi中,索引在eax中。在地址计算中索引乘4,因为每个数组元素的大小是4字节。
修改例子12.4a,使它递减而不是递增是可能的,但数据缓存是对数据前向访问,而不是后向优化的。因此,从-n递增到0更好。这可以通过使一个指针指向数组末尾,然后使用相对数组末尾的负偏移:
; Example 12.4b. For-loop with negative index from end of array
mov ecx, n ; Load n
lea esi, Array[4*ecx] ; Point to end of array
neg ecx ; i = -n
jnl LoopEnd ; Skip if (-n) >= 0
LoopTop:
; Loop body: Add 1 to all elements in Array:
add dword ptr [esi+4*ecx], 1
add ecx, 1 ; i++
js LoopTop ; Loop back if i < 0
LoopEnd:
一个稍微不同的解决方案是n乘4,然后从-4*n递增到0:
; Example 12.4c. For-loop with neg. index multiplied by element size
mov ecx, n ; Load n
shl ecx, 2 ; n * 4
jng LoopEnd ; Skip if (4*n) <= 0
lea esi, Array[ecx] ; Point to end of array
neg ecx ; i = -4*n
LoopTop:
; Loop body: Add 1 to all elements in Array:
add dword ptr [esi+ecx], 1
add ecx, 4 ; i += 4
js LoopTop ; Loop back if i < 0
LoopEnd:
在例子12,.4与12.4c之间速度没有差别,但如果数组元素的大小不是1、2、4或8,使我们不能用比例索引取址,后者更有用。
循环计数器应该总是一个整数,因为浮点比较指令比整数比较指令要低效,即使使用SSE2指令集。某些循环天然有一个浮点退出条件。一个众所周知的例子是在项变得足够小时结束的泰勒展开。在这样的情形中,总是使用最坏情形的最大重复计数可能是有帮助的(参考第105页)。通常避免在循环中计算退出条件,以及使用整数作为循环控制,足以弥补重复循环超出所需次数的代价。这个方法的另一个优势是,循环退出分支变得更可预测。即使在循环退出分支被误预测时,使用整数计数器的误预测代价更小,因为整数指令很可能在较慢的浮点指令前执行,使这个误预测早得多得到解决。
12.2. 归纳变量
如果因为其他目的,需要循环计数器的浮点值,最好同时有一个整数计数器与一个浮点计数器。考虑制作正弦表循环的例子:
// Example 12.5a. C++ loop to make sine table
double Table[100]; int i;
for (i = 0; i < 100; i++) Table[i] = sin(0.01 * i);
这可以改为:
// Example 12.5b. C++ loop to make sine table
double Table[100], x; int i;
for (i = 0, x = 0.; i < 100; i++, x += 0.01) Table[i] = sin(x);
这里,我们有一个用于循环控制及数组索引的整数计数器i,及用于替换0.01 * i的浮点计数器x。在之前值上加0.01的计算比把i转换到浮点,然后乘上0.01要快得多。汇编实现看起来像这样:
; Example 12.5c. Assembly loop to make sine table
.data
align 8
M0_01 dq 0.01 ; Define constant 0.01
_Table dq 100 dup (?) ; Define Table
.code
xor eax, eax ; i = 0
fld M0_01 ; Load constant 0.01
fldz ; x = 0.
LoopTop:
fld st(0) ; Copy x
fsin ; sin(x)
fstp _Table[eax*8] ; Table[i] = sin(x)
fadd st(0), st(1) ; x += 0.01
add eax, 1 ; i++
cmp eax, 100 ; i < n
jb LoopTop ; Loop
fcompp ; Discard st(0) and st(1)
在这个情形里无需优化循环开销,因为速度受限于浮点计算。另一个可能的优化是使用在一个XMM寄存器中一次计算2或4个正弦值的库函数。这样的函数可以在Intel与AMD的库找到,参考手册1《优化C++软件》。在例子12.5c中通过向先前值加0.01,而不是i乘0.01,计算x的方法,通常称为使用归纳变量(induction variable)。在循环中从先前值计算某个值,比从循环计数器计算它要容易时,归纳变量是有用的。归纳变量可以是整数或浮点。归纳变量最常见的使用是计算数组地址,如例子12.4c,但归纳变量也可用于更复杂的表达式。通过n个归纳变量,只需要n次加法,不需要乘法,计算任何循环计数器的n阶多项式函数。例子参考手册1《优化C++软件》。
通过归纳变量方法计算循环计数器的一个函数,导致一个循环携带依赖链。如果这个链太长,那么从一个2次或更早迭代前的值计算每个值是有利的,如第107页泰勒展开的例子。
12.3. 移动循环不变代码
任何不在循环里改变的表达式的计算,应该移出循环。
对带有循环内不改变的条件的if-else分支,这也同样适用。可以准备两个循环,与每个分支对应,并准备一个在这两个循环间选择的分支,可以避免这样一个分支。
12.4. 找出瓶颈
有若干会限制循环性能的瓶颈。最可能的瓶颈有:
- 缓存不命中与缓存竞争
- 循环携带依赖链
- 指令获取
- 指令解码
- 指令回收
- 寄存器读暂停
- 执行端口吞吐率
- 执行单元吞吐率
- μop次优的重排与调度
- 分支误预测
- 浮点异常与次优操作数
如果一个特定瓶颈限制了性能,那么优化其他别的没有帮助。因此,仔细分析循环以识别出哪些瓶颈是限制因素,是非常重要的。仅在最窄的瓶颈被成功消除后,看下一个瓶颈才是合理的。在下面的章节讨论各种瓶颈。所有这些细节都是处理器特定的。下面提到的处理器特定细节,参考手册3《Intel,AMD与VIA CPU的微架构》的解释。
有时,找出并修复限制瓶颈需要大量的实验。记住通过实验找到的解决方案是CPU特定的,在具有不同微架构的CPU上实现最优,是不大可能的。
12.5. 循环中指令获取、解码与回收
关于如何优化指令获取、解码、回收等的细节是处理器特定的,如第57页提到的。
如果代码获取是一个瓶颈,把循环入口对齐到16,减小指令大小使循环内16字节边界尽量少,是必须的。
如果指令解码是一个瓶颈,遵守CPU特定解码模式的有关规则是必要的。避免复杂指令,如LOOP,JECXZ,LODS,STOS等。
寄存器读暂停会出现在某些Intel处理器中。如果寄存器读暂停很可能出现,重排指令或刷新在循环中读多次但不写入寄存器,是必要的。
应该避免循环内的跳转与调用,因为它推迟了代码获取。
如果可能,应该内联循环内调用的子例程。
如果可能,应该避免循环内分支,因为它们干扰循环退出分支的预测。不过,不应该使用条件移动替换分支,如果这会增加循环携带依赖链的长度。
如果指令回收是一个瓶颈,可能最好使循环内的μop总数是回收率的倍数(Core2是4,其他处理器是3)。在这个情形中,需要一些实验。
12.6. 在执行单元间平均分派μop
手册4《指令表》包含了每条指令产生多少个μop,以及每个μop去往哪个执行端口的表格。当然,这个信息是CPU特定的。计算循环一共产生多少个μop,以及多少个μop去往每个执行端口与执行单元,是必要的。
回收循环中所有指令所需的时间是μop总数除回收速率。对Intel Core2与更新处理器的回收速率是每时钟周期4个μop。计算的回收时间是循环的最小执行时间。这个值作为其他可能瓶颈可以比较的标准是有用的。
在大多数Intel处理器上,一个执行单元的吞吐率是每时钟周期1个μop。一个特定执行端口上的负荷计算为,去往这个端口的μop数除该端口的吞吐率。如果这个值超过了上面计算的回收时间,这个特定的执行端口很可能是一个瓶颈。AMD处理器没有执行端口,但它们有3或4条有类似吞吐率的流水线。
在Intel处理器上,在每个执行端口上可能有多个执行单元。大多数执行单元有与执行端口相同的吞吐率。如果是这种情形,执行单元不会是比执行端口更窄的瓶颈。但在以下情形,执行单元会成为一个瓶颈:(1)如果执行单元的吞吐率低于执行端口的吞吐率,如乘法与除法;(2)如果执行单元可通过多个执行端口访问;及(3)在没有执行端口的AMD处理器上。
一个特定执行单元的负荷计算为,去往这个执行单元的μop数乘以该单元的吞吐率倒数。如果这个值超过上面计算的回收时间,这个特定的执行单元很可能是一个瓶颈。
12.7. 分析向量循环中瓶颈的一个例子
进行这些计算的方法在以下例子中展示,它是线性代数中使用的所谓DAXPY算法:
// Example 12.6a. C++ code for DAXPY algorithm
int i; const int n = 100;
double X[n]; double Y[n]; double DA;
for (i = 0; i < n; i++) Y[i] = Y[i] - DA * X[i];
下面的实现用于32位模式中,具有SSE2指令集的处理器,假定X与Y对齐到16:
; Example 12.6b. DAXPY algorithm, 32-bit mode
n = 100 ; Define constant n (even and positive)
mov ecx, n * 8 ; Load n * sizeof(double)
xor eax, eax ; i = 0
lea esi, [X] ; X must be aligned by 16
lea edi, [Y] ; Y must be aligned by 16
movsd xmm2, [DA] ; Load DA
shufpd xmm2, xmm2, 0 ; Get DA into both qwords of xmm2
; This loop does 2 DAXPY calculations per iteration, using vectors:
L1: movapd xmm1, [esi+eax] ; X[i], X[i+1]
mulpd xmm1, xmm2 ; X[i] * DA, X[i+1] * DA
movapd xmm0, [edi+eax] ; Y[i], Y[i+1]
subpd xmm0, xmm1 ; Y[i]-X[i]*DA, Y[i+1]-X[i+1]*DA
movapd [edi+eax], xmm0 ; Store result
add eax, 16 ; Add size of two elements to index
cmp eax, ecx ; Compare with n*8
jl L1 ; Loop back
现在,让我们分析这个代码在Pentium M处理器上的瓶颈,假定没有缓存不命中。我提到的CPU特定的细节在手册3《Intel,AMD与VIA CPU的微架构》中解释。
我们仅对循环感兴趣,及L1后的代码。我们需要列出循环中所有指令的μop分解,所有手册4《指令表》中的表。这个列表看起来如下:
指令 | 融合μop | 每个执行端口的μop | 执行单元 | ||||||
|
| 端口0 | 端口1 | 端口0或1 | 端口2 | 端口3 | 端口4 | FADD | FMUL |
movapd xmm1,[esi+eax] | 2 |
|
|
| 2 |
|
|
|
|
mulpd xmm1, xmm2 | 2 | 2 |
|
|
|
|
|
| 2 |
movapd xmm0,[edi+eax] | 2 |
|
|
| 2 |
|
|
|
|
subpd xmm0, xmm1 | 2 |
|
| 2 |
|
|
| 2 |
|
movapd [edi+eax],xmm0 | 2 |
|
|
|
| 2 | 2 |
|
|
add eax, 16 | 1 |
|
| 1 |
|
|
|
|
|
cmp eax, ecx | 1 |
|
| 1 |
|
|
|
|
|
jl L1 | 1 |
| 1 |
|
|
|
|
|
|
总共 | 13 | 2 | 1 | 4 | 4 | 2 | 2 | 2 | 2 |
表12.1. 在Pentium M上DAXYP循环的总μop数
去往所有端口的μop总数是15。来自movapd [edi+eax], xmm0的4个μop被融合为2个μop,因此在融合域里的μop总数是13。总的回收时间是13 μop / 每时钟周期3 μop = 4.33个时钟周期。
现在,我们将计算去往每个端口的μop数。端口0有2个μop,端口1有1个μop,4个可以去往端口0或端口1。端口0与端口1合起来总计7个μop,因而每个端口平均每次迭代得到3.5个μop。每个端口有每时钟周期1个μop的吞吐率,因此总负载对应3.5个时钟周期。这小于回收的4.33个时钟周期,因此端口0与1将不是瓶颈。端口2每时钟周期接收4个μop,因此它接近饱和。端口3与4每个仅接收2个μop。
下一个分析涉及执行单元。FADD单元有每时钟周期1个μop的吞吐率,它可以从端口0与1接收μop。总负载是2 μop / 每时钟1 μop = 2时钟周期。因此,FADD单元远不到饱和。FMUL单元有每时钟周期0.5个μop的吞吐率,它仅能从端口0接收μop。总负载是2 μop / 每时钟0.5 μop = 4时钟周期。因此,FMUL单元接近饱和。
在循环中存在一个依赖链。时延是:内存读2时钟周期,乘法5时钟周期,加法3时钟周期,内存写3时钟周期,总共13时钟周期。这3倍于回收时间,但不是循环携带的依赖,因为每次迭代的结果保存到内存,在下一次迭代中不重用。乱序执行机制与流水线使每个计算可以在之前计算完成前开始。仅有的循环携带依赖链是add eax, 16,它的时延仅是1。
指令获取所需的时间可以从指令长度算出。在循环中8条指令的总长度是30字节。如果循环入口对齐到16,处理器需要获取两个16字节块,或者在最坏情形下。三个16字节块。因此,指令获取不是瓶颈。
在PM处理器中,指令解码遵循4-1-1模式。在例子12-6.b的循环里,指令的μop模式是2-2-2-2-2-1-1-1。这不是最优的,将需要6时钟周期解码。这超过了回收时间,因此我们可以断定,指令解码是例子12.6b里的瓶颈。总执行时间是每次迭代6时钟周期,或者每次计算Y[i]值3时钟周期。通过移动其中一个1 μop指令,以及改变xmm2的符号,使movapd xmm0, [edi+eax]与subpd xmm0, xmm1可以被合并为一条指令addpd xmm1, [edi+eax],可以改善解码。因此,我们将如下改变循环的代码:
; Example 12.6c. Loop of DAXPY algorithm with improved decoding
.data
align 16
SignBit DD 0, 80000000H ; qword with sign bit set
n = 100 ; Define constant n (even and positive)
.code
mov ecx, n * 8 ; Load n * sizeof(double)
xor eax, eax ; i = 0
lea esi, [X] ; X must be aligned by 16
lea edi, [Y] ; Y must be aligned by 16
movsd xmm2, [DA] ; Load DA
xorpd xmm2, [SignBit] ; Change sign
shufpd xmm2, xmm2, 0 ; Get -DA into both qwords of xmm2
L1: movapd xmm1, [esi+eax] ; X[i], X[i+1]
mulpd xmm1, xmm2 ; X[i] * (-DA), X[i+1] * (-DA)
addpd xmm1, [edi+eax] ; Y[i]-X[i]*DA, Y[i+1]-X[i+1]*DA
add eax, 16 ; Add size of two elements to index
movapd [edi+eax-16], xmm1 ; Address corrected for changed eax
cmp eax, ecx ; Compare with n*8
jl L1 ; Loop back
循环中的μop数如前,但解码模式现在是2-2-4-1-2-1-1,解码时间从6减少到4时钟周期,因此解码不再是一个瓶颈。
实验测试表明,没有得到预期的提升。执行时间每迭代仅从6.0降到了5.6时钟周期。为什么执行时间高于每迭代预期的4.3时钟周期,有三个原因。
第一个原因是寄存器读暂停。重排缓冲每时钟周期可以处理,来自最近没有修改的寄存器,不超过3次的读。寄存器esi,edi,ecx与xmm2没有在循环里修改。Xmm2算作两个,因为它实现为两个64位寄存器。Exa也对寄存器读暂停有贡献,因为在再次使用前,有足够的时间回收它的值。这导致了3次迭代中2次出现寄存器读暂停。将指令add eax, 16移到mulpd前面,使寄存器esi、xmm2与edi的读进一步分开,因而防止在同一时钟周期内进入重排缓冲,执行时间可以降到5.4。
第二个原因是回收。循环中融合μop数不能被3整除。因此,到L1的被采用跳转不总是进入它必须进入的回收站第一个工位。这有时会花费1时钟周期。
第三个原因是,addpd的两个μop可能分别通过端口0与端口1,在同一个时钟周期发布,尽管浮点加法仅有一个执行单元。这是坏设计的结果,如手册3《Intel,AMD与VIA CPU的微架构》解释。
一个恰好工作的更好的解决方案是,使用相对数组的末尾负索引去掉cmp指令:
; Example 12.6d. Loop of DAXPY algorithm with negative indexes
.data
align 16
SignBit DD 0, 80000000H ; qword with sign bit set
n = 100 ; Define constant n (even and positive)
.code
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea esi, [X + 8 * n] ; Point to end of array X (aligned)
lea edi, [Y + 8 * n] ; Point to end of array Y (aligned)
movsd xmm2, [DA] ; Load DA
xorpd xmm2, [SignBit] ; Change sign
shufpd xmm2, xmm2, 0 ; Get -DA into both qwords of xmm2
L1: movapd xmm1, [esi+eax] ; X[i], X[i+1]
mulpd xmm1, xmm2 ; X[i] * (-DA), X[i+1] * (-DA)
addpd xmm1, [edi+eax] ; Y[i]-X[i]*DA, Y[i+1]-X[i+1]*DA
movapd [edi+eax], xmm1 ; Store result
add eax, 16 ; Add size of two elements to index
js L1 ; Loop back
这从循环消除了一个μop。我的测量显示,在PM处理器上例子12.6d的执行时间为每迭代5.0时钟周期。理论最小值是4。寄存器读暂停消失了,因为现在eax在再次使用前回收的时间更少。回收也改善了,因为循环里融合μop数现在是12,它可被回收速率3整除。浮点加法μop冲突的问题仍然存在,这对应额外的时钟周期。这个问题仅由实验捕获。我发现最优的指令次序在mulpd后紧跟着add指令:
; Example 12.6e. Loop of DAXPY. Optimal solution for PM
.data
align 16
SignBit DD 0, 80000000H ; qword with sign bit set
n = 100 ; Define constant n (even and positive)
.code
mov eax, -n * 8 ; Index = -n * sizeof(double)
lea esi, [X + 8 * n] ; Point to end of array X (aligned)
lea edi, [Y + 8 * n] ; Point to end of array Y (aligned)
movsd xmm2, [DA] ; Load DA
xorpd xmm2, [SignBit] ; Change sign
shufpd xmm2, xmm2, 0 ; Get -DA into both qwords of xmm2
L1: movapd xmm1, [esi+eax] ;
mulpd xmm1, xmm2 ;
add eax, 16 ; Optimal position of add instruction
addpd xmm1, [edi+eax-16] ; Address corrected for changed eax
movapd [edi+eax-16], xmm1 ; Address corrected for changed eax
js L1 ;
现在,执行时间减少到每迭代4.0时钟周期,是理论最小值。例子12.6e循环中瓶颈的分析,给出了如下几个:解码时间是4时钟周期。回收时间是12 / 3 = 4时钟周期。端口0与1用到了容量的75%。端口2用到100%。端口3与4用到50%。FMUL执行单元用到最大吞吐率的100%。FADD单元用到50%。结论是,速度受限于4个相同的窄瓶颈,不可能进一步改进。
我们找到消除浮点μop冲突问题的指令次序的事实,纯属运气。在更复杂的情形里,可能不存在消除这个问题的解决方案。