arrays - 如何选择二维数组 A(i,j) 的最佳配置

Question

希望你能解释一下这件事。我正在使用 Fortran 并实际编写关于 CFD 主题的代码，下面(只是为了简单起见，只是为了举例)是对我的案例的简短解释:

1. 我应该使用一个二维数组 A(i,j) 和一个一维数组 B(i)
2. 我必须做 2 次循环，第一次循环应该是 50,000 次，第二次是 5 次(不能更改)。
3. 上面的第 2 个点应该循环 10,000 次。

我用 2 个版本编写代码(我称它们为 Prog_A 和 Prog_B)。

第一个如下图:

PROGRAM PROG_A
 REAL*8, DIMENSION(50000,5):: A
 REAL*8, DIMENSION(50000)::B  
 REAL*8:: TIME1,TIME2
       !Just for initial value
           DO I=1, 50000
              A(I,1)=1.0
              A(I,2)=2.0
              A(I,3)=3.0
              A(I,4)=4.0
              A(I,5)=5.0

              B(I)=I 
           END DO
       !Computing computer's running time (start)
           CALL CPU_TIME(TIME1)
                 DO K=1, 100000
                    DO I=1, 50000             !Array should be computed first for 50,000 elements (can't be changed)
                        DO J=1, 5
                           A(I,J)=A(I,J)+SQRT(B(I))   
                        END DO
                    END DO
                 END DO
       !Computing computer's running time (finish)
           CALL CPU_TIME(TIME2)
             PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PROG_A

第二个是:

PROGRAM PROG_B
 REAL*8, DIMENSION(5,50000):: A
 REAL*8, DIMENSION(50000)::B  
 REAL*8:: TIME1,TIME2
       !Just for initial value
           DO J=1, 50000
              A(1,J)=1.0
              A(2,J)=2.0
              A(3,J)=3.0
              A(4,J)=4.0
              A(5,J)=5.0

              B(J)=J 
           END DO
       !Computing computer's running time (start)
           CALL CPU_TIME(TIME1)
                 DO K=1, 100000
                    DO J=1, 50000             !Array should be computed first for 50,000 elements (can't be changed)
                        DO I=1, 5
                           A(I,J)=A(I,J)+SQRT(B(J))   
                        END DO
                    END DO
                 END DO
       !Computing computer's running time (finish)
           CALL CPU_TIME(TIME2)
             PRINT *, 'Elapsed real time = ', TIME2-TIME1, 'second(s)'
END PROGRAM PROG_B

如您所见，不同之处在于第一个我使用了二维数组 A(50000,5)，第二个我使用了二维数组 A(5,50000)。

据我所知，由于 Fortran 是基于“主列”的，因此第二种情况会比第一种情况更快，因为我(在第二种情况中)执行了数组最内侧的循环(在这种情况下案例 i=1, ..., 5)。

但在 gfortran 上编译后(使用 -O3 优化)，我发现第二个甚至比第一个慢得多。这是结果:

1. 第一种情况:耗时 = 29.187 秒
2. 第二种情况:耗时 = 70.496 秒

谁能解释一下为什么？

PS:两种情况的结果肯定是一样的。

Answer 1

编译器可能会做这样的事情:

    DO K=1, 100000
       DO I=1, 50000
          tmp = sqrt(b(i))
          A(I,1) = A(I,1) + tmp
          A(I,2) = A(I,2) + tmp
          A(I,3) = A(I,3) + tmp
          A(I,4) = A(I,4) + tmp
          A(I,5) = A(I,5) + tmp
       END DO
    END DO

在 Prog_A 中，这为您提供了步幅为 1 的良好访问模式。如果您像在 Prog_B 中那样更改索引的顺序，您将获得步幅为5 对于此代码。其效果取决于机器，但肯定比简单的 stride-1 访问更差。