Я обнаружил некоторые NaN в массиве данных в моем коде на Фортране и, возможно, изолировал проблему в цикле DoMP OpenMP. Когда этот цикл выполняется, NaN приводит к sp% ptl (i)% ph (6):
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
if(sp%ptl(i)%ph(3) >= 2pi .or. sp%ptl(i)%ph(3)< 0D0 ) then
sp%ptl(i)%ph(3) = modulo(sp%ptl(i)%ph(3),2pi)
endif
enddo
enddo
Но если я запускаю то же самое, но добавляю строку, добавляющую 0D0 к фиктивной переменной, значения NaN уходят в sp% ptl (i)% ph (6):
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
if(sp%ptl(i)%ph(3) >= 2pi .or. sp%ptl(i)%ph(3)< 0D0 ) then
sp%ptl(i)%ph(3) = modulo(sp%ptl(i)%ph(3),2pi)
endif
tmp = tmp + 0D0
enddo
enddo
Конечно, в реальном коде есть что-то большее, и это не минимальный рабочий пример. Мой вопрос, однако, почему добавление какой-либо строки в цикл do приводит к тому, что sp% ptl (i)% ph (6) никогда не получает NaN? Является ли плохой идеей наличие в операторе OpenMP только оператора if? Сейчас это сбивает с толку вопрос, как это работает.
ОБНОВЛЕНИЕ Вот минимальный пример, который пока не совсем работает как код большего размера, так как он не имеет NaN, а имеет довольно большие числа в случайных точках массива ptl, но, по крайней мере, показывает основной рабочий процесс. Я компилирую так же, как с большей кодовой базой (компилятор Intel, 18.0.1.163), затем запускаю это с помощью srun -n1 -c24.
Я провел еще несколько тестов с большей кодовой базой и обнаружил, что перекомпиляция подпрограмм, представленных здесь как «mymod» без оптимизации (то есть -O0 -g -C), заставляет NaN уходить.
ОБНОВЛЕНИЕ 2 Не говоря уже о больших числах, я просто забыл инициализацию ptl (i)% ph (теперь добавляется в начале push), теперь, когда он добавлен, я никогда не получаю NaN или большое числа в этом минимальном примере (все еще там в большем коде с оптимизацией на).
module mymod
integer, parameter :: ptl_nphase=8
integer, parameter :: num=1000
integer, parameter :: sml_nthreads=24
type ptl_type
real(8) :: ph(ptl_nphase)
end type ptl_type
contains
logical function is_nan(a)
implicit none
real (8) :: a
is_nan = .not. ( a > 1D0 .or. a < 2D0 )
end function is_nan
subroutine split_indices(total,num_pieces,ibeg,iend)
implicit none
integer :: total
integer :: num_pieces
integer :: ibeg(num_pieces), iend(num_pieces)
integer :: itmp1, itmp2, ioffset, i
if (num_pieces > 0) then
itmp1 = total/num_pieces
itmp2 = mod(total,num_pieces)
ioffset = 0
do i=1,itmp2
ibeg(i) = ioffset + 1
iend(i) = ioffset + (itmp1+1)
ioffset = iend(i)
enddo
do i=itmp2+1,num_pieces
ibeg(i) = ioffset + 1
if (ibeg(i) > total) then
iend(i) = ibeg(i) - 1
else
iend(i) = ioffset + itmp1
ioffset = iend(i)
endif
enddo
endif
end subroutine split_indices
subroutine calc_source(ptl,icycle)
implicit none
type(ptl_type) :: ptl(num)
integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
integer :: icycle
call split_indices(num, sml_nthreads, i_beg, i_end)
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
ptl(i)%ph(6) = ptl(i)%ph(6) + 1D0
enddo
enddo
if (icycle==1) then
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
ptl(i)%ph(7) = ptl(i)%ph(6)
enddo
enddo
endif
end subroutine calc_source
subroutine push1(ptl)
implicit none
type(ptl_type) :: ptl(num)
integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
real(8) :: arr1(5)
call split_indices(num, sml_nthreads, i_beg, i_end)
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
call random_number(arr1)
ptl(i)%ph(1:5) = ptl(i)%ph(1:5) + arr1
enddo
enddo
end subroutine push1
subroutine push(ptl)
implicit none
type(ptl_type) :: ptl(num)
integer :: icycle
integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
call split_indices(num, sml_nthreads, i_beg, i_end)
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
ptl(i)%ph(:) = 0D0
enddo
enddo
do icycle=1,100
call calc_source(ptl,icycle)
call push1(ptl)
call split_indices(num, sml_nthreads, i_beg, i_end)
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
ptl(i)%ph(3) = modulo(ptl(i)%ph(3),6.28)
enddo
enddo
enddo
end subroutine push
end module mymod
program main
use mymod
implicit none
type(ptl_type) :: ptl(num)
integer :: ith, i, i_beg(sml_nthreads), i_end(sml_nthreads)
call push(ptl)
!check for nan
call split_indices(num, sml_nthreads, i_beg, i_end)
!$OMP PARALLEL DO &
!$OMP PRIVATE( ITH, I )
do ith=1, sml_nthreads
do i=i_beg(ith), i_end(ith)
! if (is_nan(ptl(i)%ph(6))) then
!print *,'is_nan',i
print *,ptl(i)%ph(6)
! endif
enddo
enddo
end program main