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

Question

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

我应该使用一个二维数组 A(i,j) 和一个一维数组 B(i)
我要循环2次，第一次循环应该是50000次，第二次是5次（不能改）。
上面的第2点应该循环10000次。

我用两个版本编写代码（我称它们为 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 基于 "column major"，因此第二种情况会比第一种情况更快，因为我（在第二种情况下）执行了数组最内侧的循环（在这种情况下，i=1, ..., 5）。

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

第一种情况：经过时间 = 29.187 秒
第二种情况：经过时间 = 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 访问差。

Answer 2

部分优化基于实施选择的启发式方法。您不会得到完全确定性的解释。当用允许完整数组操作的语言编写程序时，例如 Fortran，要点是：如果您可以将某些内容表示为完整的数组/或向量操作，请执行此操作并让编译器为您生成循环。那样的话，唯一要做的就是用最有效的编译器进行测试并坚持下去。

请记住，它依赖于编译器，因此会随编译器而变化。例如，如果我按原样使用您的程序（更改 k 迭代次数从 100000 到 10000 以加速它，并平均每次运行 5 次），这是我电脑上不同编译器的时间。我没有时间为您检查新版本的编译器。

pgf90 14
A: 5.05 s
B: 0.18 s

gfortran 4.9
A: 1.05 s
B: 4.02 s

ifort 13
A: 2.01 s
B: 1.08 s

你可以看到 gfortran 告诉你 B 不好的地方，pgfortran 完全相反，完全吹嘘 A

的结果

现在，如果我矢量化让编译器完成这项工作。这里我只修改A 并消除 I 循环得到这个：

DO K=1, 10000
    DO J=1, 5
        A(I,J)=A(I,J)+SQRT(B(I))
    END DO
END DO

然后我们得到这个（只有程序A）

pgf90 14: 5.05 s

gfortran 4.9: 5.04

ifort 13: 2.02 s

pgfortran 和 ifort 是稳定的，而 gfortran 在第一种情况下利用了一个技巧，可能是 haraldkl 的建议（见因素 5）。当我们向量化时，技巧并不明显 gfortran 表现不佳。似乎ifort和pgfortran 只需重写循环即可获得正确的顺序。

如果我变得更聪明并且也能限制 K 环

    DO J=1, 5
        A(:,J)=A(:,J)+10000*SQRT(B(:)) ! this seems to be final result
    END DO

然后我们得到这个（只有程序A）

pgf90 14: 0.001

gfortran 4.9: 0.001

ifort 13: 0.001 s

所有编译器都变得等效，因为几乎没有什么可以优化的。您会看到，只需使用数组操作即可优化所有内容。

更新 High Performance Mark 在评论中指出，如果编译器发现结果未被使用，它实际上可能会跳过所有计算，这在某些实现中可能会发生。这个答案中给出的结果说明了这种可能性，尽管我一开始并没有提到它。为了防止编译器完全跳过代码，我在计算后（计时后）打印了结果所有元素的总和；结果与小数点后的第 3 位数字相同，足以得到 ~372684326034.9146 (~10^12) 的结果。这足以确保编译器进行实际计算。我完全忘了在答案中提到它。

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

How to choose the best configuration of 2D array A(i,j)

arrays

fortran

gfortran

fortran95