gfortran或ifort编译器在计算两个数组的乘积之和时，是否会明智地使用SIMD指令？

Question

我有一些使用numpy编写的代码，并考虑将其迁移到Fortran以获得更好的性能。

我多次执行的一个操作是对两个数组的逐元素乘积进行求和：

sum(A*B)

很明显，融合乘加指令将有助于此。我的当前处理器不支持这些指令，所以我无法测试。但是，我可能会升级到一款支持FMA3的新处理器（Intel Haswell处理器）。

请问是否可以使用“-march=native”（或ifort等效项）编译程序即可让编译器（gfortran或ifort）智能地使用SIMD指令来优化该代码，还是您认为我需要小心地处理编译器或代码？

- lnmaurer

---

在这种情况下，不使用 DOT_PRODUCT(a, b) 有什么理由吗？ - Marat Dukhan

根据文档，DOT_PRODUCT 需要一维数组，但我正在使用二维数组。 - lnmaurer

大多数numpy已经用C编写 - 它应该是相当快的... - naught101

Answer 1

如果你在使用具有SIMD的机器上使用-march=native，编译器应该会生成SIMD指令，尽管我总是使用ifort的-xHost标志。但是我不太确定如何让它们“明智地”执行此操作。我感觉，在-O3级别上，ifort和gfortran都倾向于过度聚合矢量化（即，他们比他们应该使用更多的SIMD功能）。很多时候，为了获得最有效的代码，我必须关闭向量化。当然，这可能对您来说是真实的，也可能不是。

通常最好使用优化为此任务的向量库。您可以使用MKL中的vdmul或GSL中的gsl_vector_mul来完成此操作。

使用-march=NEWARCH将导致针对架构NEWARCH进行调整的代码，但无法运行在早期的架构上。您可以使用-mtune=NEWARCH标志，其中NEWARCH是您新处理器的架构。这将生成为新架构调整的代码，但仍可在旧架构上执行。由于您还没有新设备，因此目前可能需要使用-mtune。

使用ifort，您可以使用矢量化报告标志显示已矢量化程序的哪个部分。例如，ifort标志-vec-report=1将在编译过程中提供此类信息。我确信gfortran中会有等效的标志。

- Xiaolei Zhu

---

-march=NEWARCH 的提示非常有帮助。 - lnmaurer

Answer 2

sum(a*b)比dot_product(a,b)更好的向量化已经过时。您展示的代码正在使用串行AVX2 fma指令。

在没有间接索引或其他复杂性（简单循环本身）的dot_product实现中，fma可能会比simd并行乘法和加法指令的组合慢，因为乘法可以在延迟关键路径之外完成。gfortran对于dot_product的并行simd fma在更复杂的情况下可能非常有效。

您需要使用-O2 -ftree-vectorize -ffast-math -march=native或-O3 -ffast-math -march=native（以及适当的向量长度）来进行矢量化，而且gfortran可能无法在OpenMP并行区域内进行矢量化。

gfortran 4.9似乎已经取消了选项-ftree-vectorizer-verbose。 -fdump-tree-vect将矢量化传递的详细信息写入.vect文件，并为不同的主要gcc版本选择不同的名称。

- tim18

---

dot_product函数是否适用于多维数组？我看到的文档似乎只适用于一维数组。 - lnmaurer

dot_product 可以处理任意秩（rank）数组的1D部分。matmul是2D数组的相应约简运算。 - tim18

Answer 3

感谢朱晓磊的建议，我现在知道gfortran将使用融合乘加来优化sum(A*B)。例如，对于以下代码：

program test implicit none

real, dimension(7) :: a, b

a = (/ 2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0 /)

b = (/ 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0 /)

print *, sum(a*b)
endprogram

我可以使用f95 sum.f95 -o sum -O3 -march=core-avx2进行编译，objdump -d sum | grep vfmadd会显示：

40088b: c4 e2 71 99 44 24 30 vfmadd132ss 0x30(%rsp),%xmm1,%xmm0

400892: c4 e2 69 b9 44 24 34 vfmadd231ss 0x34(%rsp),%xmm2,%xmm0

400899: c4 e2 61 b9 44 24 38 vfmadd231ss 0x38(%rsp),%xmm3,%xmm0

4008a0: c4 e2 59 b9 44 24 3c vfmadd231ss 0x3c(%rsp),%xmm4,%xmm0

4008a7: c4 e2 51 b9 44 24 40 vfmadd231ss 0x40(%rsp),%xmm5,%xmm0

4008ae: c4 e2 49 b9 44 24 44 vfmadd231ss 0x44(%rsp),%xmm6,%xmm0

4008b5: c4 e2 41 b9 44 24 48 vfmadd231ss 0x48(%rsp),%xmm7,%xmm0

这是一段编程代码，使用了vfmadd指令进行浮点数运算。每行代码都对应着一个特定的内存地址和操作寄存器。

所以gfortran展开了循环并放入7个融合乘加指令。如果我创建更大、更随机的多维数组，仍然会看到vfmadd231ss弹出一次（因此它不会展开循环）。

- lnmaurer