Fortran 中的矢量化总和

Posted 2023-02-16

【中文标题】Fortran 中的矢量化总和

【英文标题】：vectorized sum in Fortran

【发布时间】：2015-08-27 18:48:32

【问题描述】：

我正在使用gfortran 和-mavx 编译我的Fortran 代码，并验证了一些指令是通过objdump 进行矢量化的，但是我没有得到我期望的速度改进，所以我想确保正在对以下参数进行矢量化（这条单条指令占运行时间的约 50%）。

我知道有些指令可以向量化，而有些则不能，所以我想确保可以：

sum(A(i1:i2,ir))

同样，由于我是在一个非常大的矩阵上执行此操作，因此这条单行占用了大约 50% 的运行时间。我可以提供更多关于为什么我这样做的信息，但只要说这是必要的就足够了，尽管我可以在必要时重组内存（例如，我可以将总和作为sum(A(ir,i1:i2))如果可以改为矢量化。

这条线是否被矢量化了？我怎么知道？如果没有被矢量化，如何强制矢量化？

编辑：感谢 cmets，我现在意识到我可以通过 -ftree-vectorizer-verbose 检查此求和的矢量化，并看到这是 not 矢量化。我已将代码重组如下：

tsum = 0.0d0
tn = i2 - i1 + 1
tvec(1:tn) = A(i1:i2, ir)
do ii = 1,tn
    tsum = tsum + tvec(ii)
enddo

当我打开-funsafe-math-optimizations 时，这ONLY 会矢量化，但由于矢量化，我确实看到速度又提高了 70%。问题仍然存在：为什么sum(A(i1:i2,ir)) 不矢量化，我怎样才能得到一个简单的sum 来矢量化？

【问题讨论】：

不是骗子，而是：***.com/questions/23644203/…

Fortran 以列优先方式访问数组，即您当前的版本比sum(A(ir,i1:i2)) 快得多。您当前的版本访问了一块连续的内存。

关于您的编辑：非常有趣。我本来希望总和被矢量化。您是否使用 sum vs loop 进行了实际基准测试？你玩过优化开关吗？更新：have you tried -ffast-math?

OpenMP 4 "simd" 子句和缩减？

您也可以尝试-fopt-info 获取更多信息（这在man gcc 中有记录）。我猜，有某种形式的循环展开。

【参考方案1】：

事实证明，除非我包含 -ffast-math 或 -funsafe-math-optimizations，否则我无法使用矢量化。

我玩过的两个代码sn-ps是：

tsum = 0.0d0
tvec(1:n) = A(i1:i2, ir)
do ii = 1,n
    tsum = tsum + tvec(ii)
enddo

和

tsum = sum(A(i1:i2,ir))

以下是我使用不同编译选项运行第一个代码 sn-p 时得到的时间：

10.62 sec ... None
10.35 sec ... -mtune=native -mavx
 7.44 sec ... -mtune-native -mavx -ffast-math
 7.49 sec ... -mtune-native -mavx -funsafe-math-optimizations

最后，通过这些相同的优化，我可以对tsum = sum(A(i1:i2,ir)) 进行矢量化以获得

 7.96 sec ... None
 8.41 sec ... -mtune=native -mavx
 5.06 sec ... -mtune=native -mavx -ffast-math
 4.97 sec ... -mtune=native -mavx -funsafe-math-optimizations

当我们将sum 和-mtune=native -mavx 与-mtune=native -mavx -funsafe-math-optimizations 进行比较时，它显示出约70% 的加速。 （请注意，这些只运行一次 - 在我们发布之前，我们将对多次运行进行真正的基准测试）。

不过，我确实受到了一点打击。当我使用 -f 选项时，我的值会略有变化。没有它们，我的变量（v1，v2）的错误是：

v1 ... 5.60663e-15     9.71445e-17     1.05471e-15
v2 ... 5.11674e-14     1.79301e-14     2.58127e-15

但是通过优化，错误是：

v1 ... 7.11931e-15     5.39846e-15     3.33067e-16
v2 ... 1.97273e-13     6.98608e-14     2.17742e-14

这表明确实发生了一些不同的事情。

【讨论】：

好吧，我猜unsafeness 的优化来源之一可能来自重新安排操作顺序。我不认为机器精度级别的错误是致命的。

【参考方案2】：

您的显式循环版本仍然会以与矢量化版本不同的顺序添加 FP。矢量版本使用 4 个累加器，每个累加器获取每 4 个数组元素。

您可以编写源代码以匹配矢量版本的功能：

tsum0 = 0.0d0
tsum1 = 0.0d0
tsum2 = 0.0d0
tsum3 = 0.0d0
tn = i2 - i1 + 1
tvec(1:tn) = A(i1:i2, ir)
do ii = 1,tn,4   ! count by 4
    tsum0 = tsum0 + tvec(ii)
    tsum1 = tsum1 + tvec(ii+1)
    tsum2 = tsum2 + tvec(ii+2)
    tsum3 = tsum3 + tvec(ii+3)
enddo

tsum = (tsum0 + tsum1) + (tsum2 + tsum3)

这可能在没有-ffast-math的情况下进行矢量化。

FP add 具有多周期延迟，但每个时钟吞吐量只有一两个，因此您需要 asm 使用多个向量累加器来使 FP add 单元饱和。 Skylake 每个时钟可以做两个 FP 添加，延迟 = 4。以前的 Intel CPU 每个时钟执行一个，延迟 = 3。所以在 Skylake 上，你需要 8 个向量累加器来使 FP 单元饱和。当然，它们必须是 256b 向量，因为 AVX 指令与 SSE 向量指令一样快，但工作量却是 SSE 向量指令的两倍。

用 8 * 8 累加器变量编写源代码会很荒谬，所以我猜你需要-ffast-math，或者一个告诉编译器不同的操作顺序都可以的 OpenMP pragma。

显式展开源意味着您必须处理不是向量宽度*展开倍数的循环计数。如果你对事物施加限制，它可以帮助编译器避免生成多个版本的循环或额外的循环设置/清理代码。

【讨论】：

如果 tn 不是 4 的倍数，那么循环是否会漏掉一些或走得太远？

@Laurbert515：是的，一个或另一个。如果我知道 Fortran，我可以告诉你哪个。 :P 这是我最后一段描述的问题：展开导致需要清理代码来处理奇数的迭代。通过使用 -ffast-math 或 OpenMP 让编译器担心这一点，将为您节省大量调试时间（或者至少在开发时间进行思考以确保您不会首先引入错误。）