怎样优化Pentium系列处理器的代码
目录
- 简介
- 文献
- 高级语言中调用汇编函数
- 调试及校验
- 内存模式
- 对齐
- Cache
- 第一次 vs 重复运行
- 地址生成互锁(AGI) (PPlain 及 PMMX)
- 配对整数指令 (PPlain 及 PMMX)
- 将复杂指令集分割为简单指令 (PPlain 及 PMMX)
- 前缀 (PPlain 及 PMMX) 以下待译...
- PPro, PII 及 PIII 流水线综述
- 指令解码 (PPro, PII 及 PIII)
- 取指令 (PPro, PII 及 PIII)
- 寄存器重命名 (PPro, PII 及 PIII)
- 乱序执行 (PPro, PII 及 PIII)
- 引退 (PPro, PII 及 PIII)
- 部分(Partial)延迟 (PPro, PII 及 PIII)
- 依赖环 (PPro, PII 及 PIII)
- 寻找瓶颈 (PPro, PII 及 PIII)
- 分支和跳转 (所有的处理器)
- 减少代码长度 (所有的处理器)
- 规划浮点代码 (PPlain 及 PMMX)
- 循环优化 (所有的处理器)
- 有问题的指令
- XCHG (所有的处理器)
- 带进位标记的循环移位 (所有的处理器)
- 字符串指令 (所有的处理器)
- 位测试 (所有的处理器)
- 整数乘法 (所有的处理器)
- WAIT 指令 (所有的处理器)
- FCOM + FSTSW AX (所有的处理器)
- FPREM (所有的处理器)
- FRNDINT (所有的处理器)
- FSCALE 及幂函数 (所有的处理器)
- FPTAN (所有的处理器)
- FSQRT (PIII)
- MOV [MEM], ACCUM (PPlain 及 PMMX)
- TEST 指令 (PPlain 及 PMMX)
- Bit scan (PPlain 及 PMMX)
- FLDCW (PPro, PII 及 PIII)
- 特别主题
- 指令速度列表 (PPlain 及 PMMX)
- 指令速度及微操作失败列表(PPro, PII 及 PIII)
- 速度测试
- 不同的微处理器间的比较
1. 简介
这本手册细致的描述了怎样写出高度优化的汇编代码, 着重于讲解Pentium®系列的微处理器.
这儿所有的信息都基于我的研究. 很多人为这本手册提供了有用的信息和错误矫正, 而我在获得任何新的重要信息后都更新它. 因此这本手册比其它类似的信息来源都更准确,详尽,精确和便于理解, 而且它还包含了许多其它地方找不到的细节描述. 这些信息使你能够用多种方法精确统计一小段代码花掉的时钟周期数. 但是,我不能保证手册里所有的信息都是精确的: 一些时间测试等是很难或者不可能精确测量的, 我看不到 Intel 手册作者拥有的内部技术文档资料.
这本手册讨论了 Pentium 处理器的下列版本:
| 缩写 | 名字 |
| PPlain | plain 老式 Pentium (没有 MMX) |
| PMMX | 有MMX的Pentium |
| PPro | Pentium Pro |
| PII | Pentium II (包括 Celeron 和 Xeon) |
| PIII | Pentium III (包括一些相当的CPU) |
这本手册中使用了MASM 5.10的汇编语法. 没有什么官方的 X86 汇编语言, 但这是最接近你能接受的标准, 因为几乎所有的汇编器都有 MASM 5.10 兼容模式. (然而我不推荐使用 MASM 的5.10 版本. 因为它在 32 位模式下有严重的 Bug. 最好是使用 TASM 或者MASM 的后续版本).
手册里的一些评语好象是对Intel的批评. 但这并不是说其它的产品会好一些. Pentium系列的微处理器能在评测分中取得更好一点的比较值,有更好的文档, 和更多的可测试特性. 由于这些原因, 不会有我或者其他人做同类商品的比较测试.
汇编语言编程比用高级语言要复杂的多. 制造 Bug 是很容易的,但是找到 Bug 却很难. 现在已经提醒你了! 我假定读者已经有汇编编程的经验. 没有的话,请在做复杂的优化前读一些汇编的书并且写些代码获得些汇编的经验.
PPlain 和 PMMX 芯片的硬件设计有许多专门为一些常用指令或者指令对设计的特性, 而不是使用那些一般的优化模块. 因此,为这个设计优化软件的规则很复杂, 且有很多的例外, 但是这样做可能获得实质性的好处. PPro, PII 和 PIII 处理器有非常不同的设计,它们会利用乱序执行来做许多的优化工作, 但是处理器的这些个设计带来了许多潜在的瓶颈, 因此为这些处理器进行手工优化将得到许多的好处. Pentium 4 处理器也用了另外一种设计, 奔腾4 的优化指导路线和前面的版本非常的不同了. 这个手册没有禳括奔腾4 - 读者请自己查阅 Intel 的手册.
在把你的代码转为汇编的之前,确认你的算法是足够优化的. 通常你可以通过优化算法来将代码效率提高的比转成汇编获得的效率多的多.
第二,你必须找到你的程序里最关键的部分. 通常 99% 的 CPU 时间花在程序最里面的循环中. 在这种情况下, 你只要优化这个循环并把其它的所有东西都用高级语言写. 一些汇编程序员将大量的精力花在了他们程序的错误的部分上, 他们努力得到的唯一结果就是程序变的更加难以调试和维护了.
如果你的程序的关键部分并不那么明显,你可以用profil来找. 如果发现瓶颈在磁盘操作, 然后你就可以试着修改程序使磁盘操作集中连续, 提高磁盘缓冲的命中率, 而不是用汇编来写代码. 如果瓶颈在图象输出,那么你就可以尝试找到一种方法来减少调用图象函数的次数. 一些高级语言编译器对于指定的处理器提供了相对好的优化, 但是手工优化将做的更好.
请不要将你的编程问题寄给我. 我不会帮你做家庭作业的!
祝你在后面的阅读中好运!
2. 文献
在 Intel 的 www 站上,打印的文本或者CD-ROM 上都有很多有用的文献和教程. 建议你研究一下这些文档来对微处理器的结构有些认识. 然而,Intel 的文档也不总是对的- 尤其是那些教程有很多错误(显然,Intel的那些人没有测试他们的例子).
这里我不给出 URL, 因为文件的位置经常的改变. 你可以利用 developer.intel.com 或者www.agner.org/assem 上的链接上的搜索工具找到你要的文档.
一些文档是.PDF格式的. 如果你没有显示或者打印PDF的工具,可以去www.adobe.com下载Acrobat 文件阅读器.
使用 MMX 和 XMM (SIMD) 指令优化专门的程序在几本使用手册里都有描述. 各种手册和教程都有描述其指令集.
VTUNE 是 Intel 用来优化代码的软件工具我没有测试它,因此这里不于评价.
还有很多比 Intel 的更多的更有用的信息. 在新闻组 comp.land.asm.x86 的 FAQ 里有把这些资源列出来. 其它的internet 上的资源在 www.agner.org/assem 上也有链接.
3. 在高级语言里调用汇编函数
你可以使用在线汇编或者用汇编写整个子程序然后再连接到你的工程中. 如果你选择后者, 建议你用可以将高级语言直接编译成汇编的编译器. 这样你可以得到函数的正确的调用原型. 所有的 C++ 编译器都能做这个工作.
传递参数的方法取决于调用形式:
| 调用方式 | 参数在堆栈里的次序 | 参数由谁来移去 |
| _cdecl | 第一个参数在低位地址 | 调用者 |
| _stdcall | 第一个参数在低位地址 | 子程序 |
| _fastcall | 编译器指定 | 子程序 |
| _pascal | 第一个参数在高位地址 | 子程序 |
函数调用原型和被编译器命名的函数名可能非常的复杂. 有很多不同的调用转换规则, 不同的编译器也互不兼容. 如果你从C++里调用汇编语言的子程序,最好的方法是将你的函数用 extern"C" 和 _cdel 定义来兼容性和一致性. 汇编代码的函数名前面必须带一个下划线 (_) 并且在外面编译时加上大小写敏感的选项 (选项 -mx). 例如:
; extern "C" int _cdecl square (int x);
_square PROC NEAR ; 整型平方函数
PUBLIC _square
MOV EAX, [ESP+4]
IMUL EAX
RET
_square ENDP
如果你需要重载函数,重载操作符,方法, 和其它 C++ 专有的东西, 就必须在 C++ 里先写好代码再用编译器编译成汇编代码以获得正确的连接信息和调用模型. 这些随着编译器的不同而不同而且很少列出文档. 如果你希望汇编函数用其它的调用原型而不是extern "C" 及 _cdecl 的, 又可以被不同的编译器调用, 那么你需要为每个编译器写一个名字. 例如重载一个 square 函数:
; int square (int x);
SQUARE_I PROC NEAR ; 整数 square 函数
@square$qi LABEL NEAR ; Borland 编译器的连接名字
?square@@YAHH@Z LABEL NEAR ; Microsoft 编译器的连接名字
_square__Fi LABEL NEAR ; Gnu 编译器的连接名字
PUBLIC @square$qi, ?square@@YAHH@Z, _square__Fi
MOV EAX, [ESP+4]
IMUL EAX
RET
SQUARE_I ENDP
; double square (double x);
SQUARE_D PROC NEAR ; 双精度浮点 square 函数
@square$qd LABEL NEAR ; Borland 编译器的连接名字
?square@@YANN@Z LABEL NEAR ; Microsoft 编译器的连接名字
_square__Fd LABEL NEAR ; Gnu 编译器的连接名字
PUBLIC @square$qd, ?square@@YANN@Z, _square__Fd
FLD QWORD PTR [ESP+4]
FMUL ST(0), ST(0)
RET
SQUARE_D ENDP
这个方法能够工作是因为所有这些编译器对重载的函数都缺省使用_cdecl 调用. 然而,对于方法(成员函数), 对于各种编译器, 甚至调用方式都不一样 (Borland 和 Gnu 编译器使用 _cdecl 方式, 'this' 指针是第一个参数, 而Microsoft 使用 _stdcall 方式, 而 'this' 指针放在 ecx 里).
通常来说, 当你使用了下列东西时, 不要指望不同的编译器在目标文件级别可以兼容: long double, 成员指针, 虚机制, new, delete, 异常, 系统函数调用, 或者标准库函数.
16 位模式 DOS 和Windows, C or C++ 的寄存器使用:
AX 是 16 位返回值, DX:AX 是 32 位返回值, ST(0) 是浮点返回值.寄存器 AX, BX, CX, DX, ES 数学标记可以被过程改变; 其它的寄存器必须保存和恢复. 一个过程要不改变SI, DI, BP, DS 和 SS 的前提下才不会影响另一个过程..
32 位模式 Windows, C++ 和其它编程语言下的寄存器使用:
整数返回值放在 EAX, 浮点返回值放在ST(0). 寄存器 EAX, ECX, EDX(没有 EBX) 可以被过程修改; 其它的寄存器必须保留和恢复.段寄存器不能被改变,甚至不能被临时改变. CS, DS, ES, 和 SS 都指向当前的段组. FS 被操作系统使用, GS 没有使用, 但是被保留. 标记位可以在下面的限制下被过程改变: 方向标志缺省是 0.方向标志可以暂时的修改, 但是必须在任何的调用或者返回前清除.中断标志必须清除. 浮点寄存器堆栈在过程入口处是空的, 返回时也应该是空的, 除了 ST(0) 如果被用于返回值则除外. MMX 寄存器可以被改变但是在返回前或者在调用可能使用浮点运算的过程前必须用EMMS 清一下. 所有的 XMM 寄存器都可以被过程修改. 在XMM 寄存器里的传递参数和返回值的描述在 Intel 的应用文档 AP 589. 一个过程可以在不改变EBX, ESI, EDI, EBP 和所有的段寄存器的前提下被另一个过程调用.
4. 调试和校验
正如你已经发现的, 调试汇编代码非常的困难和容易受到挫折. 我建议你先把你需要优化的小段代码用高级语言写成一个子程序. 然后写个小的测试程序可以充分测试你的这个子程序. 确认测试程序可以测试到所有的分支.
当高级语言的子程序可以工作了,你再把它翻译成汇编代码.
现在你可以开始优化了.每次你做了点修改都应该运行测试程序看看能不能正确工作. 将你所有的版本都标上号并保存起来,这样在发现测试程序检查不到的错误 (比如写到错误的地址)时可以回头来重新测试.
用 第30章 里提到的所有方法或者用测试程序测试最你的程序中最关键的部分.如果代码比你期望的速度慢的太多, 最可能的原因是: cache 失效 (第7章), 未对齐操作 (第6章), 第一次运行消耗(第8章), 分支预测失败(第22章)取指令问题, (第15章), 寄存器读延迟(第16章), 或者是过长的依赖环(第20章).
高度优化的代码将变得对其他人非常难读懂, 甚至对你日后在读也有困难.为了使维护代码变为可能, 将代码组织为一个个小的良好预定义宏和清楚注释的逻辑段(过程或者宏)就非常重要. 代码越复杂艰涩, 写下好的文档就越重要.
5. 内存模式
Pentium 主要为 32 位代码设计,16位代码的性能很差. 将你的代码和数据分段也会明显的降低性能,因此通常你应当使用32位平坦模式,并且使用支持这种模式的操作系统. 如果不特别注明,这本手册里所有的例子都使用32位平坦内存模式.
6. 对齐
内存里的所有数据都必须按照下表将地址对齐到可以被2,4,8 或 16 整除的位置:
|
| 对齐 | |
| 操作数据长度 | PPlain 及 PMMX | PPro, PII 及 PIII |
| 1 (byte) | 1 | 1 |
| 2 (word) | 2 | 2 |
| 4 (dword) | 4 | 4 |
| 6 (fword) | 4 | 8 |
| 8 (qword) | 8 | 8 |
| 10 (tbyte) | 8 | 16 |
| 16 (oword) | n.a. | 16 |
在 PPlain 和 PMMX 上, 访问未对齐数据在4字节边界线交错的时候将至少有 3个时钟周期的额外消耗. 当cache 边界线被交错的时候损耗更大.
在 PPro, PII 和 PIII 上, 当cache边界线交错时, 未对齐数据将消耗掉 6-12 个时钟周期. 小于 16 字节的未对齐操作数, 没有在 32 字节边界上交错时将没有额外的损耗.
在 dword 栈上以 8 或 16 对齐数据可能会有问题. 常见的情况是设置对齐的帧指针. 对齐本地数据的函数可以是这样:
_FuncWithAlign PROC NEAR
PUSH EBP ; 前续代码
MOV EBP, ESP
AND EBP, -8 ; 以 8 来对齐帧指针
FLD DWORD PTR [ESP+8] ; 函数参数
SUB ESP, LocalSpace + 4 ; 分配本地空间
FSTP QWORD PTR [EBP-LocalSpace] ; 在对齐了的空间保存一些东西
...
ADD ESP, LocalSpace + 4 ; 结束代码. 恢复 ESP
POP EBP ; (PPlain/PMMX 上有 AGI 延迟)
RET
_FuncWithAlign ENDP
虽然对齐数据永远是重要的, 但是在 PPlain 和 PMMX 上对齐代码却没有必要. PPro, PII 及 PIII 上对齐代码的原则在第15章阐述.
7. Cache
PPlain 和 PPro 有代码用的8 kb 片内 cache (一级 cache), 8 kb 数据 cache. PMMX, PII 和 PIII 则有 16 kb 代码 cache 和 16kb 数据 cache. 一级 cache 里的数据可以在一个时钟周期内读写, cache 未命中时将损失很多时钟周期. 理解 cache 是怎样工作的非常重要, 这样才能更有效的使用它.
数据 cache 由 每行 32 字节的 256 或 512 行组成. 每次你读数据未命中, 处理器将从内存读出一整条 cache 行. cache 线总是在物理地址的 32 字节对齐. 当你从一个可以被32 整除的地址读出一个字节, 下 31 字节的读写就不会有多余的消耗. 你可以对齐数据在 32 字节块里来从中获得好处. 例如, 如果你有一个循环要操作两个数组, 你就可以将两个数组合并成一个结构数组, 这样被使用的数据就储存在一起了.
如果数组或者其它数据结构是 32 字节的倍数, 你最好将其按 32 字节对齐.
cache 是联想设置的. 这就是说cache 行不能随心所欲的设置到指定内存地址. 每个cache 行有 7-bit 设定值来匹配物理地址的 5 到 11 位 (0-4 位就是cache行的 32 字节). PPlain 和PPro 可以有两条cache行对应 128 个设定值中的一个值, 因此任两条 cache 行可以设定到任何的 RAM 地址. PMMX, PII 和PIII 则可以有 4 个.
其结果是 cache 不能保留超过2个或4个的地址的 5-11 位相同的不同数据块. 你可以用用以下方法检测两个地址是否有相同的设定值: 截掉地址的低5位得到可以被 32 整除的数值. 如果两个截断地址之差是 4096 (=1000H) 的倍数, 这两个地址就有相同的设定值.
让我用下面的一小段代码来说明一下, 这里 ESI 放置了一个可以被 32 整除的地址:
AGAIN: MOV EAX, [ESI]
MOV EBX, [ESI + 13*4096 + 4]
MOV ECX, [ESI + 20*4096 + 28]
DEC EDX
JNZ AGAIN
这 3 个地址都有相同的设定值, 因为不同的截断地址都是 4096 的倍数. 这个循环在 PPlain 和PPro 上将运行的相当慢. 当你读 ECX 的时候, 没有适当设定值的空闲 cache 行,因此处理器将使用最近最少使用的两组 cache 行中的一个, 那就是EAX 使用的那个, 然后从[ESI+20*4096] 到 [ESI+20*4096+31]读出数据来填充并读入 ECX. 下一步, 当读EAX时, 你将发现为EAX保存数据的cache行已经被丢弃了, 所以将使用最近最少使用的 cache 行, 那就是保存 EBX 数值的那个了, 等等.. 这将会产生大量的 cache 丢失, 这个循环大概开销60 个时钟周期. 如果第3行改成:
MOV ECX, [ESI + 20*4096 + 32]
这样我们就会在 32 字节边界上交错, 因此我们没有和前两行相同的设定值, 这样为这三个地址分别指定cache行就没有什么问题了. 这个循环仅仅消耗3个时钟周期(除了第一次运行) - 一个相当大的提高! 如刚才提到的,PMMX, PII 和PIII 有 4 向 cache行, 因此你可以有4个相同设定值的cache行.(一些Intel文档错误的说 PII的cache是两向的).
检测你的数据地址是否有相同的设定值可能非常困难, 尤其是它们分散在不同的段里. 你能避免这种问题的最好能做的是将关键部分使用的所有数据都放在一个不超过cache大小的连续数据块里, 或者放在两个不超过cache一半大小连续数据块 (例如一个静态数据块, 一个堆栈数据块)这样你的cache行就一定会高效使用.
如果你的代码的关键部分要操作很大的数据结构或者随机数据地址, 你可能会想保存所有常用的变量(计数器,指针,控制变量等) 在一个单独的最大为 4k的连续块里面, 这样你就有一个完整的cache行集来访问随机数据. 既然你通常需要栈空间来为子程序保存参数和返回地址, 最好的做法是复制所有的常用静态数据到堆栈, 如果它们被改变, 就在关键循环外复制回去.
读一个不在一级缓存里的数据将导致从二级缓存读入整个cache行, 这大约要消耗 200ns (在100MHz 系统上是 20 时钟周期, 或是 200MHz 上的 40 个周期), 但是你最先需要的数据将在 50-100 ns后准备好. 如果数据也不在二级缓存, 你将会碰到 200-300 ns 的延迟. 如果数据在 DRAM 页边界交错, 延迟时间会更长一些. (4/8 Mb 的 72 线内存DRAM 页大小是 1Kb, 16/32 Mb 的是2kb).
当从内存读入大块的数据, 数据限制在于填充 cache 行. 有时你可以以非连续的次序读取数据来提高速度: 在你读完一个cache行之前就开始读下一个cache行的第一个数据. 这个方法可以提高从主内存或PPlain及PMMX的二级cache,和 PPro,PII,PII的二级cache读入的速度 20 - 40%. 这个方法的不利地方在于使程序代码变的非常的笨拙和难于理解. 关于这些的更多信息请参考www.intelligentfirm.com.
当你写向一个不在一级cache的地址, 在 PPlain和PMMX上这个数值将直接写到二级cache或者是 RAM(这取决于2级cache如何设置). 这大约消耗100 ns. 如果你向同一个32字节的内存块反复写8次或8次以上也没有从里面读, 而且这个块不在一级缓存, 这样就会促使对这个块的一次哑读而使其加载入一个cache行. 所有随后的向这个块的写操作就会被定向到cache里, 每次只消耗一个时钟周期. 在PPlain 和 PMMX 上, 有时会因为重复写向一个地址而在其间没有读它而带来小小的惩罚.
在 PPro, PII 和 PIII上, 一次写失误会读入一个cache行,但是可能设置一个内存区域做不同的操作, 例如显存 (见 Pentium Pro 系列开发者手册, vol. 3: 操作系统写作者指南").
提高内存读写速度的方法在下面的 27.8章节讨论.
PPlain 和 PPro 有两个写缓存,PMMX, PII 和 PIII 有四个. 在 PMMX, PII 和 PIII 上你可以做四次未完成的写操作不cache内存而在后面的操作中不产生延迟. 每个写缓存可以控制 64 位宽.
临时数据可以方便的放在堆栈里, 因为堆栈区域非常有可能在cache中. 然而,如果你的数据元素大于堆栈字大小时应该注意对齐问题.
如果两个数据结构的生命期不重叠的话, 共享相同的 RAM 区域可以提高cache效率. 一般都在堆栈为临时变量分配空间.
将临时数据保存在寄存器里将有更高的效率. 既然寄存器是一种稀有资源, 你可能想用[ESP]而不是[EBP]来定位堆栈里的数据, 这样就可以释放EBP用于其它用途. 不要忘记了 ESP 将在你每次做 PUSH 或者 POP 时都会被改变. (你不能在 16位 Windows下使用ESP, 因为时钟中断将把ESP的高字设置到你的代码中无法预测的位置.)
有一个分开的cache给代码使用, 它和数据cache是类似的. 代码cache 的大小在 PPlain和PPro上是8kb, 在 PMMX,PII和PIII上是16kb. 能让你的代码的关键部分(最里面的循环)放入代码cache是很重要的. 最常用的需要一起使用的代码或者例程最好是储存在临近的位置. 不常用的分支或者过程可以放在代码的下面或者其它的位置.
8. 第一次 vs 重复运行
一片代码往往在第一次运行时比重复运行消耗更多的时间. 原因见下:
- 从 RAM 读入代码到cache花去了比运行它更多的时间.
- 代码操作的所有数据都必须加载到cache, 这比执行那些操作更花时间. 当代码重复运行的时候, 数据几乎都在 cache 里.
- 跳转指令在第一次运行的时候并不在分支目的缓存里, 因此一般都不能正确的预定向. 见第22章.
- PPlain 上, 代码的解码是个瓶颈. 如果花掉一个时钟周期去检测指令长度, 那么就不可能在一个时钟周期解码两条指令, 因为处理器不知道第二条指令从那里开始的. PPlain 从上次运行时间记住每条保存在cache里的指令的长度来解决这个问题. 这样做的结果是, PPlain 上第一次的指令如果不是只有一个字节长的话就不会配对执行. PMMX, PPro, PII 和 PIII 在第一次解码却没有这个问题.
因为这四个原因,在循环中的一段代码通常比随后的运行花去更多的时间.
如果你使用了一个很大的循环而不能放入代码cache, 将导致效率下降, 因为它们不能在 cache 运行. 因此你应该重新组织一下循环使cache能放下它们.
如果你有非常多的跳转,调用,分支在循环里, 就会反复的产生分支目的缓存失败.
同样的, 如果循环重复操作一个比数据缓存大的数据结构, 也会在所有的时间都会未命中数据cache.
9. 地址生成互锁(AGI) (PPlain and PMMX)
指令操作内存所需要的地址需要一个时钟周期来计算. 通常在前面的指令或是指令对运行的时候它已经在流水线上计算好了. 但是如果地址的计算倚赖上个时钟周期的运行结果的话, 你就需要一个额外时钟周期来等待地址的计算. 这就叫做AGI延迟. 例如:
ADD EBX,4 / MOV EAX,[EBX]; AGI 延迟
例子里的延迟可以向 ADD EBX,4 and MOV EAX,[EBX] 间增加一些其它的指令或者重新写成 MOV EAX,[EBX+4]/ ADD EBX,4 来去掉.
当你隐性的使用ESP寻址, 比如 PUSH, POP, CALL,and RET, 且 ESP 在前个周期被MOV, ADD, 或 SUB 等修改，这样也会造成AGI延迟. PPlain 及 PMMX 有专门的电路来预测栈操作后的ESP值, 因此你在用PUSH, POP, 或 CALL 改变 ESP 后不会遇到AGI延迟. 在RET 后面, 仅仅在立即对ESP做加操作时才回产生AGI延迟.
例如:
ADD ESP,4 / POP ESI ; AGI 延迟
POP EAX / POP ESI ; 无 stall, 配对
MOV ESP,EBP / RET ; AGI 延迟
CALL L1 / L1: MOV EAX,[ESP+8] ; 无 stall
RET / POP EAX ; 无 stall
RET 8 / POP EAX ; AGI 延迟
当 LEA 指令使用了基寄存器或索引寄存器, 而它们在在前个时钟周期被改变了,同样会发生AGI延迟. 例如:
INC ESI / LEA EAX,[EBX+4*ESI] ; AGI 延迟
PPro, PII 和 PIII 在读内存和LEA上没有AGI 延迟, 但是在写内存时会有 AGI 延迟. 如果后来的的代码并不是必须等待写操作结束的话这并无大影响.
10. 整数指令配对(PPlain 及 PMMX)
10.1 完美的配对
PPlsin 及PMMX有两条流水线来运行指令, 分别叫做 U-管道和V-管道. 在特定的条件下两条指令可以一个在U-管道一个在V-管道同时运行. 这可以使速度加倍. 因此将你的指令重新组织一下次序使它们配对是很有利的.
下面这些指令可以在任意的管道内配对:
- MOV 寄存器, 内存, 或是立即数到寄存器或内存
- PUSH 寄存器或立即数, POP 寄存器
- LEA, NOP
- INC, DEC, ADD, SUB, CMP, AND, OR, XOR,
- 还有一些形式的 TEST (见 26.14章)
下面的指令只能在 U-管道配对:
- ADC, SBB
- SHR, SAR, SHL, SAL 移动立即数位
- ROR, ROL, RCR, RCL 移动立即数1位
下面的指令可以在任何管道运行,但是只能在 V-管道配对:
- near call
- short 及 near jump
- short 及 near 条件跳转.
所有的整数指令到可以在 U-管道运行, 但不一定可以配对.
两条连续的指令满足了下面的要求时就可以配对:
1. 第一条指令在 U 管道中，第二条指令在 V 管道中,且它们是可配对的.
2. 当第一条指令写一个寄存器的时候, 第二条指令不去读写它.
例如:
MOV EAX, EBX / MOV ECX, EAX ; 写后面跟着读, 不能配对
MOV EAX, 1 / MOV EAX, 2 ; 写后面跟着写, 不能配对
MOV EBX, EAX / MOV EAX, 2 ; 读后面跟着写, 可以配对
MOV EBX, EAX / MOV ECX, EAX ; 读后面跟着读, 可以配对
MOV EBX, EAX / INC EAX ; 读后面跟着读写，可以配对
3. 在第2条规则里面, 寄存器的一部分作为整个寄存器来对待, 例如:
MOV AL, BL / MOV AH, 0
写入相对寄存器的不同部分, 不能配对
4. 当两条指令同时写的是标记寄存器的不同部分时, 规则2和3都可以忽略掉. 例如:
SHR EAX, 4 / INC EBX ; 可以配对
5. 一个写标记寄存器的指令可以和一个条件跳转配对, 而忽略掉规则 2 . 例如:
CMP EAX, 2 / JA LabelBigger ; 可以配对
6. 下面的指令对, 虽然同时修改了栈指针, 但是它们依然可以配对:
PUSH + PUSH, PUSH + CALL, POP + POP
7. 对于有前缀的配对指令有一些限制. 下面列出了几种形式的前缀:
- 对非缺省段寻址的指令有段前缀.
- 当 32 位代码中使用 16 位的数据, 或16位的代码中使用 32 位数据的指令, 操作数有大小前缀.
- 16位模式中, 使用 32 位基地址或 32 位索引寄存器的指令有地址大小前缀.
- 重复性的字符串操作指令有一个重复前缀.
- 锁定指令有一个 LOCK 前缀.
- 很多在 8086 处理器中没有实现的指令有两个字节的操作码, 其中第一个字节是 0FH. 这个 0FH 字节在 PPlain 上跟前缀的表现相同, 但是后来的版本中就不会了. 最常见的带 0FH 的指令有: MOVZX, MOVSX, PUSH FS, POP FS, PUSH GS, POP GS, LFS, LGS, LSS, SETcc, BT, BTC, BTR, BTS, BSF, BSR, SHLD, SHRD, 还有带两个操作数且没有立即数的 IMUL.
在 PPlain 上, 有前缀的指令除了近距离条件跳转外只可以在 U 管道中执行.
PMMX 上, 带有操作数大小, 地址大小, 或 0FH 前缀的指令可以在任意管道运行, 但是带有段前缀, 重复前缀, 或者锁定前缀的指令还是只能在 U 管道运行 .
8. 带有一个间接数和一个立即数的指令在 PPlain 上不能配对, 而在PMMX 上只能在 U 管道配对:
MOV DWORD PTR DS:[1000], 0 ; 不能配对, 或者只能在 U 管道配对
CMP BYTE PTR [EBX+8], 1 ; 不能配对, 或者只能在 U 管道配对
CMP BYTE PTR [EBX], 1 ; 可以配对
CMP BYTE PTR [EBX+8], AL ; 可以配对
(关于一个间接数和一个立即数在 PMMX 上配对的另一个问题是 这条指令可能比 7 字节还要长, 这意味着, 只有一个时钟周期只有一条指令能够被解码, 这些放在第12章解释.)
9. 两条指令必须已经预读进来且被解码. 这些放在第 8 章解释.
10. PMMX 上, 对于 MMX 指令有特殊的配对规则:
- MMX 移位, pack 和 unpack 指令可以在任意的管道执行, 但是不能跟另外一条 MMX 移位, pack 和 unpack 指令配对.
- MMX 乘法指令可以在任意管道运行, 但是不能和另外一条 MMX 乘法指令配对. 乘法指令需要消耗 3 个时钟周期, 而后两个时钟周期可以和别的指令同时执行, 就好象浮点指令那样 (参见第 24 章).
- 一条访问内存或整数寄存器的 MMX 指令只能在 U 管道运行, 而且不能跟非 MMX 指令配对.
10.2 有缺陷配对
有几种情况下, 两条成对指令不能平行运行, 或者只是时间上部分重叠. 然而它们依然被当作是成对的, 因为第一条指令在 U 管道运行, 而第二条在 V 管道. 随后的指令要在两条有缺陷配对的指令都完成后才开始运行.
有缺陷配对发生在以下条件下:
1. 如果第二条指令遭遇了一个 AGI 延迟 (见第9章).
2. 两条指令不能同时访问内存的同一个 DWORD. 下面的例子假定 ESI 可以被 4 整除:
MOV AL, [ESI] / MOV BL, [ESI+1]
两个操作数是在同一个 DWORD 里, 因此它们不能同时执行. 这对指令需要2 个时钟周期.
MOV AL, [ESI+3] / MOV BL, [ESI+4]
这里两个操作数分别处于两个 DWORD 的边界, 因此它们完美的配对, 只需要消耗 1个时钟周期.
3. 第 2 条款可以扩展到两个地址的 2-4 位相同的情况 (cache 区冲突). 对于DWORD 地址, 这意味着两个地址差不能被 32 整除. 例如:
MOV [ESI], EAX / MOV [ESI+32000], EBX ; 有缺陷配对
MOV [ESI], EAX / MOV [ESI+32004], EBX ; 完美配对
不访问内存的可配对整数指令一个时钟周期可以执行完，但是中断预测失败的跳转例外. MOV 读或写内存的指令当数据区在 cache 里并严格对齐的时候也只需要一个时钟周期. 即使是像用索引寄存器变址的复杂寻址模式下也不会有速度上的惩罚.
一组配对整数指令, 如果需要读内存, 做一些计算后把结果保存在寄存器或标记寄存器中时, 需要消耗两个时钟周期. (读/修改 指令).
一组配对整数指令, 如果需要读内存, 做一些计算后把结果回写到内存中, 需要消耗3个时钟周期. (读/修改/写 指令).
4. 如果一条 读/修改/写 指令和一条 读/修改或 读/修改/写指令 配对,那么它们就是一个有缺陷配对.
下表展示了各种情况下需要的时钟周期数:
| 第一条指令 | 第二条指令 | ||
|
| MOV 或者 仅仅是寄存器操作 | 读/修改 | 读/修改/写 |
| MOV 或仅仅是寄存器操作 | 1 | 2 | 3 |
| 读/修改 | 2 | 2 | 3 |
| 读/修改/写 | 3 | 4 | 5 |
例如:
ADD [mem1], EAX / ADD EBX, [mem2] ; 4 个时钟周期
ADD EBX, [mem2] / ADD [mem1], EAX ; 3 个时钟周期
5. 当两条配对指令都因为cache失效,没有对齐,或跳转预测失败等情况而需要极端的时间时，一对指令消耗的时间将比其中任何一条需要的时间都长, 但是比两条指令需要的时间之和短.
6. FXCH 之后的可配对浮点指令, 当下一条指令不是一条浮点指令时 组成一个缺陷配对.
为了避免有缺陷配对, 你必须知道哪条指令进入了 U 管道, 哪条进入了 V 管道. 你可以向前看看你的代码, 找到哪条指令是不能配对的, 或者只能在一条管道中配对, 又或因为上面提及的规则而不能配对.
有缺陷配对通常可以通过重组你的指令来避免. 例如:
L1: MOV EAX,[ESI]
MOV EBX,[ESI]
INC ECX
这里两条 MOV 指令组成了一个有缺陷配对, 因为它们访问了同一内存地址, 所以这组指令需要消耗 3 个时钟周期. 你可以通过重组指令,把 INC ECX 跟其中一个 MOV 指令配对.
L2: MOV EAX,OFFSET A
XOR EBX,EBX
INC EBX
MOV ECX,[EAX]
JMP L1
INC EBX / MOV ECX,[EAX] 这对指令是一个有缺陷配对, 因为 is imperfect because the latter 因为后一条指令发生了 AGI 延迟. 这组指令消耗 4 个时钟周期. 如果你插入一条 NOP 或任意别的指令, 替换成 MOV ECX,[EAX] 跟 JMP L1 配对, 这样这组指令就只需要消耗 3 个时钟周期了.
下一个例子是16 位模式下的, 假设 SP 可以被 4 整除:
L3: PUSH AX
PUSH BX
PUSH CX
PUSH DX
CALL FUNC
这里 PUSH 指令组成了两个有缺陷配对, 因为各对指令中的两个操作数都放入了内存的同一 DWORD 中. PUSH BX 可能可以和PUSH CX 完美配对起来 (因为它们访问的是两个不同的 DWORD) 但是并不是这样, 因为它已经和 PUSH AX 配对了. 这组指令消耗了 5 个时钟周期. 如果你插入一个NOP 或者其它指令, 让PUSH BX 跟 PUSH CX 配对, 而 PUSH DX 和 CALL FUNC 配对, 这样这组指令就只需要 3 个时钟周期了. 另一个解决方案是, 让 SP 不被 4 整除. 想知道SP 是否被 4 整除在 16 位模式下是很困难的, 所以避免这个问题的最佳方案是去使用 32 位模式.
11. 将复杂指令集分割为简单指令 (PPlain 及 PMMX)
你可以把 读/修改 和 读/修改/写 指令切开来提高配对效率. 例如:
ADD [mem1],EAX / ADD [mem2],EBX ; 5 个时钟周期
这个代码可以切开, 而只需要消耗 3 个时钟周期:
MOV ECX,[mem1] / MOV EDX,[mem2] / ADD ECX,EAX / ADD EDX,EBX
MOV [mem1],ECX / MOV [mem2],EDX
同样的你可以把不能配对的指令切开让它们可以配对:
PUSH [mem1]
PUSH [mem2] ; 不能配对
切开变成:
MOV EAX,[mem1]
MOV EBX,[mem2]
PUSH EAX
PUSH EBX ; 所有的都配对了
下面还有另一些例子, 展示了一些不能配对的指令切开后变成可配对指令:
CDQ 切成: MOV EDX,EAX / SAR EDX,31
NOT EAX 改为 XOR EAX,-1
NEG EAX 切成 XOR EAX,-1 / INC EAX
MOVZX EAX,BYTE PTR [mem] 切成 XOR EAX,EAX / MOV AL,BYTE PTR [mem]
JECXZ 切成 TEST ECX,ECX / JZ
LOOP 切成 DEC ECX / JNZ
XLAT 改为 MOV AL,[EBX+EAX]
如果切开指令并不能提高速度, 你应该保持复杂指令或不能配对的指令, 这样可以减小代码的尺寸.
对于 PPro, PII and PIII, 不需要将指令切开, 除非能产生更小的代码.
英特尔多核平台编码优化大赛的优化过程
英特尔最近举办了一次多核平台编码优化大赛，我也参加了这次比赛；大赛代码提交阶段已经结束，
所以也可以公开自己的优化过程；
去年末双核处理器出货量开始超过单核处理器出货量，英特尔也发布了4核CPU；AMD今年也将发布4核，
并计划今年下半年发布8核，多核已经大势所趋；而以前的CPU性能增长，软件界都能很自然的获得相应的
性能增长，而对于这次多核方式带来的CPU性能爆炸性增长，软件界的免费午餐没有了，程序开发人员必须
手工的完成软件的并行化；而很多程序员、计算机教育界等对并行编程的准备并不充足(包括我)；英特尔
在这种背景下举办多核平台编码优化大赛是很有现实意义的；
本文章主要介绍了自己参加竞赛过程中，对代码的优化过程；
我测试和优化过程中用的 CPU：AMD64x2 3600+ (双核CPU) (开发初期还使用过AMD3000+)
操作系统：Windows XP 32bit
编译器：Visual Studio 2005
大赛公布的原始代码:
(组织者后来要求计算和输出精度到小数点后7位，这里的输出代码做了相应的调整。)
//* compute the potential energy of a collection of */
//* particles interacting via pairwise potential */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <time.h>
#define NPARTS 1000
#define NITER 201
#define DIMS 3
int rand( void );
int computePot(void);
void initPositions(void);
void updatePositions(void);
double r[DIMS][NPARTS];
double pot;
double distx, disty, distz, dist;
int main() {
int i;
clock_t start, stop;
initPositions();
updatePositions();
start=clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0;
computePot();
if (i%10 == 0) printf("%5d: Potential: %10.7f ", i, pot);
updatePositions();
}
stop=clock();
printf ("Seconds = %10.9f ",(double)(stop-start)/ CLOCKS_PER_SEC);
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[i][j] = 0.5 + ( (double) rand() / (double) RAND_MAX );
}
void updatePositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[i][j] -= 0.5 + ( (double) rand() / (double) RAND_MAX );
}
int computePot() {
int i, j;
for( i=0; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
distx = pow( (r[0][j] - r[0][i]), 2 );
disty = pow( (r[1][j] - r[1][i]), 2 );
distz = pow( (r[2][j] - r[2][i]), 2 );
dist = sqrt( distx + disty + distz );
pot += 1.0 / dist;
}
}
return 0;
}
代码执行时间 3.97秒
整个代码的代码量还是很小的，而且很容易发现computePot是一个O(N*N)的算法，运行时间上它的消耗最大，
占了总运行时间的95%以上(建议使用工具来定位程序的热点)；
对于computePot函数的优化，最容易看到的优化点就是pow函数的调用，pow(x,y)函数是一个的很耗时的操
作(很多库，会针对y为整数做优化)，把pow(x,y)改成x*x的方式；
修改后的computePot代码如下(其它代码不变)：
inline pow2(const double x) { return x*x; }
int computePot() {
int i, j;
for( i=0; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
distx = pow2( (r[0][j] - r[0][i]) );
disty = pow2( (r[1][j] - r[1][i]) );
distz = pow2( (r[2][j] - r[2][i]) );
dist = sqrt( distx + disty + distz );
pot += 1.0 / dist;
}
}
return 0;
}
代码执行时间 2.50秒 该修改加速比：158.8% 相对于原始代码加速比：158.8%
由于computePot中对数组r的访问顺序和r定义的顺序方式不一致，尝试修改r的数据组织结构;
将double r[DIMS][NPARTS] 的定义改为 double r[NPARTS][DIMS];
（后面证明这个修改使自己走了很大的弯路;原来的向量格式在使用SSE2优化的时候，代码上会稍微简单一点）
也相应的修改了其他访问r的代码:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <time.h>
#define NPARTS 1000
#define NITER 201
#define DIMS 3
int rand( void );
int computePot(void);
void initPositions(void);
void updatePositions(void);
double r[NPARTS][DIMS];
double pot;
double distx, disty, distz, dist;
int main() {
int i;
clock_t start, stop;
initPositions();
updatePositions();
start=clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0;
computePot();
if (i%10 == 0) printf("%5d: Potential: %10.7f ", i, pot);
updatePositions();
}
stop=clock();
printf ("Seconds = %10.9f ",(double)(stop-start)/ CLOCKS_PER_SEC);
getchar();
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[j][i] = 0.5 + ( (double) rand() / (double) RAND_MAX );
}
void updatePositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[j][i] -= 0.5 + ( (double) rand() / (double) RAND_MAX );
}
inline double pow2(const double x) { return x*x; }
int computePot() {
int i, j;
for( i=0; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
distx = pow2( (r[j][0] - r[i][0]) );
disty = pow2( (r[j][1] - r[i][1]) );
distz = pow2( (r[j][2] - r[i][2]) );
dist = sqrt( distx + disty + distz );
pot += 1.0 / dist;
}
}
return 0;
}
代码执行时间 2.39秒 该修改加速比：104.6% 相对于原始代码加速比：166.1%
r可以看作一个3维空间点的数组，computePot函数就是求这些点之间距离的某种势能:)
为了将每个点的内存访问对齐到16byte边界，
将double r[NPARTS][DIMS] 的定义改为 __declspec(align(16))doubler[NPARTS][DIMS+1]；
(其它代码不变)(速度略有下降，为使用SSE2优化作准备)
代码执行时间 2.42秒 该修改加速比： 98.8% 相对于原始代码加速比：164.0%
使用SSE2来优化computePot函数；SSE2是一种单指令流多数据流指令集，里面包含有能同时处理两个64bit浮点数的指令；
代码需要 #include <emmintrin.h> 文件，从而利用编译器来优化SSE2指令，就不用手工编写汇编代码了；
修改后的computePot代码如下(其它代码不变)：
int computePot() {
int i, j;
__m128d lcpot;
lcpot=_mm_setzero_pd();
for( i=0; i<NPARTS; i++ ) {
__m128d posi_l=*(__m128d*)(&r[i][0]);
__m128d posi_h=*(__m128d*)(&r[i][2]);
j=0;
for(;j+1<i;++j)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
__m128d dm0=_mm_add_pd(dm0_l,dm0_h);
__m128d dm1;
dm1=_mm_unpackhi_pd(dm0,dm0);
dm0=_mm_add_sd(dm0,dm1);
dm1=_mm_sqrt_sd(dm0,dm0);
dm0=_mm_div_sd(dm1,dm0);
lcpot=_mm_add_sd(lcpot,dm0);
}
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
pot+=*(double*)(&lcpot);
return 0;
}
代码执行时间 2.31秒 该修改加速比：104.8% 相对于原始代码加速比：171.9%
(时间提高不多，主要时间耗在了_mm_sqrt_sd、_mm_div_sd指令上，所以SSE2的优化效果没有体现出来)
SSE3中有一条_mm_hadd_pd指令(水平加)，可以用来改善一点点代码，
程序需要#include <intrin.h>文件：
int computePot() {
int i, j;
__m128d lcpot;
lcpot=_mm_setzero_pd();
for( i=0; i<NPARTS; i++ ) {
__m128d posi_l=*(__m128d*)(&r[i][0]);
__m128d posi_h=*(__m128d*)(&r[i][2]);
j=0;
for(;j+1<i;++j)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
__m128d dm0=_mm_add_pd(dm0_l,dm0_h);
__m128d dm1;
//dm1=_mm_unpackhi_pd(dm0,dm0);
//dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_hadd_pd(dm0,dm0); //SSE3
dm1=_mm_sqrt_sd(dm0,dm0);
dm0=_mm_div_sd(dm1,dm0);
lcpot=_mm_add_sd(lcpot,dm0);
}
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
pot+=*(double*)(&lcpot);
return 0;
}
代码执行时间 2.22秒 该修改加速比：104.1% 相对于原始代码加速比：178.8%
线程化computePot函数来达到并行使用两个CPU的目的；
也就是将computePot的内部运算改写为支持并行运算的形式，并合理划分成多个任务交给线程执行；
程序实现了一个比较通用的工作线程池，来利用多核CPU完成
代码需要 #include <vector>#include <process.h>
并在编译器中设定项目的属性：运行时库 为 多线程(/MT)
直接给出全部实现代码:代码执行时间 1.16秒 该修改加速比：191.4% 相对于原始代码加速比：342.2%
并行代码使程序性能得到了巨大的提高!!
//* compute the potential energy of a collection of */
//* particles interacting via pairwise potential */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <time.h>
#include <emmintrin.h> //使用SSE2
#include <intrin.h> //使用SSE3
#include <vector>
#include <process.h> //使用线程
#define NPARTS 1000
#define NITER 201
#define DIMS 3
const int DefthreadCount=2; //设为1则为单线程执行，如果有4个核可以设为4:)
int rand( void );
int computePot(void);
void initPositions(void);
void updatePositions(void);
__declspec(align(16))double r[NPARTS][DIMS+1];
double pot;
double distx, disty, distz, dist;
int main() {
int i;
clock_t start, stop;
initPositions();
updatePositions();
start=clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0;
computePot();
if (i%10 == 0) printf("%5d: Potential: %10.7f ", i, pot);
updatePositions();
}
stop=clock();
printf ("Seconds = %10.9f ",(double)(stop-start)/ CLOCKS_PER_SEC);
getchar();
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[j][i] = 0.5 + ( (double) rand() / (double) RAND_MAX );
}
void updatePositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[j][i] -= 0.5 + ( (double) rand() / (double) RAND_MAX );
}
struct TWorkData
{
long iBegin;
long iEnd;
double fResult;
};
void computePotPart(TWorkData* work_data) {
int i, j;
__m128d lcpot;
lcpot=_mm_setzero_pd();
for( i=work_data->iBegin; i<work_data->iEnd; ++i ) {
__m128d posi_l=*(__m128d*)(&r[i][0]);
__m128d posi_h=*(__m128d*)(&r[i][2]);
j=0;
for(;j+1<i;++j)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
__m128d dm0=_mm_add_pd(dm0_l,dm0_h);
__m128d dm1;
//dm1=_mm_unpackhi_pd(dm0,dm0);
//dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_hadd_pd(dm0,dm0); //SSE3
dm1=_mm_sqrt_sd(dm0,dm0);
dm0=_mm_div_sd(dm1,dm0);
lcpot=_mm_add_sd(lcpot,dm0);
}
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
work_data->fResult=*(double*)(&lcpot);
}
/
//工作线程池 TWorkThreadPool
//用于把一个任务拆分成多个线程任务
//要求每个小任务任务量相近
//todo:改成任务领取模式
class TWorkThreadPool;
typedef void (*TThreadCallBack)(void * pData);
enum TThreadState{ thrStartup=0, thrReady, thrBusy, thrTerminate, thrDeath };
class TWorkThread
{
public:
volatile HANDLE thread_handle;
volatile enum TThreadState state;
volatile TThreadCallBack func;
volatile void * pdata; //work data
volatile HANDLE waitfor_event;
TWorkThreadPool* pool;
volatile DWORD thread_ThreadAffinityMask;
TWorkThread() { memset(this,0,sizeof(TWorkThread)); }
};
void do_work_end(TWorkThread* thread_data);
void __cdecl thread_dowork(TWorkThread* thread_data) //void __stdcall thread_dowork(TWorkThread* thread_data)
{
volatile TThreadState& state=thread_data->state;
SetThreadAffinityMask(GetCurrentThread(),thread_data->thread_ThreadAffinityMask);
state = thrStartup;
while(true)
{
WaitForSingleObject(thread_data->waitfor_event, -1);
if(state == thrTerminate)
break;
state = thrBusy;
volatile TThreadCallBack& func=thread_data->func;
if (func!=0)
func((void *)thread_data->pdata);
do_work_end(thread_data);
}
state = thrDeath;
_endthread();//ExitThread(0);
}
class TWorkThreadPool
{
private:
volatile HANDLE thread_event;
volatile HANDLE new_thread_event;
std::vector<TWorkThread> work_threads;
inline int passel_count() const { return (int)work_threads.size()+1; }
void inti_threads(long work_count) {
SYSTEM_INFO SystemInfo;
GetSystemInfo(&SystemInfo);
long cpu_count =SystemInfo.dwNumberOfProcessors;
long best_count =cpu_count;
if (cpu_count>work_count) best_count=work_count;
long newthrcount=best_count - 1;
work_threads.resize(newthrcount);
thread_event = CreateSemaphore(0, 0,newthrcount , 0);
new_thread_event = CreateSemaphore(0, 0,newthrcount , 0);
for(long i = 0; i < newthrcount; ++i)
{
work_threads[i].waitfor_event=thread_event;
work_threads[i].state = thrTerminate;
work_threads[i].pool=this;
work_threads[i].thread_ThreadAffinityMask=1<<(i+1);
//DWORD thr_id;
work_threads[i].thread_handle =(HANDLE)_beginthread((void (__cdecl *)(void *))thread_dowork, 0, (void*)&work_threads[i]);
//CreateThread(0, 0, (LPTHREAD_START_ROUTINE)thread_dowork,(void*) &work_threads[i], 0, &thr_id);
}
SetThreadAffinityMask(GetCurrentThread(),0x01);
for(long i = 0; i < newthrcount; ++i)
{
while(true) {
if (work_threads[i].state == thrStartup) break;
else Sleep(0);
}
work_threads[i].state = thrReady;
}
}
void free_threads(void)
{
long thr_count=(long)work_threads.size();
long i;
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrReady) break;
else Sleep(0);
}
work_threads[i].state=thrTerminate;
}
if (thr_count>0)
ReleaseSemaphore(thread_event,thr_count, 0);
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrDeath) break;
else Sleep(0);
}
}
CloseHandle(thread_event);
CloseHandle(new_thread_event);
work_threads.clear();
}
void passel_work(TThreadCallBack work_proc,void** word_data_list,int work_count) {
//assert(work_count>=1);
//assert(work_count<=passel_count());
if (work_count==1)
{
work_proc(word_data_list[0]);
}
else
{
long i;
long thr_count=(long)work_threads.size();
for(i = 0; i < work_count-1; ++i)
{
work_threads[i].func = work_proc;
work_threads[i].pdata =word_data_list[i];
work_threads[i].state = thrBusy;
}
for(i = work_count-1; i < thr_count; ++i)
{
work_threads[i].func = 0;
work_threads[i].pdata =0;
work_threads[i].state = thrBusy;
}
if (thr_count>0)
ReleaseSemaphore(thread_event,thr_count, 0);
//current thread do a work
work_proc(word_data_list[work_count-1]);
//wait for work finish
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrReady) break;
else Sleep(0);
}
}
std::swap(thread_event,new_thread_event);
}
}
public:
explicit TWorkThreadPool(unsigned long work_count):thread_event(0),work_threads() {
//assert(work_count>=1);
inti_threads(work_count); }
~TWorkThreadPool() { free_threads(); }
long best_work_count() const { return passel_count(); }
void work_execute(TThreadCallBack work_proc,void** word_data_list,int work_count) {
while (work_count>0)
{
long passel_work_count;
if (work_count>=passel_count())
passel_work_count=passel_count();
else
passel_work_count=work_count;
passel_work(work_proc,word_data_list,passel_work_count);
work_count-=passel_work_count;
word_data_list=&word_data_list[passel_work_count];
}
}
inline void DoWorkEnd(TWorkThread* thread_data){
thread_data->waitfor_event=new_thread_event;
thread_data->func=0;
thread_data->state = thrReady;
}
};
void do_work_end(TWorkThread* thread_data)
{
thread_data->pool->DoWorkEnd(thread_data);
}
//TWorkThreadPool end;
int computePot() {
static bool is_inti_work=false;
static TWorkThreadPool g_work_thread_pool(DefthreadCount);//线程池
static TWorkData work_list[DefthreadCount];
static TWorkData* pwork_list[DefthreadCount];
int i;
if (!is_inti_work)
{
long fi=0;
for (int i=0;i<DefthreadCount;++i)
{
if (0==i)
work_list[i].iBegin=0;
else
work_list[i].iBegin=work_list[i-1].iEnd;
if (i==DefthreadCount-1)
work_list[i].iEnd=(long)NPARTS;
else
work_list[i].iEnd=(long)( (double)(NPARTS-1)*sqrt((double)(i+1)/DefthreadCount)+1+0.5 );//todo:更准确的任务分配方式
pwork_list[i]=&work_list[i];
}
is_inti_work=true;
}
g_work_thread_pool.work_execute((TThreadCallBack)computePotPart,(void **)(&pwork_list[0]),DefthreadCount);
for (i=0;i<DefthreadCount;++i)
pot+=work_list[i].fResult;
return 0;
}
在computePot函数中,开方(Sqrt)和倒数(Div)占用了大部分的计算时间;其实把这两个操作和并在一起的操作(invSqrt)反而有很快的计算方法;
这就是牛顿迭代法,牛顿迭代法原理、invSqrt迭代公式的推导过程、收敛性证明 可以参看我的一篇blog文章:
《代码优化－之－优化除法》： http://blog.csdn.net/housisong/archive/2006/08/25/1116423.aspx
invSqrt(a)=1/sqrt(a) 的牛顿迭代公式: x_next=(3-a*x*x)*x*0.5;
先给出一个invSqrt(a)的近似解然后把它带入迭代式，就可以得到更精确的解；
产生invSqrt(a)的一个近似解的方法有很多，SSE中提供了一条_mm_rsqrt_ps指令用来求invSqrt(a)的近似解，可以把它的结果作为一个初始值，
然后进行多次迭代，从而达到自己想要的解的精度，实现代码如下（其它代码不变）:代码执行时间 1.05秒 该修改加速比：110.5% 相对于原始代码加速比：378.1%
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5) };
const __m128d xmmd1_5={ (1.5),(1.5) };
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5) };
const __m128d xmmd_0_5={ (-0.5),(-0.5) };
void computePotPart(TWorkData* work_data) {
int i, j;
__m128d lcpot;
lcpot=_mm_setzero_pd();
for( i=work_data->iBegin; i<work_data->iEnd; ++i ) {
__m128d posi_l=*(__m128d*)(&r[i][0]);
__m128d posi_h=*(__m128d*)(&r[i][2]);
j=0;
for(;j+1<i;++j)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
__m128d dm0=_mm_add_pd(dm0_l,dm0_h);
__m128d dm1;
//dm1=_mm_unpackhi_pd(dm0,dm0);
//dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_hadd_pd(dm0,dm0); //SSE3
//dm1=_mm_sqrt_sd(dm0,dm0);
//dm0=_mm_div_sd(dm1,dm0);
//用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128 sm0=_mm_cvtpd_ps(dm0);
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
dm0=_mm_cvtps_pd(sm0);
//牛顿迭代，提高开方精度
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
//再次迭代，加倍精度 这样和原表达式比就几乎没有误差了
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
}
lcpot=_mm_add_sd(lcpot,dm0);
}
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
work_data->fResult=*(double*)(&lcpot);
}
对computePot的内部j循环作4路展开，从而可以充分利用SSE的计算带宽，代码如下：
(更高路(比如8路)的展开的效果沒有测试过，应该会有利于指令的并行发射和执行)代码执行时间 0.70秒 该修改加速比：150.0% 相对于原始代码加速比：567.1%
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5) };
const __m128d xmmd1_5={ (1.5),(1.5) };
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5) };
const __m128d xmmd_0_5={ (-0.5),(-0.5) };
void computePotPart(TWorkData* work_data) {
int i, j;
__m128d lcpot;
lcpot=_mm_setzero_pd();
for( i=work_data->iBegin; i<work_data->iEnd; ++i ) {
__m128d posi_l=*(__m128d*)(&r[i][0]);
__m128d posi_h=*(__m128d*)(&r[i][2]);
j=0;
//* 把这个循环做4次展开
for(;j+4<i;j+=4)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j ][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j ][2]),posi_h);
__m128d dm1_l=_mm_sub_pd(*(__m128d*)(&r[j+1][0]),posi_l);
__m128d dm1_h=_mm_sub_pd(*(__m128d*)(&r[j+1][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
dm1_l=_mm_mul_pd(dm1_l,dm1_l);
dm1_h=_mm_mul_pd(dm1_h,dm1_h);
__m128d dm0,dm1,dm2;
dm1=_mm_add_pd(dm0_l,dm0_h);
dm2=_mm_add_pd(dm1_l,dm1_h);
dm0=_mm_hadd_pd(dm1,dm2); //SSE3
dm0_l=_mm_sub_pd(*(__m128d*)(&r[j+2][0]),posi_l);
dm0_h=_mm_sub_pd(*(__m128d*)(&r[j+2][2]),posi_h);
dm1_l=_mm_sub_pd(*(__m128d*)(&r[j+3][0]),posi_l);
dm1_h=_mm_sub_pd(*(__m128d*)(&r[j+3][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
dm1_l=_mm_mul_pd(dm1_l,dm1_l);
dm1_h=_mm_mul_pd(dm1_h,dm1_h);
__m128d dm5;
dm1=_mm_add_pd(dm0_l,dm0_h);
dm2=_mm_add_pd(dm1_l,dm1_h);
dm5=_mm_hadd_pd(dm1,dm2); //SSE3
//用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
{
__m128 sm0=_mm_cvtpd_ps(dm0);
__m128 sm1=_mm_cvtpd_ps(dm5);
sm0=_mm_movelh_ps(sm0,sm1);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); //近似计算1/sqrt(a)
//牛顿迭代，提高开方精度
sm1=_mm_mul_ps(sm0,sm0);
sm1=_mm_mul_ps(sm1,sma);
sm1=_mm_add_ps(sm1,xmms1_5);
sm0=_mm_mul_ps(sm0,sm1);
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb=_mm_mul_pd(dm5,xmmd_0_5);
sm1=_mm_movehl_ps(sm0,sm0);
dm0=_mm_cvtps_pd(sm0); //
dm5=_mm_cvtps_pd(sm1); //
//再次迭代，加倍精度
dm1=_mm_mul_pd(dm0,dm0);
dm1=_mm_mul_pd(dm1,dma);
dm1=_mm_add_pd(dm1,xmmd1_5);
dm0=_mm_mul_pd(dm0,dm1);
//
dm2=_mm_mul_pd(dm5,dm5);
dm2=_mm_mul_pd(dm2,dmb);
dm2=_mm_add_pd(dm2,xmmd1_5);
dm5=_mm_mul_pd(dm5,dm2);
}
dm0=_mm_add_pd(dm0,dm5);
lcpot=_mm_add_pd(lcpot,dm0);
}//*/
for(;j+1<i;++j)
{
__m128d dm0_l=_mm_sub_pd(*(__m128d*)(&r[j][0]),posi_l);
__m128d dm0_h=_mm_sub_pd(*(__m128d*)(&r[j][2]),posi_h);
dm0_l=_mm_mul_pd(dm0_l,dm0_l);
dm0_h=_mm_mul_pd(dm0_h,dm0_h);
__m128d dm0=_mm_add_pd(dm0_l,dm0_h);
__m128d dm1;
//dm0=_mm_hadd_pd(dm0,dm0); //SSE3
dm1=_mm_unpackhi_pd(dm0,dm0);
dm0=_mm_add_sd(dm0,dm1);
//dm1=_mm_sqrt_sd(dm0,dm0);
//dm0=_mm_div_sd(dm1,dm0);
//用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
{
__m128 sm0=_mm_cvtpd_ps(dm0);
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
dm0=_mm_cvtps_pd(sm0);
//牛顿迭代，提高开方精度
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
//再次迭代，加倍精度 这样和原表达式比就几乎没有误差了
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
}
lcpot=_mm_add_sd(lcpot,dm0);
}
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
work_data->fResult=*(double*)(&lcpot);
}
将上面的computePotPart函数用汇编实现:代码执行时间 0.59秒 该修改加速比：118.6% 相对于原始代码加速比：672.9%
#define _IS_CPU_SSE3_ACTIVE
//取消该标志则使用不使用SSE3
#define asm __asm
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5) };
const __m128d xmmd1_5={ (1.5),(1.5) };
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5) };
const __m128d xmmd_0_5={ (-0.5),(-0.5) };
void __declspec(naked) computePotPart(TWorkData* work_data) {
//work_data: [ebp+8]
asm
{
push ebp
mov ebp, esp
and esp, -16 //16byte 对齐
push edi
sub esp, 4+8+16 //for a temp_m128d
mov edi, [ebp+8]//work_data
mov ecx, [edi]TWorkData.iBegin
mov edx, [edi]TWorkData.iEnd
cmp ecx, edx
pxor xmm7, xmm7 //fResult
jge $for_end
mov [esp+16], esi //save esi
mov eax, ecx
mov [esp+20], ebx //save ebx
shl eax, 5
$for_begin:
movapd xmm5, [r+eax]
movapd xmm6, [r+eax+16]
xor ebx, ebx
cmp ecx, 4
jle $border
mov esi, 4
mov edi, offset r
align 4
$for4_begin:
movapd xmm0, [edi]
movapd xmm4, [edi+16]
movapd xmm2, [edi+32]
movapd xmm3, [edi+48]
movapd xmm1, [edi+96]
subpd xmm0, xmm5
mulpd xmm0, xmm0
subpd xmm4, xmm6
mulpd xmm4, xmm4
addpd xmm0, xmm4
movapd xmm4, [edi+64]
subpd xmm2, xmm5
mulpd xmm2, xmm2
subpd xmm3, xmm6
mulpd xmm3, xmm3
addpd xmm2, xmm3
movapd xmm3, [edi+80]
subpd xmm4, xmm5
mulpd xmm4, xmm4
#ifdef _IS_CPU_SSE3_ACTIVE
//SSE3
haddpd xmm0, xmm2
#else
movapd [esp],xmm1 //temp_m128d
movapd xmm1, xmm0
unpckhpd xmm0, xmm2
unpcklpd xmm1, xmm2
addpd xmm0,xmm1
movapd xmm1,[esp]
#endif
movapd xmm2, [edi+112]
subpd xmm3, xmm6
mulpd xmm3, xmm3
addpd xmm4, xmm3
subpd xmm1, xmm5
mulpd xmm1, xmm1
subpd xmm2, xmm6
mulpd xmm2, xmm2
addpd xmm1, xmm2
movaps xmm2, xmms_0_5
#ifdef _IS_CPU_SSE3_ACTIVE
//SSE3
haddpd xmm4, xmm1
#else
movapd xmm3, xmm4
unpckhpd xmm4, xmm1
unpcklpd xmm3, xmm1
addpd xmm4, xmm3
#endif
cvtpd2ps xmm1, xmm0
mov ebx, esi
cvtpd2ps xmm3, xmm4
add edi, 4*32
movlhps xmm1, xmm3
mulps xmm2, xmm1
rsqrtps xmm3, xmm1
movaps xmm1, xmm3
mulps xmm1, xmm3
mulps xmm1, xmm2
movapd xmm2, xmmd_0_5
mulpd xmm0, xmm2
mulpd xmm4, xmm2
addps xmm1, xmms1_5
mulps xmm3, xmm1
cvtps2pd xmm1, xmm3
add esi, 4
movaps xmm2, xmm3
movhlps xmm2, xmm3
cvtps2pd xmm3, xmm2
cmp ecx, esi
movapd xmm2, xmm1
mulpd xmm2, xmm1
mulpd xmm2, xmm0
movapd xmm0, xmmd1_5
addpd xmm2, xmm0
mulpd xmm1, xmm2
movapd xmm2, xmm3
mulpd xmm2, xmm3
mulpd xmm2, xmm4
addpd xmm2, xmm0
mulpd xmm3, xmm2
addpd xmm1, xmm3
addpd xmm7, xmm1
jg $for4_begin
$border:
lea esi, [ebx+1]
cmp ecx, esi
jle $for_border_end
movapd xmm4, xmmd1_5
movapd xmm3, xmmd_0_5
$for_border_begin:
shl ebx, 5
movapd xmm2, [r+ebx]
movapd xmm1, [r+ebx+16]
subpd xmm2, xmm5
mulpd xmm2, xmm2
subpd xmm1, xmm6
mulpd xmm1, xmm1
addpd xmm2, xmm1
#ifdef _IS_CPU_SSE3_ACTIVE
//SSE3
haddpd xmm2, xmm2
#else
movapd xmm0, xmm2
unpckhpd xmm0, xmm0
addsd xmm2, xmm0
#endif
cvtpd2ps xmm0, xmm2
mov ebx, esi
add esi, 1
mulsd xmm2, xmm3
rsqrtss xmm0, xmm0
cvtps2pd xmm0, xmm0
cmp ecx, esi
movapd xmm1, xmm0
mulsd xmm1, xmm0
mulsd xmm1, xmm2
addsd xmm1, xmm4
mulsd xmm0, xmm1
movapd xmm1, xmm0
mulsd xmm1, xmm0
mulsd xmm1, xmm2
addsd xmm1, xmm4
mulsd xmm0, xmm1
addsd xmm7, xmm0
jg $for_border_begin
$for_border_end:
add ecx, 1
add eax, 32
cmp ecx, edx
jl $for_begin
mov esi, [esp+16] //resume esi
mov ebx, [esp+20] //resume ebx
$for_end:
#ifdef _IS_CPU_SSE3_ACTIVE
//SSE3
haddpd xmm7, xmm7
#else
movapd xmm0, xmm7
unpckhpd xmm0, xmm0
addpd xmm7, xmm0
#endif
mov edi,[ebp+8]//work_data
movsd qword ptr [edi]TWorkData.fResult, xmm7
add esp, 16+4+8
pop edi
mov esp, ebp
pop ebp
ret
}
}
由于修改了r数组的数据组织方式，虽然有利于内存访问，但也由此带来了一定的代码复杂度；
但这个程序的数据量本来就不大，所以尝试保留原来的向量化数据，然后代码也会简单一些，
并不再使用SSE3,所以有了新的代码，保留double r[DIMS][NPARTS]的定义，并且对整个程序能
够加快速度的其他地方都进行了一些改进。
代码执行时间 0.54秒 该修改加速比：109.3% 相对于原始代码加速比：735.2%
完整代码如下：
/* compute the potential energy of a collection of */
/* particles interacting via pairwise potential */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <xmmintrin.h>
#include <process.h>
#include <vector>
//#include <iostream>
#include <emmintrin.h> //使用SSE2
//作者：侯思松 HouSisong@263.net
//计算结果的精度为double双精度浮点 版本
//#define _MY_TIME_CLOCK
#ifdef _MY_TIME_CLOCK
#define clock_t __int64
clock_t CLOCKS_PER_SEC=0;
inline clock_t clock() {
__int64 result;
if (CLOCKS_PER_SEC==0)
{
QueryPerformanceFrequency((LARGE_INTEGER *)&result);
CLOCKS_PER_SEC=(clock_t)result;
}
QueryPerformanceCounter((LARGE_INTEGER *)&result);
return (clock_t)result;
}
#else
#include <time.h>
#endif
#define _IS_USES_MY_RAND
//单线程执行rand函数，所以使用自定义rand是安全的
const long DefthreadCount=2; //1,2,4,8,16,.. 把计算任务分成多个任务并行执行
inline double& m128d_value(__m128d& x,const long index) { return ((double*)(&x))[index]; }
#define NPARTS 1000
#define NITER 201
#define DIMS 3
#ifdef _IS_USES_MY_RAND
class CMyRand
{
private:
unsigned long _my_holdrand;
public:
CMyRand():_my_holdrand(1){}
inline int _my_rand (void)
{
unsigned long result=_my_holdrand * 214013 + 2531011;
_my_holdrand =result;
return ( (result>>16) & RAND_MAX );
}
};
CMyRand _MyRand;
inline int _my_rand (void){ return _MyRand._my_rand(); }
#else
#define _my_rand rand
#endif
int rand( void );
int computePot();
void initPositions(void);
void updatePositions(void);
__declspec(align(16)) double r[DIMS][(NPARTS+1)/2*2]; //16byte对齐
enum{ vector_byte_size=sizeof(r)/DIMS };
double pot;
int main() {
int i;
clock_t start, stop;
//char ctmp; std::cin>>ctmp;
initPositions();
updatePositions();
start=clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0;
computePot();
if (i%10 == 0) printf("%5d: Potential: %20.7f", i, pot);
updatePositions();
}
pot = 0.0;
stop=clock();
printf ("Seconds = %10.9f",(double)(stop-start)/ CLOCKS_PER_SEC);
return 0;
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ ){
for( j=0; j<NPARTS; ++j )
{
r[i][j]=( 0.5 + _my_rand() * (1.0/RAND_MAX) );
}
}
}
void updatePositions() {
int i,j;
for (i=0;i<DIMS;++i)
{
for( j=0; j<NPARTS; ++j )
{
r[i][j] -=( 0.5 + _my_rand() *(1.0/RAND_MAX) );
}
}
}
struct TWorkData
{
long iBegin;
long iEnd;
double fResult;
};
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5) };
const __m128d xmmd_0_5={ (-0.5),(-0.5) };
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5) };
const __m128d xmmd1_5={ (1.5),(1.5) };
void computePotPart_forj(int i,__m128d* pResult)
{
__m128d lcpot=_mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0][i]);
__m128d _mmi1=_mm_set1_pd(-r[1][i]);
__m128d _mmi2=_mm_set1_pd(-r[2][i]);
int j=0;
//*
for(;j+4<i;j+=4)
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0][j],_mmi0);
dm0=_mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1][j],_mmi1);
dm1=_mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2][j],_mmi2);
dm2=_mm_mul_pd(dm2,dm2);
dm0=_mm_add_pd(dm0,dm1);
dm0=_mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2],_mmi0);
dm5=_mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2],_mmi1);
dm6=_mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2],_mmi2);
dm2=_mm_mul_pd(dm2,dm2);
dm5=_mm_add_pd(dm5,dm6);
dm5=_mm_add_pd(dm5,dm2);
//用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128 sm0=_mm_cvtpd_ps(dm0);
__m128 sm1=_mm_cvtpd_ps(dm5);
sm0=_mm_movelh_ps(sm0,sm1);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); //计算1/sqrt(a)
//牛顿迭代，提高开方精度
sma=_mm_mul_ps(sma,sm0);
sma=_mm_mul_ps(sma,sm0);
sma=_mm_add_ps(sma,xmms1_5);
sm0=_mm_mul_ps(sm0,sma);
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb=_mm_mul_pd(dm5,xmmd_0_5);
sm1=_mm_movehl_ps(sm1,sm0);
dm0=_mm_cvtps_pd(sm0); //
dm5=_mm_cvtps_pd(sm1); //
//再次迭代，加倍精度 这样和原表达式比就几乎没有任何误差了
dma=_mm_mul_pd(dma,dm0);
dmb=_mm_mul_pd(dmb,dm5);
dma=_mm_mul_pd(dma,dm0);
dmb=_mm_mul_pd(dmb,dm5);
dma=_mm_add_pd(dma,xmmd1_5);
dmb=_mm_add_pd(dmb,xmmd1_5);
dm0=_mm_mul_pd(dm0,dma);
dm5=_mm_mul_pd(dm5,dmb);
}
lcpot=_mm_add_pd(lcpot,dm0);
lcpot=_mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++j)
{
__m128d dm0=_mm_set_sd(r[0][j]);
dm0=_mm_add_pd(dm0,_mmi0);
dm0=_mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1][j]);
dm1=_mm_add_sd(dm1,_mmi1);
dm1=_mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2][j]);
dm2=_mm_add_sd(dm2,_mmi2);
dm2=_mm_mul_sd(dm2,dm2);
dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_add_sd(dm0,dm2);
{
__m128 sm0=_mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
//再次迭代
dma=_mm_mul_sd(dma,dm0);
dma=_mm_mul_sd(dma,dm0);
dma=_mm_add_sd(dma,xmmd1_5);
dm0=_mm_mul_sd(dm0,dma);
}
lcpot=_mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(*pResult,lcpot);
}
void computePotPart(TWorkData* work_data) {
int i;
__m128d lcpot=_mm_setzero_pd();
for( i=work_data->iBegin; i<work_data->iEnd; ++i ) {
computePotPart_forj(i,&lcpot);
}
__m128d dm1;
dm1=_mm_unpackhi_pd(lcpot,lcpot);
lcpot=_mm_add_sd(lcpot,dm1);
work_data->fResult=m128d_value(lcpot,0);
}
/
//工作线程池 TWorkThreadPool
//用于把一个任务拆分成多个线程任务
//要求每个小任务任务量相近
//todo:改成任务领取模式
class TWorkThreadPool;
typedef void (*TThreadCallBack)(void * pData);
enum TThreadState{ thrStartup=0, thrReady, thrBusy, thrTerminate, thrDeath };
class TWorkThread
{
public:
volatile HANDLE thread_handle;
volatile enum TThreadState state;
volatile TThreadCallBack func;
volatile void * pdata; //work data
volatile HANDLE waitfor_event;
TWorkThreadPool* pool;
volatile DWORD thread_ThreadAffinityMask;
TWorkThread() { memset(this,0,sizeof(TWorkThread)); }
};
void do_work_end(TWorkThread* thread_data);
void __cdecl thread_dowork(TWorkThread* thread_data) //void __stdcall thread_dowork(TWorkThread* thread_data)
{
volatile TThreadState& state=thread_data->state;
//SetThreadAffinityMask(GetCurrentThread(),thread_data->thread_ThreadAffinityMask);
state = thrStartup;
while(true)
{
WaitForSingleObject(thread_data->waitfor_event, -1);
if(state == thrTerminate)
break;
state = thrBusy;
volatile TThreadCallBack& func=thread_data->func;
if (func!=0)
func((void *)thread_data->pdata);
do_work_end(thread_data);
}
state = thrDeath;
_endthread();//ExitThread(0);
}
class TWorkThreadPool
{
private:
volatile HANDLE thread_event;
volatile HANDLE new_thread_event;
std::vector<TWorkThread> work_threads;
inline int passel_count() const { return (int)work_threads.size()+1; }
void inti_threads(long work_count) {
SYSTEM_INFO SystemInfo;
GetSystemInfo(&SystemInfo);
long cpu_count =SystemInfo.dwNumberOfProcessors;
long best_count =cpu_count;
if (cpu_count>work_count) best_count=work_count;
long newthrcount=best_count - 1;
work_threads.resize(newthrcount);
thread_event = CreateSemaphore(0, 0,newthrcount , 0);
new_thread_event = CreateSemaphore(0, 0,newthrcount , 0);
for(long i = 0; i < newthrcount; ++i)
{
work_threads[i].waitfor_event=thread_event;
work_threads[i].state = thrTerminate;
work_threads[i].pool=this;
work_threads[i].thread_ThreadAffinityMask=1<<(i+1);
//DWORD thr_id;
work_threads[i].thread_handle =(HANDLE)_beginthread((void (__cdecl *)(void *))thread_dowork, 0, (void*)&work_threads[i]);
//CreateThread(0, 0, (LPTHREAD_START_ROUTINE)thread_dowork,(void*) &work_threads[i], 0, &thr_id);
}
//SetThreadAffinityMask(GetCurrentThread(),0x01);
for(long i = 0; i < newthrcount; ++i)
{
while(true) {
if (work_threads[i].state == thrStartup) break;
else Sleep(0);
}
work_threads[i].state = thrReady;
}
}
void free_threads(void)
{
long thr_count=(long)work_threads.size();
long i;
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrReady) break;
else Sleep(0);
}
work_threads[i].state=thrTerminate;
}
if (thr_count>0)
ReleaseSemaphore(thread_event,thr_count, 0);
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrDeath) break;
else Sleep(0);
}
}
CloseHandle(thread_event);
CloseHandle(new_thread_event);
work_threads.clear();
}
void passel_work(TThreadCallBack work_proc,void** word_data_list,int work_count) {
//assert(work_count>=1);
//assert(work_count<=passel_count());
if (work_count==1)
{
work_proc(word_data_list[0]);
}
else
{
long i;
long thr_count=(long)work_threads.size();
for(i = 0; i < work_count-1; ++i)
{
work_threads[i].func = work_proc;
work_threads[i].pdata =word_data_list[i];
work_threads[i].state = thrBusy;
}
for(i = work_count-1; i < thr_count; ++i)
{
work_threads[i].func = 0;
work_threads[i].pdata =0;
work_threads[i].state = thrBusy;
}
if (thr_count>0)
ReleaseSemaphore(thread_event,thr_count, 0);
//current thread do a work
work_proc(word_data_list[work_count-1]);
//wait for work finish
for(i = 0; i <thr_count; ++i)
{
while(true) {
if (work_threads[i].state == thrReady) break;
else Sleep(0);
}
}
std::swap(thread_event,new_thread_event);
}
}
public:
explicit TWorkThreadPool(unsigned long work_count):thread_event(0),work_threads() {
//assert(work_count>=1);
inti_threads(work_count); }
~TWorkThreadPool() { free_threads(); }
long best_work_count() const { return passel_count(); }
void work_execute(TThreadCallBack work_proc,void** word_data_list,int work_count) {
while (work_count>0)
{
long passel_work_count;
if (work_count>=passel_count())
passel_work_count=passel_count();
else
passel_work_count=work_count;
passel_work(work_proc,word_data_list,passel_work_count);
work_count-=passel_work_count;
word_data_list=&word_data_list[passel_work_count];
}
}
inline void DoWorkEnd(TWorkThread* thread_data){
thread_data->waitfor_event=new_thread_event;
thread_data->func=0;
thread_data->state = thrReady;
}
};
void do_work_end(TWorkThread* thread_data)
{
thread_data->pool->DoWorkEnd(thread_data);
}
//TWorkThreadPool end;
static TWorkThreadPool g_work_thread_pool(DefthreadCount);//线程池
int computePot() {
static bool is_inti_work=false;
static TWorkData work_list[DefthreadCount];
static TWorkData* pwork_list[DefthreadCount];
if (!is_inti_work)
{
long fi=0;
for (int i=0;i<DefthreadCount;++i)
{
if (0==i)
work_list[i].iBegin=0;
else
work_list[i].iBegin=work_list[i-1].iEnd;
if (i==DefthreadCount-1)
work_list[i].iEnd=(long)NPARTS;
else
work_list[i].iEnd=(long)( (double)(NPARTS-1)*sqrt((double)(i+1)/DefthreadCount)+1);
pwork_list[i]=&work_list[i];
}
is_inti_work=true;
}
g_work_thread_pool.work_execute((TThreadCallBack)computePotPart,(void **)(&pwork_list[0]),DefthreadCount);
for (long i=0;i<DefthreadCount;++i)
pot+=work_list[i].fResult;
return 0;
}
由于新的实现方案用编译器优化出来的速度就超过了前面的手工汇编的实现，所以自己的推想不错；
接着用汇编实现了computePotPart_forj函数，并细调了汇编码(这方面的经验不多，很可能有这样的风险，就是某个调整在AMD上执行加快，但在酷睿2上执行得更慢:)，代码如下(其它部分代码不变)：
代码执行时间 0.453秒 该修改加速比：119.2% 相对于原始代码加速比：876.4%
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5) };
const __m128d xmmd_0_5={ (-0.5),(-0.5) };
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5) };
const __m128d xmmd1_5={ (1.5),(1.5) };
#define asm __asm
void __declspec(naked) computePotPart_forj(int i,__m128d* pResult)
{
//i -> [ebp+ 8]
//pResult -> [ebp+12]
asm
{
push ebp
mov ebp, esp
and esp, -16 //16byte对齐
sub esp, 16 //[esp]一个临时__m128d空间
mov eax, [ebp+8] //i
movsd xmm1, [r+eax*8+vector_byte_size*2]
mov edx, offset r
pxor xmm0, xmm0
movsd xmm2, [r+eax*8+vector_byte_size]
mov ecx, 4
pxor xmm5, xmm5
movsd xmm3, [r+eax*8]
cmp eax, ecx
pxor xmm6, xmm6
subsd xmm0, xmm1
subsd xmm5, xmm2
subsd xmm6, xmm3
pxor xmm7, xmm7 //result
unpcklpd xmm0, xmm0
unpcklpd xmm5, xmm5
unpcklpd xmm6, xmm6
movapd [esp], xmm0
jle $for4_end
align 4
$for4_begin:
movapd xmm0, [edx]
movapd xmm4, [edx+vector_byte_size]
movapd xmm1, [edx+vector_byte_size*2]
movapd xmm3, [esp]
addpd xmm0, xmm6
movapd xmm2, [edx+vector_byte_size*2+16]
mulpd xmm0, xmm0
addpd xmm4, xmm5
mulpd xmm4, xmm4
addpd xmm1, xmm3
addpd xmm0, xmm4
mulpd xmm1, xmm1
movapd xmm4, [edx+16]
addpd xmm0, xmm1
addpd xmm4, xmm6
movapd xmm1, [edx+vector_byte_size+16]
mulpd xmm4, xmm4
addpd xmm1, xmm5
addpd xmm2, xmm3
mulpd xmm1, xmm1
mulpd xmm2, xmm2
addpd xmm4, xmm1
cvtpd2ps xmm3, xmm0
addpd xmm4, xmm2
movapd xmm2, xmmd_0_5
cvtpd2ps xmm1, xmm4
add ecx, 4
add edx, 4*8
mulpd xmm0, xmm2
mulpd xmm4, xmm2
movaps xmm2, xmms_0_5
movlhps xmm3, xmm1
mulps xmm2, xmm3 //-0.5*a
rsqrtps xmm1, xmm3
mulps xmm2, xmm1
movaps xmm3, xmms1_5
mulps xmm2, xmm1
addps xmm2, xmm3 //+1.5
mulps xmm1, xmm2
movhlps xmm2, xmm1
cvtps2pd xmm3, xmm1
cvtps2pd xmm1, xmm2
mulpd xmm0,xmm3
mulpd xmm4,xmm1
movapd xmm2, xmmd1_5
cmp eax, ecx
mulpd xmm0,xmm3
mulpd xmm4,xmm1
addpd xmm0,xmm2
addpd xmm4,xmm2
mulpd xmm0,xmm3
mulpd xmm4,xmm1
addpd xmm7, xmm0
addpd xmm7, xmm4
jg $for4_begin
$for4_end:
sub ecx, 3
cmp eax, ecx
jle $for1_end
movsd xmm3, [esp]
movapd xmm4, xmmd1_5
$for1_begin:
movsd xmm0, [edx]
movsd xmm2, [edx+vector_byte_size]
addsd xmm0, xmm6
addsd xmm2, xmm5
mulsd xmm0, xmm0
mulsd xmm2, xmm2
addsd xmm0, xmm2
movsd xmm2, [edx+vector_byte_size*2]
inc ecx
addsd xmm2, xmm3
mulsd xmm2, xmm2
movsd xmm1, xmmd_0_5
addsd xmm0, xmm2
cvtpd2ps xmm2, xmm0
add edx, 8
cmp eax, ecx
mulsd xmm0, xmm1
rsqrtss xmm2, xmm2
movapd xmm1,xmm0
cvtps2pd xmm2, xmm2
mulsd xmm0, xmm2
mulsd xmm0, xmm2
addsd xmm0, xmm4
mulsd xmm0, xmm2
mulsd xmm1, xmm0
mulsd xmm1, xmm0
addsd xmm1, xmm4
mulsd xmm0, xmm1
addsd xmm7, xmm0
jg $for1_begin
$for1_end:
mov edx,[ebp+12] //pResult
addpd xmm7, [edx]
movapd [edx], xmm7
add esp, 16
mov esp, ebp
pop ebp
ret
}
}
这个版本也是自己提交的最终版本，相对于原始代码加速比：876.4%
把结果汇总，优化过程说明:
初始代码 3.97s
pow(x,2) 改为 (x*x) 2.50s
将double r[DIMS][NPARTS] 改为 doubler[NPARTS][DIMS] 2.39s
接着改为__declspec(align(16))doubler[NPARTS][DIMS+1]； 2.42s
具体说明：这两步主要是增加内存顺序访问(增加缓存命中)和对齐内存
使用SSE2优化computePot函数 2.31s
具体说明： 引入<emmintrin.h>文件
使用SSE3优化computePot函数 2.22s
具体说明： 引入<intrin.h>文件，使用_mm_hadd_pd等函数指令优化SEE2版本的computePot函数
建立一个工作线程池来并行执行computePot函数 1.16s
具体说明：利用线程池来减少线程的创建和销毁开销
使用_mm_rsqrt_ps开方和牛顿迭代两次来替换SSE2的_mm_sqrt_pd和_mm_div_pd 1.05s
对computePot的内部j循环作4路展开 0.70s
使用线程绑定CPU来减少线程切换 0.71s
具体说明：效果不明显
使用内联汇编来优化computePot函数 0.59s
保留double r[DIMS][NPARTS]的定义 0.54s
具体说明：因为数据量比较小，以前的优化内存顺序访问带来了算法的复杂，
不如保留原数据格式，简化算法的效果好
内联汇编来优化computePot函数 0.453s
这些介绍是我这次优化过程的主线，其实还实现了一些其他的版本：以牺牲计算精度换取速度的好几个中间版本,
和一些不成功的尝试；这些可以参考我的另一篇文章的补充说明(提供一个完整的float实现版本,double到float的手工转换、手工得到invSqrt的粗略起始迭代值等其它几个不算成功的实现的简要介绍，也许还有那么点参考价值)
96
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
