Mic*_*ler 7 x86 fortran simd vectorization avx512
我对于在分支方面理论上掩蔽可以做什么感到困惑.假设我有一个Skylake-SP(哈,我希望......),我们忽略了编译器功能,理论上可能的是:
如果分支条件依赖于静态标志,并且所有分支都将数组设置为计算结果,假设编译器不将其优化为两个单独的循环,它是否可以向量化？
do i = 1, nx
if (my_flag .eq. 0) then
a(i) = b(i) ** 2
else
a(i) = b(i) ** 3
end if
end do
如果仅作为分支的子集设置有问题的值,它可以矢量化吗？
do i = 1, nx
if (my_flag .eq. 0) then
a(i) = b(i) ** 2
end if
end do
如果分支条件本身依赖于矢量数据,它可以矢量化吗？
do i = 1, nx
if (c(i) > 0) then
a(i) = b(i) ** 2
else
a(i) = b(i) ** 3
end if
end do
是的,高效的ASM实现是可能的任何SSE2/SSE4.1(为blendps)/ AVX/AVX-512,为您所有的循环,和编译器做自动向量化的做法,但gcc7.2/clang5.0/ICC18都错过了优化.
根据对SKYLAKE微架构-AVX512(见下文),静态分析的有效展开实现最终环路可以(取决于你有多少解开加循环开销)以每1.25个时钟周期的结果一个64字节的矢量运行.实际上,如果您的数据在L1D高速缓存中很热,则每个向量的1.33或1.5个时钟周期可能是可实现的.否则你很容易在L2带宽上遇到瓶颈,因为每个存储矢量64B存储加载2x 64B.
对于你的循环的C版本,gcc,clang和ICC都会像我手工做的那样自动向量化:请参阅Godbolt编译器资源管理器中的 source + asm .
我不得不使用-ffast-mathgcc来自动矢量化.IDK为什么没有意识到它可以安全地自动矢量化而不破坏严格的FP规则.
Clang似乎正在评估tmp*tmp和tmp*tmp*tmp分开,并将这两个结果混合而不是有条件地进行第二次乘法.
gcc两者都相乘并使用单独的movaps来合并另一种方式,因为它没有弄清楚如何反转条件.
ICC用于KNOTW反转条件,但是然后第二次使用合并屏蔽,就像我一样.
更改代码以在分支而不是分支中进行额外的乘法(**3而不是**2)使得所有3个编译器生成更好的代码,ifelse而不会从其他方式分支错误优化.(gcc仍然缺少优化,但ICC和clang看起来很稳定,两者基本上都是我手写代码所做的事情.)
ICC选择仅使用256b向量自动向量化.也许它默认情况下是为了避免降低最大涡轮时钟速度？也许有一个选项可以使用全宽矢量？gcc 8.0快照也可以,但gcc7.2使用ZMM向量.
AVX-512掩码寄存器和合并掩码使其更加高效,但是长时间使用SIMD(甚至非SIMD无分支代码)进行两种方式然后进行混合.例如,为了基于矢量比较结果有条件地添加,使用该矢量比较结果作为AND掩码以保持一些元素不变,并使其他元素为零.
0是附加的身份:x + 0 = x.x + (y&mask)如果掩码是全零,那么是无操作,或者x+y如果掩码是全1 ,则是无操作.请参阅如何在内在函数中使用if条件.(有趣的技巧:使用打包比较结果作为整数-1或0,因此您可以计算匹配但减去比较掩码).
乘法不太简单,因为1是乘法身份,但你可以通过混合来解决这个问题.
假设编译器没有将它优化为两个单独的循环,它可以矢量化吗？
在第一种情况下,如果它没有将条件提升出循环并进行两次循环,那么您应该对编译器不满意.特别是在第二种情况下,它只需要一个循环,因为如果条件为假,则不修改数组.
我们来谈谈第三种情况,因为它只是编译器不应该只提升条件的情况.(如果你的编译器感觉很愚蠢,它可以使用这个版本的循环不变掩码全0或全部为其他版本).
if (c(i) > 0)
所以我们需要加载一个元素向量c并与零进行比较.AVX512可以为16个单精度的向量执行此float操作,其中一个指令具有屏蔽寄存器目标和存储器源操作数.
; with zmm0 = 0.0 in all elements, from vxorps xmm0,xmm0,xmm0 outside the loop.
vcmpps k1, zmm0, [rdx], _CMP_NLT_UQ ; !(0 < c(i))
我知道(从已经写下一部分开始)我想要k1对于c(i) > 0条件为假的元素是真的.只有第二个向量操作数可以是内存而不是寄存器,所以我不得不将其反转并使用不小于而不是大于.(我不能随便用>=的,而不是<,因为那会把无序的情况下,(一个或两个NAN)在错误的类别FP比较有4个可能的结果:上/下/等于/无序的,所以你必须选择一个对于所有4种情况,谓词都可以执行您想要的操作(即源代码所说的,如果您是编译器).如果使用编译-ffast-math器,则允许编译器忽略NaN的可能性.
如果需要将两个条件链接在一起,AVX512比较屏蔽指令可以屏蔽写入屏蔽的操作,具有零屏蔽或合并屏蔽.
vcmpltps k1, zmm1, zmm2 ; k1 = zmm1<zmm2
vcmpltps k2{k1}{z}, zmm3, zmm4 ; k2 = (zmm3<zmm4) & (zmm1<zmm2)
k2zmm3k1为零时到处都是0,因为我们用作k1零掩码.
if (c(i) > 0) then
a(i) = b(i) ** 2
else
a(i) = b(i) ** 3
end if
这里常见的子表达式是b(i) * b(i).我们可以b(i)**3通过乘以b(i)一个额外的时间来得到它.
vmovups zmm1, [rsi] ; load a vector from b(i)
vmulps zmm2, zmm1, zmm1 ; zmm2 = zmm1*zmm1 = b(i)**2
AVX-512可以基于掩码进行合并,作为(几乎)任何其他指令的一部分.
vmulps zmm2{k1}, zmm2, zmm1 ; zmm2 *= zmm1 for elements where k1 is true
vmovups [rdi], zmm2 ; store all 16 elements into a(i)
BTW,AVX512为商店提供合并屏蔽功能.以前的SIMD指令集将加载[rdi],混合,然后存储回来[rdi].这意味着您可以a(i)比使用AVX1/AVX2更有效地实现每个元素条件的第二个循环(有时保持未修改).
把这一切放在一起:( NASM语法)
; x86-64 System V calling convention
; args: rdi = a() output array.
; rsi = b() input array
; rdx = c() array to be tested for positive numbers
; rcx = count (in elements)
; preferably all 64-byte aligned, but will work slowly if some aren't
; rcx must be >= 16, and a multiple of 16, because I didn't write any cleanup code
global square_or_cube
square_or_cube:
vxorps xmm0, xmm0,xmm0
.loop: ; do {
vcmpps k1, zmm0, [rdx], 21 ; _CMP_NLT_UQ ; !(0 < c(i))
vmovups zmm1, [rsi] ; load a vector from b(i)
vmulps zmm2, zmm1, zmm1 ; zmm2 = zmm1*zmm1 = b(i)**2
vmulps zmm2{k1}, zmm2, zmm1 ; zmm2 *= zmm1 for elements where k1 is true, otherwise unmodified.
vmovups [rdi], zmm2 ; store all 16 elements into a(i)
; TODO: unroll some and/or use indexed addressing mode tricks to save instructions
add rdi, 64 ; pointer increments
add rsi, 64
add rdx, 64
sub rcx, 16 ; count -= 16
ja .loop ; } while(count>0);
我用IACA分析了这一点(省略指针增量指令以模拟展开和更聪明的asm技巧).根据IACA的说法,即使是合并掩码vmulps也只是一个uop,而内存源指令微型融合到前端的单个uop.(商店也是如此.)这是我所希望的,IACA的输出在这种情况下看起来是正确的,尽管我无法访问SKL-SP硬件上的性能计数器来检查.
$ iaca.sh -arch SKX avx512-conditional
Intel(R) Architecture Code Analyzer Version - 2.3 build:246dfea (Thu, 6 Jul 2017 13:38:05 +0300)
Analyzed File - avx512-conditional
Binary Format - 64Bit
Architecture - SKX
Analysis Type - Throughput
Throughput Analysis Report
--------------------------
Block Throughput: 1.50 Cycles Throughput Bottleneck: FrontEnd
Port Binding In Cycles Per Iteration:
---------------------------------------------------------------------------------------
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 |
---------------------------------------------------------------------------------------
| Cycles | 1.5 0.0 | 0.0 | 1.0 1.0 | 1.0 1.0 | 1.0 | 1.5 | 1.0 | 1.0 |
---------------------------------------------------------------------------------------
N - port number or number of cycles resource conflict caused delay, DV - Divider pipe (on port 0)
D - Data fetch pipe (on ports 2 and 3), CP - on a critical path
F - Macro Fusion with the previous instruction occurred
* - instruction micro-ops not bound to a port
^ - Micro Fusion happened
# - ESP Tracking sync uop was issued
@ - SSE instruction followed an AVX256/AVX512 instruction, dozens of cycles penalty is expected
X - instruction not supported, was not accounted in Analysis
| Num Of | Ports pressure in cycles | |
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | |
---------------------------------------------------------------------------------
| 2^ | | | 1.0 1.0 | | | 1.0 | | | CP | vcmpps k1, zmm0, zmmword ptr [rdx], 0x15
| 1 | | | | 1.0 1.0 | | | | | | vmovups zmm1, zmmword ptr [rsi]
| 1 | 1.0 | | | | | | | | CP | vmulps zmm2, zmm1, zmm1
| 1 | 0.5 | | | | | 0.5 | | | CP | vmulps zmm2{k1}, zmm2, zmm1
| 2^ | | | | | 1.0 | | | 1.0 | | vmovups zmmword ptr [rdi], zmm2
| 1 | | | | | | | 1.0 | | | sub rcx, 0x10
| 0F | | | | | | | | | | jnbe 0xffffffffffffffdd
Total Num Of Uops: 8
AVX-512实际上有vfpclassps(C/C++内在的[_mm512_fpclass_ps_mask)4,asm文档与相关的表vfpclasspd(打包双))根据你选择的谓词对FP值进行分类.它可能比使用恰好为零的另一个寄存器的完全比较稍微有效.
(实际上,根据IACA,事实并非如此.两者都被列为3周期的延迟InstLatx64电子表格.昂纳雾的测量AVX2 cmpps上SKYLAKE微架构-S(非AVX512桌面芯片)示出了4个周期,所以它的奇怪的是,AVX512在生成掩码寄存器结果而不是向量时,版本的延迟较低.
我希望结果只对正数有效,我认为vfpclassps可以通过设置几乎所有的谓词位来获得-Inf,有限负,安静和信令NaN,-0.0和+0.0.
vfpclassps k1, [rdx], 0x1 | 0x2 | 0x4 | 0x10 | 0x40 | 0x80 ; QNaN | -0.0 | +0.0 | -Infinity | Negative (finite) | SNaN
; k1 = a 16-bit bitmap of which elements (from memory at [rdx]) need an extra multiply
vpfclassps有趣的是,因为它可以让你+0.0 -0.0和区分,比如你可以通过检查二进制表示的符号位(比如,你可以用AVX2 vblendps使用符号位作为混合控制,而无需首先做一个比较).
此外,在这种情况下,它在循环外部保存一条指令,设置全零的寄存器.
相关:AVX512具有乘以2**floor(x)(vscalefpd)的指令,但不是将数字乘以任意幂(整数或其他). Xeon Phi拥有AVX512ER,可以为您提供快速近似2**x(无地板x),但我们也不能直接使用指数函数,而SKL-SP无论如何都没有AVX512ER.
IACA_start/end的NASM宏:
我是根据iaca_marks.hC/C++标题编写的.
%if 1
%macro IACA_start 0
mov ebx, 111
db 0x64, 0x67, 0x90
%endmacro
%macro IACA_end 0
mov ebx, 222
db 0x64, 0x67, 0x90
%endmacro
%else
%define IACA_start
%define IACA_end
%endif
将它们包裹在您要分析的任何代码周围.
编译器可以在循环内分支.IDK,如果有的话会像这样编写代码,但他们当然可以.
; rdi = destination
; rsi = source
; edx = condition
; rcx = element count
global square_or_cube
square_or_cube:
.loop: ; do {
vmovups zmm1, [rsi] ; load a vector from b(i)
vmulps zmm2, zmm1, zmm1 ; zmm2 = zmm1*zmm1 = b(i)**2
test edx,edx
jz .only_square ; test-and-branch to conditionally skip the 2nd multiply
vmulps zmm2, zmm2, zmm1 ; zmm2 *= zmm1
.only_square:
vmovups [rdi], zmm2 ; store all 16 elements into a(i)
add rdi, 64 ; pointer increments
add rsi, 64
sub rcx, 16 ; count -= 16
ja .loop ; } while(count>0);
注意：这个答案主要讨论向量化时一个非常具体的内存访问问题，它主要在概念层面上应用于将对数组的一系列标量访问转换为向量化访问，而不假设底层数组的哪些部分被映射。在像 Fortran 这样的语言中，语言本身的语义可以保证数组是连续映射的，或者在进入循环之前进行边界检查可能足以避免下面提到的问题。
一般来说，这个答案不应该被视为对矢量化的良好处理，尤其是在 Fortran 中。另一个答案中出现了对矢量化问题的更全面的处理，该答案也专门解决了 AVX-512。
向量化条件经常被忽视的一个问题是，编译器可以通过混合或其他逐元素预测技术对您感兴趣的类型的条件循环进行向量化，前提是它们可以证明向量化访问的元素与在向量化中访问的元素相同。标量逐元素实现。如果指令集不提供按元素方式执行向量加载（尊重此条件），或者编译器无法使用它们，则这可能会有效地阻止向量化。
换句话说，如果循环体的所有路径访问相同的元素，编译器通常只能使用普通向量加载进行完全向量化。
根本原因是编译后的代码不得访问原始代码语义未访问的元素，即使它们后来被“混合掉”，因为这样做可能会导致错误！如果指令集不提供有条件地访问内存中的元素并抑制未选择的元素的错误的指令，那么这对于优化来说是一个重大障碍。
在您给出的示例中，这意味着 (1) 和 (3) 可以“不提升条件”进行矢量化，而 (2) 则不能，因为 (2) 访问a[i]且b[i]仅在if主体中，但如果if不执行则不会。当然，真正的编译器只会将一个简单的标志检查提升到循环之外，并且在这种myflag == false情况下根本不执行循环，所以这并不是一个很好的例子。
让我们看几个包含您所有示例的案例。首先，我们需要一个无法悬挂的标志 - 让我们只使用一个bool值数组。a因此，带有一个输出数组、两个输入数组b和c一个标志数组的有趣的通用循环f可能如下所示：
do i = 1, nx
if (f(i) > 0) then
a(i) = g(b(i), c(i));
else
a(i) = h(b(i), c(i));
end if
end do
根据f(i)与每个元素对应的标志，我们将函数g或应用于h输入元素b(i)和c(i)。根据我上面的条件，只有当和实际访问和的相同元素时g，我们才能进行矢量化。hbc
让我们继续看上面的两个实际工作示例：
void example1(bool* f, int* __restrict__ a, int* __restrict__ b, int* __restrict__ c, size_t n) {
for (size_t i = 0; i < n; i++) {
if (f[i]) {
a[i] = b[i];
} else {
a[i] = c[i];
}
}
}
void example2(bool* f, int* __restrict__ a, int* __restrict__ b, int* __restrict__ c, size_t n) {
for (size_t i = 0; i < n; i++) {
if (f[i]) {
a[i] = b[i] + c[i] ;
} else {
a[i] = b[i] - c[i] * 2 + 1 ;
}
}
}
两者具有相同的基本形式，但哪个更难矢量化？第一个是根据标志对b[i]或进行简单的直接分配。第二个是两者c[i]的更复杂的函数，并且在两条路径上显着不同。 b[i]c[i]
嗯，第二个更容易矢量化，因为它b[i]无条件访问c[i]。事实上，gcc由于某种原因，无法对其中任何一个进行矢量化。clang只对第二个进行向量化。有点令人惊讶的是，它成功地对两者icc进行了矢量化——因为它足够聪明，可以使用屏蔽加载来抑制卸载元素的错误。vpmaskmovd
我最初开始这个答案的想法是，访问不同的数组元素目前是当前编译器矢量化不可逾越的障碍，但那是因为我通常不检查icc. icc以这种方式使用蒙面动作对我来说实际上是新闻。所以障碍是存在的，但至少一些编译器可以克服它2。
作为开发人员，您通常知道这两个数组都是完全可访问的，因此可以安全地访问该范围内的所有元素，b并且c最好[0, n)将其传达给编译器。我尝试添加无条件虚拟语句，例如b[i] = b[i]; c[i] = c[i];or ... + c[i] * 0，它应该编译为空，但至少允许编译器看到语义上所有元素都被访问。确实“编译了”，但代码生成没有得到改进：不会发生额外的矢量化。可能在矢量化分析完成之前，它们已经在编译过程的早期被消除了，因此矢量化器会丢失信息。
除了非免费且不完全通用的屏蔽移动指令之外，还有其他方法可以改善这种情况吗？编译器可以利用其对平台内存保护模型的了解。例如，一旦访问了 x86 上 4K 页中的任何字节，就可以自由读取该页上的所有其他字节。人们可以想象一种复杂的实现，它以安全的标量代码开始，但一旦“注意到”对两个数组的写入，就会切换到页面其余部分的矢量化循环。
如果数组访问对齐，则可以使用类似的技巧：矢量化循环可以检查标志数组是否一致为 0 或一致 1，如果不是，则可以安全地使用简单的无条件无屏蔽读取实现，否则它将回退到更多认真执行。显然，只有当掩模很少统一或几乎总是统一3 时，这种转变才会有利可图，因此可能不太可能在实践中实施。
2至少如果 AVX 可用：icc如果将第一个示例限制为 AVX 之前的指令，则仍然无法矢量化第一个示例，因为 和 就是在那时vpmaskmovd/q引入vmaskmovps/pd的。
3因为在这种情况下，如果您已经确定掩模是均匀的，则可以通过仅执行选定的一侧而无条件地实现操作，而if无需根据它是统一的0还是统一的进行任何掩蔽/混合1。因此，您最终会得到内部实现的三个循环：全零标志情况、全一标志情况和混合标志情况，当下一个标志向量与当前循环不同时，它们之间会跳转。