Как сохранить точность для программы Fortran MPI портативным способом?

У меня есть программа 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 всегда будут иметь одинаковую точность, независимо от того, компилятор и машина?

UPDATE:

Я создал простую программу, которая демонстрирует проблему, которую я вижу, когда мои внутренние значения имеют более высокую точность, чем то, что предоставляется типом 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

Если я создам и запускаю с двумя ядрами, я получаю:

       0   1.12345683574676513672      
       1   4.71241976735884452383E-3998

Теперь измените значение 16 на 15 в selected_real_kind, и я получаю:

       0   1.1234568357467651     
       1   1.1234568357467651

Всегда ли безопасно использовать selected_real_kind(15) с MPI_DOUBLE_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 и планируется установить в MPICH на 3.1.1. См. Билеты OpenMPI 3432 и 3435, а также как билеты MPICH 1769 и 1770.

Ответ 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 с этого момента.