MPI_REDUCE 返回 1000x1000x6 REAL 数组的错误答案

Question

我已经与 gfortran 大气传输模型代码中的 MPI_REDUCE 调用进行了几天的斗争，具有良好的输入参数，但在 master 的 recvbuf 中返回非常不合理的结果。我已经能够在一个简单的示例中复制该问题，如下所示：

PROGRAM TEST

    USE mpi

    IMPLICIT NONE

    INTEGER my_rank, size, ierror
    INTEGER, PARAMETER :: nx=1000, ny=1000, nz=6
    INTEGER :: buffsize

    REAL, DIMENSION(nx, ny, nz) :: u, v

    call MPI_INIT(ierror)
    call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierror)
    call MPI_COMM_RANK(MPI_COMM_WORLD, my_rank, ierror)

    PRINT *, 'my_rank, size: ', my_rank, size

    buffsize = nx*ny*nz

    u = my_rank + 1

    PRINT *, 'PE: ', my_rank, ', Before reduce, SHAPE(u): ', SHAPE(u)
    PRINT *, 'PE: ', my_rank, ', Before reduce, SUM(u): ', SUM(u)

    CALL MPI_REDUCE(u, v, buffsize, MPI_REAL, &
&                   MPI_SUM, 0, MPI_COMM_WORLD, ierror)

    CALL MPI_BARRIER(MPI_COMM_WORLD, ierror)

    PRINT *, 'PE: ', my_rank, ', after reduce, ierror: ', ierror
    PRINT *, 'PE: ', my_rank, ', after reduce, SUM(u): ', SUM(u)
    PRINT *, 'PE: ', my_rank, ', after reduce, SUM(v): ', SUM(v)

    CALL MPI_FINALIZE(ierror)

END PROGRAM test

返回：

mpirun -np 2 ./test3
 my_rank, size:            0           2
 my_rank, size:            1           2
 PE:            1 , Before reduce, SHAPE(u):         1000        1000           6
 PE:            0 , Before reduce, SHAPE(u):         1000        1000           6
 PE:            0 , Before reduce, SUM(u):    6000000.00    
 PE:            1 , Before reduce, SUM(u):    12000000.0    
 PE:            0 , after reduce, ierror:            0
 PE:            1 , after reduce, ierror:            0
 PE:            1 , after reduce, SUM(u):    12000000.0    
 PE:            0 , after reduce, SUM(u):    6000000.00    
 PE:            1 , after reduce, SUM(v):    0.00000000    
 PE:            0 , after reduce, SUM(v):    18407592.0    

PE0“应该”在最后一行显示 18000000.0 作为 SUM(v)。

如果我将代码中的 nz 参数从 6 设置为 5，则运行会产生正确的结果。真正令人困惑的是，它以这种方式运行，在 a) 具有 gfortran 5.3 和 openmpi 的 AWS EC2 实例上返回相同的减少值总和，b) 我的笔记本电脑的 gfortran 5.4 和 mpich，以及 c) 工作站的 gfortran 4.4 和 openmpi。

如果我将数组的类型更改为 DOUBLE PRECISION（以及在 MPI_REDUCE 调用中指定）它可以正常工作，即使对于更大的数组也是如此。如果我使用 REAL4 而不是 REAL，它会产生同样糟糕的结果。

我知道这必须很简单，而且我在这里是一个真正的白痴，但我就是不明白这一点。我读过一些建议，我的缓冲区大小需要是小于 2^31-1 的整数值，但这里肯定是这种情况。

Answer 1

这与 MPI 无关，只是求和精度问题：

PROGRAM TEST
    IMPLICIT NONE
    INTEGER, PARAMETER :: nx=1000, ny=1000, nz=6
    REAL, DIMENSION(nx, ny, nz) :: u
    u = 3
    PRINT *, SUM(u)
END PROGRAM test

返回相同的结果。如果您将一个大数字添加到一个小数字中，可能会出现舍入问题，在许多小数字的总和中，这种影响可能会累积成一个重大错误。有求和算法可以防止这种影响，比如Kahan summation，显然Fortran的SUM不是这样实现的。