¿Cómo conservo la precisión para un programa Fortran MPI de manera portátil?

precision

Question

Tengo un programa Fortran donde especifico el kind de tipos de datos numéricos en un intento de mantener un nivel mínimo de precisión, independientemente de qué compilador se use para construir el programa. Por ejemplo:

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

El problema es que he usado MPI para paralelizar el código y necesito asegurarme de que las comunicaciones MPI especifiquen el mismo tipo con la misma precisión. Estaba usando el siguiente enfoque para mantener la coherencia con el enfoque de mi programa:

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)

Sin embargo, he descubierto que esta rutina de MPI no es particularmente compatible con las diferentes implementaciones de MPI, por lo que en realidad hace que mi programa no sea portátil. Si MPI_Type_create rutina MPI_Type_create , entonces me queda confiar en los tipos de datos MPI_REAL y MPI_DOUBLE_PRECISION estándar, pero ¿qué MPI_DOUBLE_PRECISION si ese tipo no es coherente con lo que selected_real_kind elige como el tipo real que finalmente pasará MPI? ¿Estoy atascado simplemente usando la declaración real estándar para un tipo de datos, sin ningún atributo de kind y, si lo hago, tengo garantizado que MPI_REAL y real siempre tendrán la misma precisión, independientemente del compilador y la máquina?

ACTUALIZAR:

Creé un programa simple que demuestra el problema que veo cuando mis reales internos tienen una precisión mayor que la que ofrece el tipo 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

Si construyo y corro con 2 núcleos, obtengo:

 0 1.12345683574676513672
1 4.71241976735884452383E-3998

Ahora cambie el 16 a un 15 en selected_real_kind y obtengo:

 0 1.1234568357467651
1 1.1234568357467651

¿Siempre será seguro usar selected_real_kind(15) con MPI_DOUBLE_PRECISION sin importar qué máquina / compilador se use para construir?

Answer 1

Qué tal si:

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

y use MPIREAL en lugar de MPI_DOUBLE_PRECISION a partir de ese momento.

Answer 2

Realmente no es una respuesta, pero tenemos el mismo problema y utilizamos algo como esto:

!> 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)

Y luego tenemos una rutina de inicialización donde hacemos:

!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

Si bien esto no es agradable, funciona razonablemente bien y parece ser utilizable en Cray, IBM BlueGene y clusters Linux convencionales. Lo mejor que puede hacer es presionar a los sitios y proveedores para que lo admitan correctamente en MPI. Por lo que yo sé, se ha corregido en OpenMPI y estaba previsto que se arreglara en MPICH por 3.1.1. Vea las entradas OpenMPI 3432 y 3435 , así como las entradas MPICH 1769 y 1770 .

Answer 3

Utilice el STORAGE_SIZE intrínseco de Fortran 2008 para determinar el número de bytes que requiere cada número y enviarlos como bytes. Tenga en cuenta que STORAGE_SIZE devuelve el tamaño en bits, por lo que deberá dividir por 8 para obtener el tamaño en bytes.

Esta solución funciona para mover datos pero no ayuda a usar reducciones. Para eso, tendrá que implementar una operación de reducción definida por el usuario. Si eso es importante para ti, actualizaré mi respuesta con los detalles.

Por ejemplo:

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

Confirmé que este cambio corrige el problema y la salida que veo es:

 0 1.12345683574676513672
1 1.12345683574676513672