如何以便携方式保留Fortran MPI程序的精度？

Question

我有一个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

现在将16改为15 in selected_real_kind我得到:

       0   1.1234568357467651     
       1   1.1234568357467651  

难道总是要安全使用selected_real_kind(15)与MPI_DOUBLE_PRECISION不管什么机器/编译器是用来做构建？

Answer 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  

Answer 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。