循环优化

4

我正在尝试理解在源代码中可以进行哪些缓存或其他优化来使这个循环更快。我认为它相当友好,但是是否有任何专家能够挤出更多性能来调整这段代码呢?

DO K = 1, NZ 
    DO J = 1, NY
        DO I = 1, NX

             SIDEBACK  = STEN(I-1,J-1,K-1) + STEN(I-1,J,K-1) + STEN(I-1,J+1,K-1) + &
                         STEN(I  ,J-1,K-1) + STEN(I  ,J,K-1) + STEN(I  ,J+1,K-1) + &
                         STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1) 

             SIDEOWN   = STEN(I-1,J-1,K)   + STEN(I-1,J,K)   + STEN(I-1,J+1,K) + &
                         STEN(I  ,J-1,K)   + STEN(I  ,J,K)   + STEN(I  ,J+1,K) + &
                         STEN(I+1,J-1,K)   + STEN(I+1,J,K)   + STEN(I+1,J+1,K)

             SIDEFRONT = STEN(I-1,J-1,K+1) + STEN(I-1,J,K+1) + STEN(I-1,J+1,K+1) + &
                         STEN(I  ,J-1,K+1) + STEN(I  ,J,K+1) + STEN(I  ,J+1,K+1) + &
                         STEN(I+1,J-1,K+1) + STEN(I+1,J,K+1) + STEN(I+1,J+1,K+1)

             RES(I,J,K) = ( SIDEBACK + SIDEOWN + SIDEFRONT ) / 27.0

        END DO
    END DO
END DO

2
只是为了确认一下,这是Fortran语言,而STEN是一个列主序数组? - szmate1618
@szmate1618 正确 - Manolete
1
我认为你可以做得更好。当前的代码会多次计算每个部分的总和。我的意思是,一个迭代的 STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1) 将成为下一个迭代的 STEN(I,J-1,K-1) + STEN(I,J,K-1) + STEN(I,J+1,K-1),因此无需再次获取和计算,可以存储这些部分结果。我已经在进行完整实现,很快就会发布。 - szmate1618
1
等等,我有点蠢。如果Fortran编译器展开循环,它可能能够完成我正在尝试做的同样的优化。如果STEN被传递为输入参数,它可以非常容易地实现。没有别名,极端循环展开是Fortran的招牌。 - szmate1618
1
这些都不是指针,对吧? - Ross
显示剩余5条评论
1个回答

5

好的,我认为我已经尽力了,但不幸的是我的结论是,除非你愿意进入并行化,否则优化的空间不太大。让我们看看为什么,看看你能做什么和不能做什么。

编译器优化

现今的编译器在优化代码方面非常出色,比人类要好得多。依靠编译器进行优化的另一个好处是它们不会破坏您源代码的可读性。无论您做什么(在优化速度时),始终尝试使用每种合理的编译器标志组合。您甚至可以尝试多个编译器。个人只使用了gfortran(包含在GCC中)(操作系统为64位Windows),我相信它具有高效且正确的优化技术。

-O2 几乎总能显著提高速度,但即使是 -O3 也是一个安全的选择(其中包括美味的循环展开)。对于这个问题,我还尝试了 -ffast-math-fexpensive-optimizations,它们没有任何可衡量的影响,但 -march-corei7(针对Core i7的特定CPU架构调整)有影响,所以我用 -O3 -march-corei7 进行了测试。

到底有多快?

我编写了以下代码来测试您的解决方案,并使用 -O3 -march-corei7 进行编译。通常它会在0.78-0.82秒内运行。

program benchmark
implicit none

real :: start, finish
integer :: I, J, K
real :: SIDEBACK, SIDEOWN, SIDEFRONT

integer, parameter :: NX = 600
integer, parameter :: NY = 600
integer, parameter :: NZ = 600
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: STEN
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: RES
call random_number(STEN)

call cpu_time(start)
DO K = 1, NZ
    DO J = 1, NY
        DO I = 1, NX

            SIDEBACK =  STEN(I-1,J-1,K-1) + STEN(I-1,J,K-1) + STEN(I-1,J+1,K-1) + &
                        STEN(I  ,J-1,K-1) + STEN(I  ,J,K-1) + STEN(I  ,J+1,K-1) + &
                        STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1)

            SIDEOWN =   STEN(I-1,J-1,K)   + STEN(I-1,J,K)   + STEN(I-1,J+1,K) + &
                        STEN(I  ,J-1,K)   + STEN(I  ,J,K)   + STEN(I  ,J+1,K) + &
                        STEN(I+1,J-1,K)   + STEN(I+1,J,K)   + STEN(I+1,J+1,K)

            SIDEFRONT = STEN(I-1,J-1,K+1) + STEN(I-1,J,K+1) + STEN(I-1,J+1,K+1) + &
                        STEN(I  ,J-1,K+1) + STEN(I  ,J,K+1) + STEN(I  ,J+1,K+1) + &
                        STEN(I+1,J-1,K+1) + STEN(I+1,J,K+1) + STEN(I+1,J+1,K+1)

            RES(I,J,K) = ( SIDEBACK + SIDEOWN + SIDEFRONT ) / 27.0

        END DO
    END DO
END DO
call cpu_time(finish)
!Use the calculated value, so the compiler doesn't optimize away everything.
!Print the original value as well, because one can never be too paranoid.
print *, STEN(1,1,1), RES(1,1,1)
print '(f6.3," seconds.")',finish-start

end program

好的,编译器只能带我们走到这里。接下来是什么?

存储中间结果?

正如你从问号中可以猜测到的那样,这个方法并没有真正起作用。抱歉。但是让我们不要急于求成。 如评论所述,你当前的代码会多次计算每个部分的和,这意味着一个迭代的STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1)将成为下一个迭代的STEN(I,J-1,K-1) + STEN(I,J,K-1) + STEN(I,J+1,K-1),因此无需再次获取和计算,可以存储这些部分结果。 问题在于,我们不能存储太多的部分结果。正如你所说,你的代码已经相当友好,每个部分和你存储意味着你可以在L1缓存中存储少量数组元素。我们可以存储一些值,从I的最后几次迭代中(索引I-2I-3等)选择一些值,但编译器几乎肯定已经这样做了。我有两个证据支持这种怀疑。首先,我的手动循环展开使程序变慢了,大约5%。

DO K = 1, NZ
    DO J = 1, NY
        DO I = 1, NX, 8

            SIDEBACK(0) = STEN(I-1,J-1,K-1) + STEN(I-1,J,K-1) + STEN(I-1,J+1,K-1)
            SIDEBACK(1) = STEN(I  ,J-1,K-1) + STEN(I  ,J,K-1) + STEN(I  ,J+1,K-1)
            SIDEBACK(2) = STEN(I+1,J-1,K-1) + STEN(I+1,J,K-1) + STEN(I+1,J+1,K-1)
            SIDEBACK(3) = STEN(I+2,J-1,K-1) + STEN(I+2,J,K-1) + STEN(I+2,J+1,K-1)
            SIDEBACK(4) = STEN(I+3,J-1,K-1) + STEN(I+3,J,K-1) + STEN(I+3,J+1,K-1)
            SIDEBACK(5) = STEN(I+4,J-1,K-1) + STEN(I+4,J,K-1) + STEN(I+4,J+1,K-1)
            SIDEBACK(6) = STEN(I+5,J-1,K-1) + STEN(I+5,J,K-1) + STEN(I+5,J+1,K-1)
            SIDEBACK(7) = STEN(I+6,J-1,K-1) + STEN(I+6,J,K-1) + STEN(I+6,J+1,K-1)
            SIDEBACK(8) = STEN(I+7,J-1,K-1) + STEN(I+7,J,K-1) + STEN(I+7,J+1,K-1)
            SIDEBACK(9) = STEN(I+8,J-1,K-1) + STEN(I+8,J,K-1) + STEN(I+8,J+1,K-1)

            SIDEOWN(0) = STEN(I-1,J-1,K) + STEN(I-1,J,K) + STEN(I-1,J+1,K)
            SIDEOWN(1) = STEN(I  ,J-1,K) + STEN(I  ,J,K) + STEN(I  ,J+1,K)
            SIDEOWN(2) = STEN(I+1,J-1,K) + STEN(I+1,J,K) + STEN(I+1,J+1,K)
            SIDEOWN(3) = STEN(I+2,J-1,K) + STEN(I+2,J,K) + STEN(I+2,J+1,K)
            SIDEOWN(4) = STEN(I+3,J-1,K) + STEN(I+3,J,K) + STEN(I+3,J+1,K)
            SIDEOWN(5) = STEN(I+4,J-1,K) + STEN(I+4,J,K) + STEN(I+4,J+1,K)
            SIDEOWN(6) = STEN(I+5,J-1,K) + STEN(I+5,J,K) + STEN(I+5,J+1,K)
            SIDEOWN(7) = STEN(I+6,J-1,K) + STEN(I+6,J,K) + STEN(I+6,J+1,K)
            SIDEOWN(8) = STEN(I+7,J-1,K) + STEN(I+7,J,K) + STEN(I+7,J+1,K)
            SIDEOWN(9) = STEN(I+8,J-1,K) + STEN(I+8,J,K) + STEN(I+8,J+1,K)

            SIDEFRONT(0) = STEN(I-1,J-1,K+1) + STEN(I-1,J,K+1) + STEN(I-1,J+1,K+1)
            SIDEFRONT(1) = STEN(I  ,J-1,K+1) + STEN(I  ,J,K+1) + STEN(I  ,J+1,K+1)
            SIDEFRONT(2) = STEN(I+1,J-1,K+1) + STEN(I+1,J,K+1) + STEN(I+1,J+1,K+1)
            SIDEFRONT(3) = STEN(I+2,J-1,K+1) + STEN(I+2,J,K+1) + STEN(I+2,J+1,K+1)
            SIDEFRONT(4) = STEN(I+3,J-1,K+1) + STEN(I+3,J,K+1) + STEN(I+3,J+1,K+1)
            SIDEFRONT(5) = STEN(I+4,J-1,K+1) + STEN(I+4,J,K+1) + STEN(I+4,J+1,K+1)
            SIDEFRONT(6) = STEN(I+5,J-1,K+1) + STEN(I+5,J,K+1) + STEN(I+5,J+1,K+1)
            SIDEFRONT(7) = STEN(I+6,J-1,K+1) + STEN(I+6,J,K+1) + STEN(I+6,J+1,K+1)
            SIDEFRONT(8) = STEN(I+7,J-1,K+1) + STEN(I+7,J,K+1) + STEN(I+7,J+1,K+1)
            SIDEFRONT(9) = STEN(I+8,J-1,K+1) + STEN(I+8,J,K+1) + STEN(I+8,J+1,K+1)

            RES(I    ,J,K) = (  SIDEBACK(0) + SIDEOWN(0) + SIDEFRONT(0) + &
                                SIDEBACK(1) + SIDEOWN(1) + SIDEFRONT(1) + &
                                SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) ) / 27.0
            RES(I + 1,J,K) = (  SIDEBACK(1) + SIDEOWN(1) + SIDEFRONT(1) + &
                                SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) + &
                                SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) ) / 27.0
            RES(I + 2,J,K) = (  SIDEBACK(2) + SIDEOWN(2) + SIDEFRONT(2) + &
                                SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) + &
                                SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) ) / 27.0
            RES(I + 3,J,K) = (  SIDEBACK(3) + SIDEOWN(3) + SIDEFRONT(3) + &
                                SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) + &
                                SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) ) / 27.0
            RES(I + 4,J,K) = (  SIDEBACK(4) + SIDEOWN(4) + SIDEFRONT(4) + &
                                SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) + &
                                SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6)   ) / 27.0
            RES(I + 5,J,K) = (  SIDEBACK(5) + SIDEOWN(5) + SIDEFRONT(5) + &
                                SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6) + &
                                SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) ) / 27.0
            RES(I + 6,J,K) = (  SIDEBACK(6) + SIDEOWN(6) + SIDEFRONT(6) + &
                                SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) + &
                                SIDEBACK(8) + SIDEOWN(8) + SIDEFRONT(8) ) / 27.0
            RES(I + 7,J,K) = ( SIDEBACK(7) + SIDEOWN(7) + SIDEFRONT(7) + &
                                SIDEBACK(8) + SIDEOWN(8) + SIDEFRONT(8) + &
                                SIDEBACK(9) + SIDEOWN(9) + SIDEFRONT(9) ) / 27.0

        END DO
    END DO
END DO

更糟糕的是,很容易证明我们已经非常接近理论最小执行时间。为了计算所有这些平均值,我们至少需要访问每个元素一次,并将它们除以27.0。因此,您永远无法比以下代码更快,在我的计算机上执行时间在0.48-0.5秒之间。

program benchmark
implicit none

real :: start, finish
integer :: I, J, K

integer, parameter :: NX = 600
integer, parameter :: NY = 600
integer, parameter :: NZ = 600
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: STEN
real, dimension (0 : NX + 2, 0 : NY + 2, 0 : NZ + 2) :: RES
call random_number(STEN)

call cpu_time(start)
DO K = 1, NZ
    DO J = 1, NY
        DO I = 1, NX

            !This of course does not do what you want to do,
            !this is just an example of a speed limit we can never surpass.
            RES(I, J, K) = STEN(I, J, K) / 27.0

        END DO
    END DO
END DO
call cpu_time(finish)
!Use the calculated value, so the compiler doesn't optimize away everything.
print *, STEN(1,1,1), RES(1,1,1)
print '(f6.3," seconds.")',finish-start

end program

但是,即使是负面结果也是一种结果。如果仅仅访问每个元素一次(并除以27.0)就占用了大部分执行时间,那就意味着内存访问是瓶颈。然后你可以优化

更少的数据

如果您不需要64位双精度浮点数的完全精度,可以使用real(kind=4)类型声明数组。但是,也许您的实数已经是4字节。在这种情况下,我相信一些Fortran实现支持非标准16位双精度浮点数,或者根据您的数据,您可以使用整数(可能是浮点数乘以一个数字,然后舍入为整数)。基本类型越小,您就可以将更多元素放入缓存中。当然,最理想的是integer(kind=1),与real(kind=4)比较,在我的机器上导致了超过2倍的加速。但这取决于您需要的精度。

更好的局部性

当您需要来自相邻列的数据时,列主数组很慢,而对于相邻行的数据,行主数组很慢。 幸运的是,有一种奇妙的存储数据的方式,称为Z-顺序曲线,它在计算机图形学中具有类似于您使用情况的应用。 我不能保证它会有所帮助,也许它会非常低效,但也可能不会。抱歉,老实说我没有感觉要自己实现它。

并行化

说到计算机图形学,这个问题在GPU上非常易于并行化,但是如果您不想走得太远,您可以使用普通的多核CPU。 Fortran Wiki似乎是搜索Fortran并行化库的好地方。


3
real(kind=4) 并不一定是 4 字节,有些编译器中甚至不存在它。类似地,real(kind=1) 在某些编译器中存在,但可能只是通常的默认实数类型。 - Vladimir F Героям слава
展开可能不是有效的,因为您没有使用括号或其他方法来确保编译器识别相邻赋值之间的公共子表达式。此外,您应该尝试在多个级别上进行少量展开。内部循环中的展开是在-funroll-loops和--param选项的控制下自动完成的,因此写入它几乎没有意义。 - tim18
自动展开可能是无效的,手动展开也可能是不必要的。但是自动展开无效并且手动展开不必要是不可能的,我认为这是矛盾的。其中一个一定比另一个快。 我认为我的手动展开和手动公共子表达式优化比自动展开慢,这表明编译器已经执行了这些优化,并且比我做得更好。 - szmate1618
编译器成功将内部循环向量化,这正是您所期望的。我尝试阻塞外部循环,但这会减慢运行时间。 - Manolete
我也同意,内存访问是导致该循环非常缓慢的原因,但是除非尽可能多地使用L1,否则我没有更好的加速内存访问的想法。 - Manolete
1
这个例子可能不需要循环阻塞,它可能可以一直通过内部循环进入下一个中间循环,而不会丢失任何预取数据。硬件事件分析工具旨在帮助加速此类实验。 - tim18

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接