如何以可移植的方式保持 Fortran MPI 程序的精度？答案

【问题标题】：How do I retain precision for a Fortran MPI program in a portable way?如何以可移植的方式保持 Fortran MPI 程序的精度？
【发布时间】：2013-11-12 15:20:04
【问题描述】：

我有一个 Fortran 程序，我在其中指定了数字数据类型的 kind，以尝试保持最低级别的精度，而不管使用什么编译器来构建程序。例如：

integer, parameter :: rsp = selected_real_kind(4)
...
real(kind=rsp) :: real_var

问题是我使用 MPI 来并行化代码，我需要确保 MPI 通信指定具有相同精度的相同类型。我使用以下方法与我的程序中的方法保持一致：

call MPI_Type_create_f90_real(4,MPI_UNDEFINED,rsp_mpi,mpi_err)
...
call MPI_Send(real_var,1,rsp_mpi,dest,tag,MPI_COMM_WORLD,err)

但是，我发现这个 MPI 例程对不同的 MPI 实现并不是特别好的支持，所以它实际上使我的程序不可移植。如果我省略了MPI_Type_create 例程，那么我只能依赖标准的MPI_REAL 和MPI_DOUBLE_PRECISION 数据类型，但是如果该类型与selected_real_kind 选择的最终类型不一致怎么办？被 MPI 传递？我是否坚持只对数据类型使用标准的real 声明，没有kind 属性，如果我这样做，我是否保证MPI_REAL 和real 总是具有相同的精度，不管编译器和机器？

更新：

我创建了一个简单的程序来演示当我的内部实数比 MPI_DOUBLE_PRECISION 类型提供的精度更高时看到的问题：

program main

   use mpi

   implicit none

   integer, parameter :: rsp = selected_real_kind(16)
   integer :: err
   integer :: rank

   real(rsp) :: real_var

   call MPI_Init(err)
   call MPI_Comm_rank(MPI_COMM_WORLD,rank,err)

   if (rank.eq.0) then
      real_var = 1.123456789012345
      call MPI_Send(real_var,1,MPI_DOUBLE_PRECISION,1,5,MPI_COMM_WORLD,err)
   else
      call MPI_Recv(real_var,1,MPI_DOUBLE_PRECISION,0,5,MPI_COMM_WORLD,&
         MPI_STATUS_IGNORE,err)
   end if

   print *, rank, real_var

   call MPI_Finalize(err)

end program main

如果我使用 2 个内核构建和运行，我会得到：

       0   1.12345683574676513672      
       1   4.71241976735884452383E-3998

现在将 selected_real_kind 中的 16 更改为 15，我得到：

       0   1.1234568357467651     
       1   1.1234568357467651

无论使用什么机器/编译器进行构建，使用selected_real_kind(15) 和MPI_DOUBLE_PRECISION 总是安全的吗？

【问题讨论】：

我使用的代码将其工作精度定义为wp = selected_real_kind(14,40)，然后在MPI 调用中使用MPI_DOUBLE_PRECISION。据我所知，没有问题。
我添加了一个简单的例子来说明我需要避免什么问题。所以我想知道，如果我明确定义他们的kind，我的reals 是否总是与MPI_DOUBLE_PRECISION 兼容？
上次我检查过，使用 16 表示四倍精度（因此是 -3998 范围），这不适用于 MPI_DOUBLE_PRECISION。只要您使用 7 到 15（含）之间，您应该可以使用 MPI_DOUBLE_PRECISION。
是的，你是对的。但是MPI_DOUBLE_PRECISION 是否总是意味着无论系统/编译器如何，MPI 通信都保持 7 到 15 位精度？换句话说，如果我将我的reals 指定为 14 位精度，然后与 MPI_DOUBLE_PRECISION 进行 MPI 通信，我是否会担心得到像我发布的那样的结果，其中号码在接收端作为垃圾出现？
不使用selected_real_kind，您可以使用use ISO_FORTRAN_ENV 并使用类型real32 和real64（位数）可移植地指定单精度和双精度实数。

标签： fortran mpi precision

【解决方案1】：

使用 Fortran 2008 内在函数 STORAGE_SIZE 确定每个数字所需的数字字节数并以字节形式发送。请注意，STORAGE_SIZE 以位为单位返回大小，因此您需要除以 8 才能获得以字节为单位的大小。

此解决方案适用于移动数据，但不能帮助您使用归约。为此，您必须实现用户定义的归约操作。如果这对您很重要，我会用详细信息更新我的答案。

例如：

program main

   use mpi

   implicit none

   integer, parameter :: rsp = selected_real_kind(16)
   integer :: err
   integer :: rank

   real(rsp) :: real_var

   call MPI_Init(err)
   call MPI_Comm_rank(MPI_COMM_WORLD,rank,err)

   if (rank.eq.0) then
      real_var = 1.123456789012345
      call MPI_Send(real_var,storage_size(real_var)/8,MPI_BYTE,1,5,MPI_COMM_WORLD,err)
   else
      call MPI_Recv(real_var,storage_size(real_var)/8,MPI_BYTE,0,5,MPI_COMM_WORLD,&
         MPI_STATUS_IGNORE,err)
   end if

   print *, rank, real_var

   call MPI_Finalize(err)

end program main

我确认此更改解决了问题，我看到的输出是：

   0   1.12345683574676513672      
   1   1.12345683574676513672

【讨论】：

【解决方案2】：

不是一个真正的答案，但我们有同样的问题并使用这样的东西：

!> Number of digits for single precision numbers
integer, parameter, public :: single_prec = 6
!> Number of digits for double precision numbers
integer, parameter, public :: double_prec = 15
!> Number of digits for extended double precision numbers
integer, parameter, public :: xdble_prec = 18
!> Number of digits for quadruple precision numbers
integer, parameter, public :: quad_prec = 33

integer, parameter, public :: rk_prec = double_prec

!> The kind to select for default reals
integer, parameter, public :: rk = selected_real_kind(rk_prec)

然后有一个我们做的初始化例程：

!call mpi_type_create_f90_real(rk_prec, MPI_UNDEFINED, rk_mpi, iError)
!call mpi_type_create_f90_integer(long_prec, long_k_mpi, iError)
! Workaround shitty MPI-Implementations.
select case(rk_prec)
case(single_prec)
  rk_mpi = MPI_REAL
case(double_prec)
  rk_mpi = MPI_DOUBLE_PRECISION
case(quad_prec)
  rk_mpi = MPI_REAL16
case default
  write(*,*) 'unknown real type specified for mpi_type creation'
end select
long_k_mpi = MPI_INTEGER8

虽然这不是很好，但它工作得相当好，并且似乎可以在 Cray、IBM BlueGene 和传统的 Linux 集群上使用。最好的办法是推动站点和供应商在 MPI 中正确支持这一点。据我所知，它已在 OpenMPI 中修复，并计划在 3.1.1 之前在 MPICH 中修复。请参阅 OpenMPI 票证 3432 和 3435 以及 MPICH 票证 1769 和 1770。

【讨论】：

你说得对，我的代码使用 OpenMPI 1.4.3 运行良好，但使用 MPICH 时出现问题。我想我最终会回到使用默认的 MPI Fortran 数据类型，而不是尝试手动设置精度。

【解决方案3】：

怎么样：

integer, parameter :: DOUBLE_PREC = kind(0.0d0)
integer, parameter :: SINGLE_PREC = kind(0.0e0)

integer, parameter :: MYREAL = DOUBLE_PREC


if (MYREAL .eq. DOUBLE_PREC) then
   MPIREAL = MPI_DOUBLE_PRECISION
else if (MYREAL .eq. SINGLE_PREC) then
   MPIREAL = MPI_REAL
else
   print *, "Erorr: Can't figure out MPI precision."
   STOP
end if

从那时起使用 MPIREAL 而不是 MPI_DOUBLE_PRECISION。

【讨论】：